PRO reducals,yr,mo,day,cam,R_s_per_DN,verbose=vb,notv=notv $
				     ,nointeractive=nointeractive

; 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

; 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
    h=gethd(hb)
    seq(j)=(h.misc.exp_scale*h.exp.exposure GT 25)
    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

; 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(512,768)			; array to hold dark frames
ccdt=FLTARR(2,3)			; array to hold CCD temps

; find last 1-s dark frame taken before first calibration image
bef1s=MAX(WHERE((t(dk) LT firstcal) AND (STRPOS(f(dk),'001.0') NE -1)))
IF bef1s NE -1 THEN BEGIN
    ; read in the header bytes and image
    rdkimg,f(dk(bef1s)),hb,dkim
    ; decode the header bytes
    h=gethd(hb)
    ; and store the CCD temperature for later reference
    ccdt(0,0)=h.misc.temp.tcf.mean/10.
    ENDIF ELSE rdmeandk,cam,1,dkim,nointeractive=nointeractive
; store the image for later use
dkfr(0,0)=stdsize(MEDIAN(dkim,3))

; now repeat the process for the first 1-s dark frame after calibration
aft1s=MIN(WHERE((t(dk) GT lastcal) AND (STRPOS(f(dk),'001.0') NE -1)))
IF aft1s NE -1 THEN BEGIN
    rdkimg,f(dk(aft1s)),hb,dkim
    h=gethd(hb)
    ccdt(1,0)=h.misc.temp.tcf.mean/10.
    ENDIF ELSE rdmeandk,cam,1,dkim,nointeractive=nointeractive
dkfr(256,0)=stdsize(MEDIAN(dkim,3))

; now 10-s dark frames
bef10s=MAX(WHERE((t(dk) LT firstcal) AND (STRPOS(f(dk),'010.0') NE -1)))
IF bef10s NE -1 THEN BEGIN
    rdkimg,f(dk(bef10s)),hb,dkim
    h=gethd(hb)
    ccdt(0,1)=h.misc.temp.tcf.mean/10.
    ENDIF ELSE rdmeandk,cam,10,dkim,nointeractive=nointeractive
dkfr(0,256)=stdsize(MEDIAN(dkim,3))

aft10s=MIN(WHERE((t(dk) GT lastcal) AND (STRPOS(f(dk),'010.0') NE -1)))
IF aft10s NE -1 THEN BEGIN
    rdkimg,f(dk(aft10s)),hb,dkim
    h=gethd(hb)
    ccdt(1,1)=h.misc.temp.tcf.mean/10.
    ENDIF ELSE rdmeandk,cam,10,dkim,nointeractive=nointeractive
dkfr(256,256)=stdsize(MEDIAN(dkim,3))

; and finally 60-s dark frames
bef60s=MAX(WHERE((t(dk) LT firstcal) AND (STRPOS(f(dk),'060.0') NE -1)))
IF bef60s NE -1 THEN BEGIN
    rdkimg,f(dk(bef60s)),hb,dkim
    h=gethd(hb)
    ccdt(0,2)=h.misc.temp.tcf.mean/10.
    ENDIF ELSE rdmeandk,cam,60,dkim,nointeractive=nointeractive
dkfr(0,512)=stdsize(MEDIAN(dkim,3))

aft60s=MIN(WHERE((t(dk) GT lastcal) AND (STRPOS(f(dk),'060.0') NE -1)))
IF aft60s NE -1 THEN BEGIN
    rdkimg,f(dk(aft60s)),hb,dkim
    h=gethd(hb)
    ccdt(1,2)=h.misc.temp.tcf.mean/10.
    ENDIF ELSE rdmeandk,cam,60,dkim,nointeractive=nointeractive
dkfr(256,512)=stdsize(MEDIAN(dkim,3))

; 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+'608c'
tppath(0,2)=tproot+dd+'630c'
tppath(0,3)=tproot+dd+'668c'
tppath(0,4)=tproot+dd+'520c'
tppath(1,0)=tproot+dd+'835c'
tppath(1,1)=tproot+dd+'820c'
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
    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(0,j2)-tcf) LT ABS(ccdt(1,j2)-tcf) THEN j1=0 ELSE j1=1

    ; Convert image to standard size and attempt to correct
    ; binning-induced overflow
    im=stdsize(im)
    nsat=TOTAL(im GE 16383)

    IF (NOT null_image) AND (nsat LT 100) THEN BEGIN

        ; subtract appropriate dark frame
        im=(im-dkfr(256*j1:256*(j1+1)-1,256*j2:256*(j2+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 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
	        wtsum(filt)=wtsum(filt)+(texp+0.045)*sens
	        sumwt(filt)=sumwt(filt)+(texp+0.045)
                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
        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 100) $
	    THEN MESSAGE,STRING(f(cal(j)),nsat $
		               ,FORMAT='(/,A,": ",I0.0," saturated pixels")') $
			       ,/INFORMATIONAL
	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))

    IF ok THEN BEGIN
        savename=STRING(calsroot,dd,yr MOD 100,mo,day,dc,cam,cal_suffix $
            ,FORMAT='(A,A1,3I2.2,A1,I1,A)')
        OPENW,unit,savename,/GET_LUN
        PRINTF,unit,R_s_per_DN,FORMAT='(5F6.1)'
        CLOSE,unit
        FREE_LUN,unit
        ENDIF ELSE BEGIN
	PRINT,yr MOD 100,mo,day,cam $
	     ,FORMAT='("Wonky calibration data on ",3I2.2," for camera ",I1)'
	ENDELSE
    ENDIF

RETURN
END