PRO smoothff,cam,filt,ffimean,ffisd,smffi,fit

rdffimg,cam,filt,ffimg,az,el
PRINT,'Apply 5x5 median smoothing to the flat-field image first? (Y/[N])'
resp=''
READ,resp
IF STRUPCASE(resp) EQ 'Y' THEN ffi=MEDIAN(ffimg,5) ELSE ffi=ffimg
loadct,13,/silent
WINDOW,0
TVSCL,ffi,0
TVSCL,ffi,0

elmin=0. & elmax=90. & elres=0.25
; Create the array 'r' which contains pointers back to all pixels,
; ordered by increasing elevation.  This will allow us to average
; flat-field image pixel values over all pixels at the same elevation
; angle.
elh=HISTOGRAM(el,MIN=elmin,MAX=elmax,BINSIZE=elres,REVERSE_INDICES=r)

ffimean=FLTARR((elmax-elmin)/elres)
ffisd=ffimean
smffi=INTARR(256,256)
FOR j=0,N_ELEMENTS(elh)-1 DO BEGIN
    ; first collect all the values from the given elevation
    IF elh(j) NE 0 THEN BEGIN
        ffivals=ffi(r(r(j):r(j+1)-1))
        IF N_ELEMENTS(ffivals) GT 0 THEN BEGIN
	    ; if there are some values, take mean and st. dev.
	    IF MAX(ffivals,MIN=fm)-fm GT 0 THEN BEGIN
                ffisd(j)=stdev(ffivals,mm)
                ffimean(j)=mm
		ENDIF ELSE ffimean(j)=fm	; s.d.=0 since all values equal
	    ; then set the appropriate pixels in the smoothed flat-field
	    ; image to the mean value
            smffi(r(r(j):r(j+1)-1))=ROUND(ffimean(j))
	    ENDIF ELSE BEGIN; otherwise set them to 
	    ffisd(j)=-999.	; dummy values,
	    ffimean(j)=-999.; obviously!
	    ENDELSE
        ENDIF
    ENDFOR

TVSCL,smffi,1

; restrict plot to lower half of window
pp=[!X.MARGIN(0)*!D.X_CH_SIZE,!Y.MARGIN(0)*!D.Y_CH_SIZE $
   ,!D.X_SIZE-!X.MARGIN(1)*!D.X_CH_SIZE,0.48*!D.Y_SIZE]
nl=N_ELEMENTS(ffimean)
; construct mean profiles across whole field of view for ...
allels=FINDGEN(2*nl-1)*elres+elmin
; the mean flat-field response, and ...
ffmerid=INDGEN(2*nl-1)
ffmerid(0:nl-1)=ffimean
ffmerid(nl:2*nl-2)=reverse(ffimean(0:nl-2))
; the standard deviation in the above
sdmerid=ffmerid
sdmerid(0:nl-1)=ffisd
sdmerid(nl:2*nl-2)=reverse(ffisd(0:nl-2))
; now plot them versus elevation angle
PLOT,allels,ffmerid,POSITION=pp,/DEVICE,XTITLE='Elevation (deg)' $
    ,YTITLE='Mean pixel value',/NOERASE,XRANGE=[0,180],XSTYLE=1 $
    ,XTICKS=6,XTICKV=INDGEN(7)*30
OPLOT,allels,sdmerid,psym=3
; now set up to fit a sum of cosines to the mean response vs el.
fit=1
cs=1
PRINT,'How many terms to fit? (min. 2, recommend 4)'
READ,nterms
; select only angles where the response is reasonably certain
p=WHERE((sdmerid>2)/(ffmerid>1) LT 0.1)
; in effect, fit versus zenith angle so cosines are symmetric
pfit=svdfit(!DTOR*(allels(p)-90.),ffmerid(p),nterms,weight=1./(sdmerid(p)>2) $
	   ,funct='cosines',yfit=fit,chisq=cs)
; show best-fit curve
OPLOT,allels(p),fit,THICK=2,COLOR=!D.N_COLORS/2
PRINT,"Fit parameters:"
PRINT,pfit
PRINT,"Chi squared = ",STRTRIM(cs,2)

; show difference of actual montage and best smooth fit
WINDOW,2
tmp=smffi
FOR j=MIN(p),MIN(p)+N_ELEMENTS(p)/2 DO $
    tmp(r(r(j):r(j+1)-1))=ROUND(fit(j-MIN(p)))
PRINT,'Now displaying relative error between actual montage and fit'
tvsi,(((ffi-tmp)/(tmp>1))<0.3)>(-0.3),0,0,255,255,window=2 $
    ,contour=FINDGEN(7)*0.05-0.15

PRINT,'Adopt this fit? (Y/[N])'
READ,resp
IF STRUPCASE(resp) EQ 'Y' THEN BEGIN
    smffi=tmp
    outname=STRING(cam,filt,FORMAT='("~cnsr3/flat/ffimg.",2I1,".smooth")')
    OPENW,out,outname,/GET_LUN,/XDR
    ss=SIZE(smffi)
    WRITEU,out,ss(1),ss(2)
    WRITEU,out,smffi
    CLOSE,out
    FREE_LUN,out
    WSET,0
    TVSCL,smffi,1
    ENDIF

RETURN
END

; The following lines are almost straight out of ~steele/idl/image_cont.pro
;CONTOUR,[[0,0],[1,1]],/NODATA, XSTYLE=4, YSTYLE = 4
;px = !X.WINDOW * !D.X_VSIZE     ;Get size of window in device units
;py = !Y.WINDOW * !D.Y_VSIZE
;dpy=py(1)-py(0)			;height in device pixels of plot window
;a=CONGRID((((ffi-tmp)/(tmp>1))<0.5)>(-0.5),dpy,dpy)
;TVSCL,a,px(0),py(0)
;sz = SIZE(a)                    ;Size of image
;six = FLOAT(sz(1))              ;Image sizes
;siy = FLOAT(sz(2))
;swx=six
;swy=siy

;levs=FINDGEN(7)*0.1-0.3		; -0.3, -0.2, ..., 0.2, 0.3
;mx = !D.N_COLORS-1		;Brightest color
;colors = [mx,mx,mx,0,0,0]	;color vectors
;IF !ORDER EQ 1 THEN ar=ROTATE(a,7) ELSE ar=a    ; DPS, Apr, 1993
;CONTOUR,ar,/NOERASE,/XSTYLE,/YSTYLE $	;Do the contour ; 'a' changed to 'ar'
;       ,POSITION=[px(0),py(0),px(0)+swx,py(0)+swy],/DEVICE $
;       ,LEVELS=levs,/FOLLOW;$
;;      ,C_COLOR=colors