PRO reducaln,yr,mo,day,cam,R_s_per_DN $ ,verbose=vb,notv=notv,nointeractive=nointeractive $ ,Inspect = inspect ; find calibration images and dark images for given day ; read calibration images and dark images ; perform dark subtraction ; read flat-field images and do flat-field corrections ; determine center and approximate radius of PLBS images ; average camera signal over central 90% of images ; read in filter throughput values ; calculate Rayleigh-seconds/DN for various channels, using all ; available images for each channel, and adding 0.045 s to nominal ; exposure times ; return array of Rayleigh-seconds/DN values COMMON ffstuff,ffimg IF N_PARAMS() LT 5 THEN BEGIN PRINT,'Usage: reducals,yr,mo,day,cam,R_s_per_DN $' PRINT,' ,verbose=vb,notv=notv,nointeractive=nointeractive' RETURN ENDIF ; Initialize OS-specific items @isitdos ; Load flat-field images that will be needed later getffimg,yr,mo,day ; find calibration images and dark images for given day rdilist,cam,yr,mo,day,0,0,0,n,t,f IF n EQ 0 THEN RETURN cal=WHERE(STRPOS(f,calname) NE -1,ncal) dk=WHERE(STRPOS(f,dkname) NE -1,ndk) IF ncal EQ 0 THEN RETURN ; first get header bytes and identify the real calibration images seq=BYTARR(ncal) FOR j=0,ncal-1 DO BEGIN rdkihd,f(cal(j)),hb IF N_Elements(hb) LT 512 THEN seq(j)=0 ELSE BEGIN h=gethd(hb) seq(j)=(h.misc.exp_scale*h.exp.exposure GT 25) ENDELSE ENDFOR IF TOTAL(seq) EQ 0 THEN RETURN cal=cal(WHERE(seq)) ; calibration images take longer than 25 ms ncal=N_ELEMENTS(cal) ; display calibration image file names PRINT,'Camera '+STRTRIM(cam,2)+' calibration files:' FOR j=0,ncal-1 DO PRINT,f(cal(j)) ; not enough files for a real calibration sequence IF ncal LT 8 THEN RETURN ; Find out exposure times used caltimes=Fix(exposure(f(cal))) dktimes=Fix(exposure(f(dk))) exptimes=caltimes(Uniq(caltimes,Sort(caltimes))) n_exptimes=N_Elements(exptimes) ; read calibration images and dark images ; Now load the dark frames nearest in time to the calibration images firstcal=MIN(t(cal),MAX=lastcal) ; times of first and last cal images dkfr=INTARR(n_exptimes*256,512) ; array to hold dark frames ccdt=FLTARR(n_exptimes,2) ; array to hold CCD temps ;***************************** FOR j=0,n_exptimes-1 DO BEGIN ; find last 1-s dark frame taken before first calibration image bef=Max(Where((t(dk) LT firstcal) AND (dktimes EQ exptimes(j)))) IF bef NE -1 THEN BEGIN ; read in the header bytes and image rdkimg,f(dk(bef)),hb,dkim ; decode the header bytes h=gethd(hb) ; and store the CCD temperature for later reference ccdt(j,0)=h.misc.temp.tcf.mean/10. ENDIF ELSE BEGIN rdmeandk,cam,exptimes(j),dkim,/nointeractive Message,StrTrim(exptimes(j),2)+'-s dark frame synthesized',/INFORMATIONAL ENDELSE ; store the image for later use dkfr(256*j:256*j+255,0:255)=stdsize(Median(dkim,3)) ; now repeat the process for the first 1-s dark frame after calibration aft=Min(Where((t(dk) GT lastcal) AND (dktimes EQ exptimes(j)))) IF aft NE -1 THEN BEGIN rdkimg,f(dk(aft)),hb,dkim h=gethd(hb) ccdt(j,1)=h.misc.temp.tcf.mean/10. ENDIF ELSE BEGIN rdmeandk,cam,exptimes(j),dkim,/nointeractive Message,StrTrim(exptimes(j),2)+'-s dark frame synthesized',/INFORMATIONAL ENDELSE dkfr(256*j:256*j+255,256:511)=stdsize(Median(dkim,3)) ENDFOR ;***************************** ; Before actually processing, load the path names for the throughput ; data tppath=STRARR(2,5) tppath(0,0)=tproot+dd+'829c' tppath(0,1)=tproot+dd+'609c' tppath(0,2)=tproot+dd+'630c' tppath(0,3)=tproot+dd+'668c' tppath(0,4)=tproot+dd+'520c' tppath(1,0)=tproot+dd+'834c' tppath(1,1)=tproot+dd+'821c' tppath(1,2)=tproot+dd+'785c' tppath(1,3)=tproot+dd+'558c' tppath(1,4)=tproot+dd+'714c' ; Now go through the real calibration data files one at a time wtsum=FLTARR(5) sumwt=wtsum R_s_per_DN=wtsum alltfw=FLTARR(ncal) FOR j=0,ncal-1 DO BEGIN rdkimg,f(cal(j)),hb,im h=gethd(hb) null_image=MAX(im) EQ 0 texp=h.misc.exp_scale*h.exp.exposure/1000 j2=Where(exptimes EQ texp) & j2=j2(0) ; CASE texp OF ; 1: j2=0 ; 10: j2=1 ; 60: j2=2 ; ENDCASE alltfw(j)=h.misc.temp.tfw.mean/10. tcf=h.misc.temp.tcf.mean/10. IF ABS(ccdt(j2,0)-tcf) LT ABS(ccdt(j2,1)-tcf) THEN j1=0 ELSE j1=1 ; Convert image to standard size and attempt to correct ; binning-induced overflow im=stdsize(im) wsat=Where(im GE 16383,nsat) ; IF nsat GT 0 THEN Hlp,im(wsat) IF (NOT null_image) AND (nsat LT 3000) THEN BEGIN ; subtract appropriate dark frame im=(im-dkfr(256*j2:256*(j2+1)-1,256*j1:256*(j1+1)-1))>0 ; load and apply appropriate flat-field correction image cam=h.misc.cam filt=h.misc.fw col1=filt*256L row1=cam*256L ffim=ffimg(col1:col1+255,row1:row1+255) im=FIX(ROUND(im*FLOAT(ffim)/100.)) ; determine position and size of PLBS image meancent,im,clbs,rlbs,npix,radius=radlbs,notv=notv ; IF Keyword_Set(vb) THEN Help,radlbs,!Pi*radlbs*radlbs ; IF the PLBS radius is reasonable, carry on IF radlbs GT 10 THEN BEGIN ; restrict averaging to central 90% of PLBS image rad90=SQRT(0.9)*radlbs one2two,LINDGEN(65536L),im,cmask,rmask mk,mask mask(WHERE(SQRT((cmask-clbs)^2+(rmask-rlbs)^2) LE rad90))=1 avgdn=TOTAL(mask*im)/TOTAL(mask) ; If average DN is zero, report the problem IF avgdn GT 1 THEN BEGIN ; read in throughput data for this filter tpfile=getfile(tppath(cam,filt)+dd+'thruput') ndata=N_ELEMENTS(tpfile)-1 tptemp=FLTARR(ndata) tpkr=tptemp FOR k=1,ndata DO BEGIN tptemp(k-1)=FLOAT(getwrd(tpfile(k),0)) tpkr(k-1)=FLOAT(getwrd(tpfile(k),1)) ENDFOR tfw=h.misc.temp.tfw.mean/10. IF tfw EQ 0 THEN BEGIN wtfw=WHERE(alltfw,nwtfw) IF nwtfw GT 0 THEN tfw=alltfw(MAX(wtfw)) ELSE BEGIN tfw=15. PRINT,'Nominal filter wheel temperature of 15 C assumed!' ENDELSE ENDIF tp=spline(tptemp,tpkr,tfw) ; Now calculate sensitivity sens=1000.*tp*(texp+0.045)/avgdn sat_frac=nsat/4150. wtsum(filt)=wtsum(filt)+(texp+0.045)*sens*(1.-sat_frac) sumwt(filt)=sumwt(filt)+(texp+0.045)*(1.-sat_frac) IF KEYWORD_SET(vb) THEN PRINT,cam,filt,texp,tp,avgdn $ ,FORMAT='("cam ",I1," filt ",I1,F6.1' $ +'," s ",F7.2," kR ",F7.1," DN")' IF KEYWORD_SET(vb) THEN PRINT,sens,FORMAT='(F6.1," Rayleigh-s/DN")' ENDIF ELSE BEGIN MESSAGE,STRING(f(cal(j)) $ ,FORMAT='(A,": average PLBS DN < 1; skip")') $ ,/INFORMATIONAL ENDELSE ENDIF ELSE MESSAGE,f(cal(j))+': image model invalid' $ ,/INFORMATIONAL ENDIF ELSE BEGIN IF null_image THEN MESSAGE,f(cal(j))+' signal too low: ' $ +STRTRIM(h.misc.pixel_max,2)+'; skip' $ ,/INFORMATIONAL IF KEYWORD_SET(vb) AND (nsat GE 3000) THEN BEGIN MESSAGE,STRING(f(cal(j)),nsat $ ,FORMAT='(/,A,": ",I0.0," saturated pixels")') $ ,/INFORMATIONAL ; Window,/Free,XSize=256,YSize=256 ; mim=im ; mim(wsat)=2*mim(wsat) ; TVScl,mim ; Wait,5 ; Stop ; WDelete,!D.Window ENDIF ENDELSE ENDFOR FOR j=0,4 DO IF sumwt(j) GT 0. THEN R_s_per_DN(j)=wtsum(j)/sumwt(j) ; Save results to a file if they mean anything IF N_ELEMENTS(WHERE(R_s_per_DN)) GE 2 THEN BEGIN ; First see whether results are reasonable, i.e., within 25% ; of other results nzR=WHERE(R_s_per_DN) IF cam EQ 0 THEN frac=R_s_per_DN/[128.8,61.6,71.3,61.9,90.0] $ ELSE frac=R_s_per_DN/[112.1,102.8,86.7,68.4,53.9] ok=(TOTAL(ABS(frac-1.) LT 0.25) EQ N_ELEMENTS(nzR)) savename=STRING(calsroot,dd,yr MOD 100,mo,day,dc,cam,cal_suffix $ ,FORMAT='(A,A1,3I2.2,A1,I1,A)') IF NOT ok THEN BEGIN ; tell the user and change the filename PRINT,yr MOD 100,mo,day,cam $ ,FORMAT='("Wonky calibration data on ",3I2.2," for camera ",I1)' savename=savename+'w' ENDIF OPENW,unit,savename,/GET_LUN PRINTF,unit,R_s_per_DN,FORMAT='(5F6.1)' FREE_LUN,unit ENDIF RETURN END