;!order=1
;loadct,13
; define origin of time: 0 UT, 1 Jan 1993
jsjan0=ymds2js(1993,1,1,0)
; load list of all January calibration file names
; openr,1,'/jasper/cnsr3_data1/calstats/cal9301'
; load list of all February calibration file names
openr,1,'/jasper/cnsr3_data1/calstats/cal9302'
cf=strarr(572)
readf,1,cf
close,1
; trim whitespace off the ends of the file names
cf=strtrim(cf,2)
; define mask to select image region showing PLBS screen
prad=42
d256=SHIFT(dist(256),128,128)
mask=bytarr(256,256)
wprad=WHERE(d256 LE prad) & mask(wprad)=1
; load all mean dark images for both cameras
dk01=bytarr(256,256) & openr,1,'dk/c0/dk0.1s.mean' & readu,1,dk01 & close,1
dk010=dk01 & openr,1,'dk/c0/dk0.10s.mean' & readu,1,dk010 & close,1
dk060=dk01 & openr,1,'dk/c0/dk0.60s.mean' & readu,1,dk060 & close,1
dk11=dk01 & openr,1,'dk/c1/dk1.1s.mean' & readu,1,dk11 & close,1
dk110=dk01 & openr,1,'dk/c1/dk1.10s.mean' & readu,1,dk110 & close,1
dk160=dk01 & openr,1,'dk/c1/dk1.60s.mean' & readu,1,dk160 & close,1
; semiarbitrarily reassign values to all pixels that might be hot spots
;spike=where(dk01 gt 158) & dk01(spike)=dk01(spike)-45
;spike=where(dk010 gt 176) & dk010(spike)=dk010(spike)-43
;spike=where(dk060 gt 213) & dk060(spike)=dk060(spike)-47
;spike=where(dk11 gt 190) & dk11(spike)=dk11(spike)-71
;spike=where(dk110 gt 195) & dk110(spike)=dk110(spike)-71
;spike=where(dk160 gt 240) & dk160(spike)=dk160(spike)-61
; create plotting windows
 ;window,2
 ;window,0
; open output file
openw,outf,'/jasper/cnsr3_data1/calstats/cal9302.stats',/GET_LUN
printf,outf,'/jasper/cnsr3_data1/cal/c?/f?/File Name   MeanCt  StDev   TCCD   Text JulSec93'
for i=0,571 do begin
    ; read i'th calibration image file
    rdkimg,cf(i),hb,img
    kih=gethd(hb)
    ; extract CCD and external temperatures
    tccd=kih.misc.temp.tcf.mean/10.
    text=kih.misc.temp.tee.mean/10.
    ; extract image date, time and express as offset Julian Seconds 
    ftime=kih.misc.tm
    futsec=(ftime.hour*60.+ftime.minute)*60.+ftime.second
    fdate=kih.misc.dt
    fyear=fdate.year
    fmonth=fdate.month
    fday=fdate.day
    fjs=ymds2js(fyear,fmonth,fday,futsec)-jsjan0
    ; generate name of appropriate dark image to subtract
    fname=strmid(cf(i),strlen(cf(i))-12,12)
    texp=strtrim(fix(strmid(fname,5,3)),2)
    cam=strmid(fname,3,1)
    dkimg='dk'+cam+texp
    ; subtract appropriate dark image
    retval=execute('img=img-'+dkimg)
    ; display image in window 0, autoscaled to 0-255 full scale
     ;wset,0
     ;tvscl,byte(img/64),i mod 4
    mmin=min(img(wprad),max=mmax)
    hbin=0.5
    repeat begin
        hbin=FIX(hbin*2)
        nx=((mmax-mmin)/hbin)+1
        hmin=mmax-nx*hbin+1
        ; define abscissa array for subimage histogram
        xx=indgen(nx)*hbin+hbin/2+hmin
        ; generate histogram of subimage showing PLBS screen
        tt=histogram(img(wprad),min=hmin,max=mmax,bin=hbin)
        print,hbin,nx,hmin,mmax,max(tt)
    endrep until max(tt) gt 50 
    ; plot histogram of selected region of image in window 2
    ;WSET,2
    ;PLOT,xx,tt,XRANGE=[hmin,mmax],XSTY=1
    ; calculate "center of gravity" of histogram
    cg=TOTAL(xx*tt)/TOTAL(tt)
    ; get indices of abscissa values greater than cg
    wgcg=WHERE(xx GT cg)
    ; total # of pixels greater than cg; need this often
    ttwgcg=TOTAL(tt(wgcg))
    ; "mean DN" = "center of gravity" of histogram above cg -
    ; this should be very close to the mean DN of the screen image if
    ; the count rate is high enough
    ww=WHERE(img*mask GT cg)
    nww=N_ELEMENTS(ww)
    mdn=TOTAL(LONG(img(ww)))/nww
    ; index to abscissa value closest to mdn
    mindiff=MIN(ABS(mdn-xx),wmdn)
    ; index to abscissa value closest to cg
    mindiff=MIN(ABS(cg-xx),wcg)
    ; sum of 1% of histogram width centered on cg
    tcg=TOTAL(tt(wcg-nx/200:min([wcg+nx/200,nx-1])))/ $
             N_ELEMENTS(tt(wcg-nx/200:min([wcg+nx/200,nx-1])))
    ; sum of 1% of histogram width centered on mdn, making sure we don't
    ; run off the end of the histogram array
    tmdn=TOTAL(tt(wmdn-nx/200:min([wmdn+nx/200,nx-1])))/ $
             N_ELEMENTS(tt(wmdn-nx/200:MIN([wmdn+nx/200,nx-1])))
    ; if average # of pixels per DN with values near mdn is significantly
    ; more than # of pixels near cg, then there is a peak at high DN,
    ; which corresponds to the image of the PLBS screen; output stats
    ; for this peak
    IF tmdn GT 2.*tcg THEN BEGIN
        ; calculate statistics of high peak, starting with mean square DN
        msdn=TOTAL(LONG(img(ww))^2)/nww
        ; then calculate standard deviation in mean DN
        sdmdn=SQRT((msdn-mdn*mdn)*nww/(nww-1))
        ;OPLOT,[mdn,mdn],!Y.CRANGE
        PRINT,cf(i),mdn,' (high peak, statistics)',sdmdn,MAX(tt(wgcg))
        PRINTF,outf,cf(i),mdn,sdmdn,tccd,text,fjs,FORMAT='(A,2F7.0,2F7.1,I9)'
    ENDIF ELSE BEGIN
        ; otherwise there is only one peak in the histogram, at fairly
        ; low DN, close to cg - output stats about this peak (from
        ; gaussfit)
        gf=gaussfit(xx,tt,aa)
        cgf=aa(1)
        sdcgf=aa(2)
        ;OPLOT,xx,gf,PSYM=10
        ;OPLOT,[cgf,cgf],!Y.CRANGE
        PRINT,cf(i),cgf,' (middle peak, gaussfit)',sdcgf
        PRINTF,outf,cf(i),cgf,sdcgf,tccd,text,fjs,FORMAT='(A,2F7.0,2F7.1,I9)'
    ENDELSE
ENDFOR
CLOSE,outf
FREE_LUN,outf
END