; DKMEAN.PRO
; calculate mean and standard deviation of a series of dark frames
; set up display
!ORDER=1
loadct,13
; create statistics arrays
si=LONARR(256,256) & sis=si
; get user input
PRINT,'Enter the day number, camera number and exposure time (s)'
READ,day,cam,texp
; find the right images
path=STRING(cam,day,texp $
     ,FORMAT='("/jasper/cnsr3_data1/dk/c",I1,"/d",I2.2,"*",I3.3,".*")')
df=FINDFILE(path)
; get the filename extensions
slen=MEDIAN(STRLEN(df))
; sort image file names in time order
ss=SORT(FIX(STRMID(df,slen-3,3)))
df=df(ss)
nn=N_ELEMENTS(df)
ni=0
st=0.0
nt=0
resp=''
WINDOW,0
WINDOW,2
; process each image
FOR i=0,nn-1 DO BEGIN
;   read image
    rdkimg,df(i),hb,img
;   get header
    kih=gethd(hb)
;   get CCD temperature
    ccdt=kih.misc.temp.tcf.mean/10.
;   remove 'hot spots'
    img=cleanmod(img,2,2,2)
;   show image ...
    WSET,0
    TVSCL,img<(160+2*texp),i MOD 4
;   ... and DN histogram
    WSET,2
    PLOT,HISTOGRAM(img,min=0,max=(160+2*texp))
    PRINT,'Accept this dark frame in mean/s.d. computation? ([Y]/N)'
    READ,resp
    IF STRUPCASE(resp) NE 'N' THEN BEGIN
;       Image ok, update statistics
        IF ccdt LE -5.0 THEN BEGIN
            st=st+ccdt
            nt=nt+1
        ENDIF
        si=si+img
        sis=sis+LONG(img)*LONG(img)
        ni=ni+1
        PRINT,df(i),ccdt,MIN(img,MAX=mxi),mxi,MIN(si,MAX=mxsi),mxsi,MIN(sis,MAX=mxsis),mxsis $
             ,FORMAT='(A,F6.1,2I6,2I7,2I12)'
    ENDIF
ENDFOR
; calculate mean and s.d. images
mt=ROUND(st/FLOAT(nt))
mi=si/FLOAT(ni)
m=FLOAT(ni*DOUBLE(sis)-DOUBLE(si)*DOUBLE(si))
sdi=SQRT(m/(ni*(ni-1)))
; save to appropriately-named file
; mean image first ...
mname=STRING(cam,cam,STRTRIM(FIX(texp),2),FIX(mt) $
      ,FORMAT='("/jasper/cnsr3_data1/dk/c",I1,"/dk",I1,".",A,"s.mean.",I3,"C")')
OPENW,unit,mname,/GET_LUN
gb=N_ELEMENTS(WHERE(mi GT 255))
IF gb GT 0.01*N_ELEMENTS(mi) THEN bmi=FIX(mi+0.5) $
                             ELSE bmi=BYTE((mi+0.5)<255)
WRITEU,unit,bmi
CLOSE,unit
FREE_LUN,unit
; ... then standard deviation image
sdname=STRING(cam,cam,STRTRIM(FIX(texp),2),FIX(mt) $
      ,FORMAT='("/jasper/cnsr3_data1/dk/c",I1,"/dk",I1,".",A,"s.10stdev.",I3,"C")')
OPENW,unit,sdname,/GET_LUN
bsdi=BYTE((10.*sdi+0.5)<255)
WRITEU,unit,bsdi
CLOSE,unit
FREE_LUN,unit
END