PRO merslice,img,azimg,zaimg,meraz,elev,slice,sd,swidth=swidth,debug=debug IF NOT KEYWORD_SET(swidth) THEN swidth=1. ; default to a 1-deg slice ximg=SIN(azimg*!DTOR)*SIN(zaimg*!DTOR) ; X points E, yimg=COS(azimg*!DTOR)*SIN(zaimg*!DTOR) ; Y points N, zimg=COS(zaimg*!DTOR) ; and Z points up. ; ; Take dot product of (X,Y,Z) with vector 90 deg CW from azimuth ; tmp=ximg*COS(meraz*!DTOR)-yimg*SIN(meraz*!DTOR) ; ; Arccosine gives angles between image pixels and vector 90 deg CW ; merdist=!RADEG*ACOS(tmp) ; ; Now find pixels within swidth/2 of 90 deg from vector 90 deg CW ; These, if any, are the pixels in the desired meridian slice ; merslice=WHERE(ABS(merdist-90.) LT swidth/2.,nmer) IF nmer GT 0 THEN BEGIN ; ; Now find angles w.r.t. horizontal vector pointing to azimuth ; These are really 'elevation angles' lying in [0..180] deg ; cosel=ximg*SIN(meraz*!DTOR)+yimg*COS(meraz*!DTOR) sinel=zimg alpha=!RADEG*ATAN(sinel,cosel) ; ; Now select pixels in meridian slice ; slice=img(merslice) IF KEYWORD_SET(debug) THEN BEGIN imax=pcentile(img,.99) imin=pcentile(img,.01) mk,mask mask(merslice)=1 timg=(img(-0.5) hmax=(ROUND(elmax)+0.5)<179.5 IF KEYWORD_SET(debug) THEN PRINT,elmin,elmax,hmin,hmax,FORMAT='(4F7.2)' h=HISTOGRAM(elev,MIN=hmin,MAX=hmax,REVERSE_INDICES=r) eslice=FINDGEN(181) mslice=FLTARR(181) sdslice=mslice FOR e=hmin,hmax-1 DO BEGIN i=FIX(e-hmin) IF (r(i) LT (r(i+1)-1)) THEN BEGIN epop=slice(r(r(i):r(i+1)-1)) mslice(CEIL(e))=mean(epop) IF N_ELEMENTS(epop) GT 1 THEN sdslice(CEIL(e))=stdev(epop) ENDIF ENDFOR elev=eslice slice=mslice sd=sdslice IF !D.WINDOW EQ 2 THEN BEGIN PRINT,'Press any key to continue' k=GET_KBRD(1) WDELETE,2 ENDIF ENDIF ELSE BEGIN MESSAGE,'No pixels found on selected meridian!',/INFORMATIONAL ENDELSE RETURN END