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<imax)-imin
	timg=timg-mask*timg+mask*(imax-imin)
	WINDOW,2,XSIZE=256,YSIZE=256
	TVSCL,timg
	ENDIF
    ;
    ; Select elevation angles corresponding to meridian slice
    ;
    elev=alpha(merslice)
    ;
    ; Now average all measurements over 1-degree bins to smooth
    ; the profile, and optionally supply the standard deviation to
    ; the calling routine
    ;
    elmin=MIN(elev,MAX=elmax)
    hmin=(ROUND(elmin)-0.5)>(-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)
	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)
	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