; TELTEST.PRO
;
; DPS, August 1992
;
; Program to test the telecentricity of the Polar Camera optics.
;
; Extract mean camera signal from two sequences of images taken
; viewing a diffusing screen monochromatically illuminated at
; wavelengths throughout the passband of a filter.  Sequences of
; images obtained viewing the screen at the center and the periphery
; of the field of view are processed, and the wavelength variation
; of the screen brightnesses in the two sequences are overplotted to
; show any shift in the filter passband from center to periphery.
; 
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,'lite0	containing (in order) the names of the image files acquired at the'
PRINT,'	wavelengths specified in wl.nnn, at zenith angle 0 degrees'
PRINT,'lite57	containing (in order) the names of the image files acquired at the'
PRINT,'	wavelengths specified in wl.nnn, at zenith angle 57 degrees'
PRINT,'dark0	containing the names of the dark frames acquired at a zenith angle of'
PRINT,'	0 degrees'
PRINT,'dark57	containing the names of the dark frames acquired at a zenith angle of'
PRINT,'	57 degrees'
PRINT,'mask0	containing a list of the image pixels representing the peak screen'
PRINT,'	brightness in the images acquired at 0 degrees zenith angle'
PRINT,'mask57	containing a list of the image pixels representing the peak screen'
PRINT,'	brightness in the images acquired at 57 degrees zenith angle'
PRINT,''
PRINT,'If you don''t have these files in place, quit now and generate them!'
STOP
;
; Find out name of ball park
;
PRINT,'Enter the nominal center wavelength (nm) of the filter'
READ,ncw
defpath=STRING(ncw,FORMAT='(I3,"/")')	; create default path
; read file of wavelengths, common to both zenith angles
wlfile=defpath+'wl.'+STRING(ncw,FORMAT='(I3)')
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 for 0 deg z.a. images
OPENR,unit,defpath+'lite0',/GET_LUN
READF,unit,n0
IF (n0 NE nw) THEN BEGIN	; check the ther are the same number of
    PRINT,STRING(n0,nw, $	; images as there are wavelengths
     FORMAT='(I2," images at 0 deg doesn''t match ",I2," wavelengths!")')
    STOP
ENDIF
za0f=STRARR(n0)
FOR i=0,n0-1 DO BEGIN
    READF,unit,filename
    za0f(i)=STRCOMPRESS(filename,/REMOVE_ALL)	; lose all white space
ENDFOR
CLOSE,unit & FREE_LUN,unit
; read file of image file names for 57 deg z.a. images
OPENR,unit,defpath+'lite57',/GET_LUN
READF,unit,n57
IF (n57 NE nw) THEN BEGIN
    PRINT,STRING(n57,nw, $
     FORMAT='(I2," images at 57 deg doesn''t match ",I2," wavelengths!")')
    STOP
ENDIF
za57f=STRARR(n57)
FOR i=0,n57-1 DO BEGIN
    READF,unit,filename
    za57f(i)=STRCOMPRESS(filename,/REMOVE_ALL)
ENDFOR
CLOSE,unit & FREE_LUN,unit
; read file of image file names for 0 deg z.a. dark frames
OPENR,unit,defpath+'dark0',/GET_LUN
READF,unit,nd0
za0d=STRARR(nd0) 
FOR i=0,nd0-1 DO BEGIN
    READF,unit,filename
    za0d(i)=STRCOMPRESS(filename,/REMOVE_ALL)
ENDFOR
CLOSE,unit & FREE_LUN,unit
; read file of image file names for 57 deg z.a. dark frames
OPENR,unit,defpath+'dark57',/GET_LUN
READF,unit,nd57
za57d=STRARR(nd57)
FOR i=0,nd57-1 DO BEGIN
    READF,unit,filename
    za57d(i)=STRCOMPRESS(filename,/REMOVE_ALL)
ENDFOR
CLOSE,unit & FREE_LUN,unit
; read file of pixels showing illuminated screen at 0 deg z.a.
OPENR,1,defpath+'mask0' & READF,1,nm0 & wm0=INTARR(nm0)
READF,1,wm0 & CLOSE,1
; read file of pixels showing illuminated screen at 57 deg z.a.
OPENR,1,defpath+'mask57' & READF,1,nm57 & wm57=INTARR(nm57)
READF,1,wm57 & CLOSE,1
; define minimum (i.e., spikefree) dark frame for 0 deg z.a.
FOR i=0,nd0-1 DO BEGIN
    rdkimg,defpath+za0d(i),h,dkf
    IF (i EQ 0) THEN dk0=dkf $	; first dark frame gives baseline
    ELSE BEGIN
        wl=WHERE(dkf LT dk0)	; wherever second frame < first,
        dk0(wl)=dkf(wl)		; adopt value from second
    ENDELSE
ENDFOR
; define minimum (i.e., spikefree) dark frame for 0 deg z.a.
FOR i=0,nd57-1 DO BEGIN
    rdkimg,defpath+za57d(i),h,dkf
    IF (i EQ 0) THEN dk57=dkf $
    ELSE BEGIN
        wl=WHERE(dkf LT dk57)
        dk57(wl)=dkf(wl)
    ENDELSE
ENDFOR
; calculate mean and standard deviation of screen brightnesses versus
; wavelength at 0 deg z.a.
m0=FLTARR(n0) & sd0=m0
FOR i=0,n0-1 DO BEGIN
    rdkimg,defpath+za0f(i),h,im0	; read in image showing screen
    im0=im0-dk0				; subtract dark frame
    m0(i)=TOTAL(im0(wm0))/N_ELEMENTS(wm0)	; mean and ...
    sd0(i)=STDEV(im0(wm0))	; standard deviation of pixels on screen
ENDFOR
; calculate mean and standard deviation of screen brightnesses versus
; wavelength at 57 deg z.a.
m57=FLTARR(n57) & sd57=m57
FOR i=0,n57-1 DO BEGIN
    rdkimg,defpath+za57f(i),h,im57
    im57=im57-dk57
    m57(i)=TOTAL(im57(wm57))/N_ELEMENTS(wm57)
    sd57(i)=STDEV(im57(wm57))
ENDFOR
OPENW,7,defpath+'signal.za'
out=[TRANSPOSE(w),TRANSPOSE(m0),TRANSPOSE(sd0),TRANSPOSE(m57), $
     TRANSPOSE(sd57)]	; concatenate columns into 5 by nw array
PRINTF,7,5,nw		; write out # of columns and # of rows
PRINTF,7,out		; write out array
CLOSE,7
PRINT,STRING(MIN(w,MAX=wmx),wmx, $
      FORMAT='("Data span wavelengths from ",F5.1," nm to ",F5.1," nm;")')
PRINT,'Enter limits for wavelength axis'
READ,wlo,whi		; get user input
!p.title=STRING(ncw, $
         FORMAT='("Passband of ",I3," nm Filter in PoCa Optics Module")')
!x.title="Wavelength ( nm )"
!y.title="Average Signal from Center of Illuminated Screen ( DN )"
PLOT,w,m0,XRANGE=[wlo,whi],YRANGE=[0,MAX([m0+sd0,m57+sd57])],XSTYLE=1
ERRPLOT,w,m0-sd0,m0+sd0		; plot error bars on 0 deg z.a. curve
OPLOT,w,m57,line=2
ERRPLOT,w,m57-sd57,m57+sd57
top=!Y.CRANGE(1)		; upper limit of Y (camera signal) axis
wrange=whi-wlo			; length of X (wavelength) axis
; put descriptive legend on graph
PLOTS,[wlo+0.05*wrange,wlo+0.10*wrange],[0.95*top,0.95*top]
XYOUTS,wlo+0.12*wrange,0.95*top,'0 !9%!x zenith angle'
PLOTS,[wlo+0.05*wrange,wlo+0.1*wrange],[0.9*top,0.9*top],linestyle=2
XYOUTS,wlo+0.12*wrange,0.9*top,'57 !9%!x zenith angle'
; make a list of all 'old' in-optics filter calibration data files
; pbold=FINDFILE('/wesson/user/steele/PAC/filter/calib/*/*.spl',COUNT=nold)
; PRINT,'You can overplot one of the following old passband curves for comparison.'
; FOR i=0,nold-1 DO PRINT,STRING(i,pbold(i),FORMAT='(I1,2X,A)')
; rfilter,wold,sold		; read in the selected data file
; nsold=sold*(MAX(m0)/MAX(sold))	; normalize the old data to the new data
; OPLOT,wold,nsold,PSYM=4		; and overplot it for comparison
; PLOTS,[wlo+0.05*wrange,wlo+0.10*wrange],[0.85*top,0.85*top],PSYM=4
; XYOUTS,wlo+0.12*wrange,0.85*top,'0 !9%!x z. a. (old)'	; update legend
END