; PBINCAM.PRO
;
; DPS, August 1992
;
; Program to test the in-camera passbands of the Polar Camera filters.
;
; Modified 1993/05/11 from TELTEST.PRO.
;
; Extract mean camera signal from a sequence of images taken
; viewing a diffusing screen monochromatically illuminated at
; wavelengths throughout the passband of a filter.  A sequence of
; images obtained viewing the screen at the center 
; of the field of view are processed, and the wavelength variation
; of the screen brightnesses in the sequence is plotted to
; show the filter passband.
; 
PRO pbincam,ncw,cw,rpb,hw,tw,noplot=no

; First, find out name of ball park
ncw=STRING(ncw,FORMAT='(I3,"c")')
path=STRING(ncw,FORMAT='("/jasper/cnsr3_data1/passband/",A,"/")')

; Look for a 'signal.za' file,  If found, we can save a bunch of
; time!
sig=FINDFILE(path+'signal.za',COUNT=nsig)
fromimage=(nsig EQ 0)
IF fromimage THEN BEGIN
    PRINT,'To use this program you need the following files:'
    PRINT,'--------------------------------------------------------------------------------'
    PRINT,'wl.nnn	(where nnn is the nominal center wavelength in nm) containing the'
    PRINT,'	wavelengths at which images were acquired'
    PRINT,'lite	containing (in order) the names of the image files acquired at the'
    PRINT,'	wavelengths specified in wl.nnn'
    PRINT,'mask	containing a list of the image pixels representing the peak screen'
    PRINT,'	brightness in the images'
    PRINT,''
    PRINT,'If you don''t have these files in place, quit now and generate them!'

; read file of wavelengths
    wlfile=path+STRING(ncw,FORMAT='("wl.",A)')
    PRINT,wlfile
    wavelength=627.0 & filename=''

; the following is a standard pattern for reading in data from files
    OPENR,unit,wlfile,/GET_LUN
    READF,unit,nw		; first read number of wavelength entries
    w=FLTARR(nw) 		; create an array to hold the wavelengths
    FOR i=0,nw-1 DO BEGIN
        READF,unit,wavelength	; read the wavelength, and ...
        w(i)=wavelength		; store in in the array, since we can't
    ENDFOR			; read directly into an array element
    CLOSE,unit & FREE_LUN,unit

; read file of image file names
    OPENR,unit,path+'lite',/GET_LUN
    READF,unit,n
    IF (n NE nw) THEN BEGIN	; check that there are the same number of
        PRINT,STRING(n,nw $	; images as there are wavelengths
       ,FORMAT='(I2," images doesn''t match ",I2," wavelengths!")')
        PRINT,'Aborting...'
	GOTO,ENDPOINT
    ENDIF
    f=STRARR(n)
    FOR i=0,n-1 DO BEGIN
        READF,unit,filename
        f(i)=STRCOMPRESS(filename,/REMOVE_ALL)	; lose all white space
    ENDFOR
    CLOSE,unit & FREE_LUN,unit

; read file of pixels showing illuminated screen
    OPENR,1,path+'mask' & READF,1,nm & wm=LONARR(nm)
    READF,1,wm & CLOSE,1

; Determine minimum pixel values in 'illuminated' region over whole
; set of images; subtract these rather than dark frame
    illim=INTARR(N_ELEMENTS(wm),n)
    PRINT,'Generating minimum image from '+STRTRIM(FIX(n),2)+' images...'
    rdkimg,f(0),h,im
    PRINT,0,FORMAT='(I3,$)'
    mim=im
    illim(*,0)=im(wm)	; save important pixels for later use
    FOR i=1,n-1 DO BEGIN
        IF i MOD (75/3) EQ 0 THEN PRINT,''
        PRINT,i,FORMAT='(I3,$)'
        rdkimg,f(i),h,im
        wl=WHERE(im LT mim,nwl)
        IF nwl GE 1 THEN mim(wl)=im(wl)	; adjust mim if necessary
        illim(*,i)=im(wm)	; save important pixels as above
    ENDFOR
    PRINT,''

; calculate mean and standard deviation of screen brightnesses versus
; wavelength
    m=FLTARR(n) & sd=m
    FOR i=0,n-1 DO BEGIN
        illim(*,i)=illim(*,i)-mim(wm)
        m(i)=TOTAL(illim(*,i))/N_ELEMENTS(wm)	; mean and ...
        sd(i)=STDEV(illim(*,i))	; standard deviation of pixels on screen
    ENDFOR

; save results
    OPENW,unit,path+'signal.za',/GET_LUN
    out=[TRANSPOSE(w),TRANSPOSE(m),TRANSPOSE(sd)]
    PRINTF,unit,3,nw		; write out # of columns and # of rows
    PRINTF,unit,out		; write out array
    CLOSE,unit
    FREE_LUN,unit

; This is where we get the previously-saved data from 'signal.za'.
ENDIF ELSE rdcols,path+'signal.za',nw,w,m,sd

; Remove points with artificially low signals (like 0)
good=WHERE(m)
w=w(good)
m=m(good)
sd=sd(good)

; Wherever the data came from, it's time to plot them up.
wlo=MIN(w,MAX=whi)
IF NOT KEYWORD_SET(no) THEN BEGIN
    PLOT,w,m,XRANGE=[wlo,whi],YRANGE=[0,MAX([m+sd])],XSTYLE=1 $
        ,XTITLE="Wavelength ( nm )" $
        ,YTITLE="Average Signal from Center of Illuminated Screen ( DN )" $
        ,TITLE=STRING(ncw $
        ,FORMAT='("Passband of ",A," nm Filter in PoCa Optics Module")')
    ERRPLOT,w,(m-sd)>0,m+sd		; plot error bars
    ENDIF

; determine transmission curve widths at half- and tenth-peak transmission
hw=myfwhm(w,m,/whole_curve,yzero=MIN(m),yhm=0.5,l_format='(F5.2)' $
	      ,hp_pts=hpts,line_pts=hht,noplot=no)
tw=myfwhm(w,m,/whole_curve,yzero=MIN(m),yhm=0.1,l_format='(F5.2)' $
	      ,hp_pts=tpts,line_pts=tht,noplot=no)
cw=TOTAL(hpts)/2.	; average half-transmission points to get center wavelength
ratio=hw/tw

; Locate global minima on either side of central maximum
wcmax=WHERE(m EQ MAX(m))
wcmax=wcmax(0)
lomin=MIN(m(0:wcmax-1),wlomin)
himin=MIN(m(wcmax+1:*),whimin)
whimin=whimin+wcmax+1

; Reduce values outside global minima to the minimum values
m(0:wlomin)=MIN([lomin,himin])
m(whimin:*)=MIN([lomin,himin])

; Calculate rectangular passband
w1=FINDGEN(10*(whi-wlo)+1)*0.1+wlo	; wavelengths at uniform 0.1 nm intervals
m1=spline(w,m,w1)			; mean DN values ...
m1n=MIN(m1,MAX=m1x)			; min and max values
rm1=(m1-m1n)/(m1x-m1n)			; relative transmissions
rpb=TOTAL(rm1)*0.1			; rectangular passband (nm)

; Annotate graph, if desired
IF NOT KEYWORD_SET(no) THEN BEGIN
    wrange=whi-wlo		; length of X (wavelength) axis
    top=!Y.CRANGE(1)		; upper limit of Y (camera signal) axis
    XYOUTS,wlo+0.05*wrange,0.9*top,STRING(cw,FORMAT='("Center ",F6.2," nm")')
    XYOUTS,wlo+0.05*wrange,0.85*top,STRING(ratio,FORMAT='("Shape ",F4.2)')
    ENDIF
ENDPOINT:
END