; MTELTEST.PRO
;
; DPS, October 1992, based on TELTEST.PRO by 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,'lite55	containing (in order) the names of the image files acquired at the'
PRINT,'	wavelengths specified in wl.nnn, at zenith angle 55 degrees'
PRINT,'lite-55	containing (in order) the names of the image files acquired at the'
PRINT,'	wavelengths specified in wl.nnn, at zenith angle -55 degrees'
PRINT,'dark0	containing the names of the dark frames acquired at a zenith angle of'
PRINT,'	0 degrees'
PRINT,'dark55	containing the names of the dark frames acquired at a zenith angle of'
PRINT,'	55 degrees'
PRINT,'dark-55	containing the names of the dark frames acquired at a zenith angle of'
PRINT,'	-55 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,'mask55	containing a list of the image pixels representing the peak screen'
PRINT,'	brightness in the images acquired at 55 degrees zenith angle'
PRINT,'mask-55	containing a list of the image pixels representing the peak screen'
PRINT,'	brightness in the images acquired at -55 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=FLOAT(ncw) & 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 55 deg z.a. images
OPENR,unit,defpath+'lite55',/GET_LUN
READF,unit,n55
IF (n55 NE nw) THEN BEGIN
    PRINT,STRING(n55,nw, $
     FORMAT='(I2," images at 55 deg doesn''t match ",I2," wavelengths!")')
    STOP
ENDIF
za55f=STRARR(n55)
FOR i=0,n55-1 DO BEGIN
    READF,unit,filename
    za55f(i)=STRCOMPRESS(filename,/REMOVE_ALL)
ENDFOR
CLOSE,unit & FREE_LUN,unit
; read file of image file names for -55 deg z.a. images
OPENR,unit,defpath+'lite-55',/GET_LUN
READF,unit,nn55
IF (nn55 NE nw) THEN BEGIN
    PRINT,STRING(nn55,nw, $
     FORMAT='(I2," images at -55 deg doesn''t match ",I2," wavelengths!")')
    STOP
ENDIF
zan55f=STRARR(nn55)
FOR i=0,nn55-1 DO BEGIN
    READF,unit,filename
    zan55f(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 55 deg z.a. dark frames
OPENR,unit,defpath+'dark55',/GET_LUN
READF,unit,nd55
za55d=STRARR(nd55)
FOR i=0,nd55-1 DO BEGIN
    READF,unit,filename
    za55d(i)=STRCOMPRESS(filename,/REMOVE_ALL)
ENDFOR
CLOSE,unit & FREE_LUN,unit
; read file of image file names for -55 deg z.a. dark frames
OPENR,unit,defpath+'dark-55',/GET_LUN
READF,unit,ndn55
zan55d=STRARR(ndn55)
FOR i=0,ndn55-1 DO BEGIN
    READF,unit,filename
    zan55d(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 55 deg z.a.
OPENR,1,defpath+'mask55' & READF,1,nm55 & wm55=INTARR(nm55)
READF,1,wm55 & CLOSE,1
; read file of pixels showing illuminated screen at -55 deg z.a.
OPENR,1,defpath+'mask-55' & READF,1,nmn55 & wmn55=INTARR(nmn55)
READF,1,wmn55 & 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 55 deg z.a.
FOR i=0,nd55-1 DO BEGIN
    rdkimg,defpath+za55d(i),h,dkf
    IF (i EQ 0) THEN dk55=dkf $
    ELSE BEGIN
        wl=WHERE(dkf LT dk55)
        dk55(wl)=dkf(wl)
    ENDELSE
ENDFOR
; define minimum (i.e., spikefree) dark frame for -55 deg z.a.
FOR i=0,ndn55-1 DO BEGIN
    rdkimg,defpath+zan55d(i),h,dkf
    IF (i EQ 0) THEN dkn55=dkf $
    ELSE BEGIN
        wl=WHERE(dkf LT dkn55)
        dkn55(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 55 deg z.a.
m55=FLTARR(n55) & sd55=m55
FOR i=0,n55-1 DO BEGIN
    rdkimg,defpath+za55f(i),h,im55
    im55=im55-dk55
    m55(i)=TOTAL(im55(wm55))/N_ELEMENTS(wm55)
    sd55(i)=STDEV(im55(wm55))
ENDFOR
; calculate mean and standard deviation of screen brightnesses versus
; wavelength at -55 deg z.a.
mn55=FLTARR(nn55) & sdn55=mn55
FOR i=0,nn55-1 DO BEGIN
    rdkimg,defpath+zan55f(i),h,imn55
    imn55=imn55-dkn55
    mn55(i)=TOTAL(imn55(wmn55))/N_ELEMENTS(wmn55)
    sdn55(i)=STDEV(imn55(wmn55))
ENDFOR
OPENW,7,defpath+'signal.za'
; concatenate columns into 7 by nw array
out=[TRANSPOSE(w),TRANSPOSE(m0),TRANSPOSE(sd0),TRANSPOSE(m55), $
     TRANSPOSE(sd55),TRANSPOSE(mn55),TRANSPOSE(sdn55)]
PRINTF,7,7,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,m55+sd55,mn55+sdn55])],XSTYLE=1
ERRPLOT,w,m0-sd0,m0+sd0		; plot error bars on 0 deg z.a. curve
OPLOT,w,m55,line=2
ERRPLOT,w,m55-sd55,m55+sd55
OPLOT,w,mn55,line=3
ERRPLOT,w,mn55-sdn55,mn55+sdn55
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,'55 !9%!x zenith angle'
PLOTS,[wlo+0.05*wrange,wlo+0.1*wrange],[0.85*top,0.85*top],linestyle=3
XYOUTS,wlo+0.12*wrange,0.85*top,'-55 !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.8*top,0.8*top],PSYM=4
XYOUTS,wlo+0.12*wrange,0.8*top,'0 !9%!x z. a. (old)'	; update legend
END