;+
; NAME		favgpb.pro
;
; PURPOSE	Routine to average a list of images
;
; AUTHOR	David Elmore   13 Nov 2001
;
; HISTORY	Andrew L. Stanger
;		 3 Jan 2002: Cosmetic changes only.
;		 4 Jan 2002: Plots go to postscript file.
;			     Update AVGTIME* keywords in FITS header.
;		 7 Jan 2002: Update PROCSTEP.
;-

PRO favgpb, fits_list

avgtime1 = ''
avgtime2 = ''
avgtime3 = ''
avgtime4 = ''
avgtime5 = ''
avgtime6 = ''
avgtime7 = ''
avgtime8 = ''

cavgtime1 = ''
cavgtime2 = ''
cavgtime3 = ''
cavgtime4 = ''
cavgtime5 = ''
cavgtime6 = ''
cavgtime7 = ''
cavgtime8 = ''

;--- Set up display values.

!y.title    = 'Daily Average pB'
!x.title    = 'Height (R-Sun)'
!p.title    = 'Average pB vs. Height'
!p.font     = 1
!p.thick    = 1.6
!p.ticklen  = 1
;!p.charsize = 1.6
!p.charsize = 1.2
!x.style    = 1
!y.style    = 1
!x.ticks    = 10
!x.range    = [1.10,     2.8]
!y.range    = [1.e-10, 1.e-6]

is_cal_avg  = 18
azplace     = 220

filename = 'mk4.fts'


;--- Read a header to set some variables.

CLOSE, 2
OPENR, 2, fits_list 
READF, 2, filename
CLOSE, 2

img  = readfits (filename, hdu0, /noscale)
hdr0 = readfits (filename, ext0, /noscale, exten_no=1)

numaz    = fxpar (hdu0, 'NAXIS1')	; Xdim
numpix   = fxpar (hdu0, 'NAXIS2')	; Ydim
date     = fxpar (hdu0, 'DATE-OBS')
pixoff   = fxpar (ext0, 'R0')		; r0
pixscale = fxpar (ext0, 'DR')		; dr

avgtime1 = STRTRIM (fxpar (ext0, 'AVGTIME1', cavgtime1), 2)
avgtime2 = STRTRIM (fxpar (ext0, 'AVGTIME2', cavgtime2), 2)
avgtime3 = STRTRIM (fxpar (ext0, 'AVGTIME3', cavgtime3), 2)
avgtime4 = STRTRIM (fxpar (ext0, 'AVGTIME4', cavgtime4), 2)
avgtime5 = STRTRIM (fxpar (ext0, 'AVGTIME5', cavgtime5), 2)
avgtime6 = STRTRIM (fxpar (ext0, 'AVGTIME6', cavgtime6), 2)
avgtime7 = STRTRIM (fxpar (ext0, 'AVGTIME7', cavgtime7), 2)
avgtime8 = STRTRIM (fxpar (ext0, 'AVGTIME8', cavgtime8), 2)

syear  = STRMID (date, 0, 4)
smonth = STRMID (date, 5, 2)
sday   = STRMID (date, 8, 2)
plotfile = syear + smonth + sday + '.mk4.cpb.plot.ps'

PRINT, numaz, numpix, pixscale, pixoff


;--- Heights from header of first file.

height = FINDGEN (numpix) * pixscale + pixoff

inimg  = INTARR (NUMAZ, NUMPIX)
avgimg = FLTARR (NUMAZ, NUMPIX)

CLOSE, 2
OPENR, 2, fits_list

;--- Determine number of images to average.

numimg = 0
WHILE ( NOT EOF (2) ) DO BEGIN
    numimg = numimg + 1
    READF, 2, filename
ENDWHILE

print, 'numimg: ', numimg

rays = FLTARR (numimg, NUMPIX)

;--- Open list file.

CLOSE, 2
OPENR, 2, fits_list

WINDOW, 0, xs=numaz, ys=numpix
WINDOW, 1, xs=numaz, ys=numpix
imgnum = 0

;--- File loop.

WHILE ( NOT EOF (2) ) DO BEGIN
    READF, 2, filename
    PRINT, filename

    ;--- Read FITS image.

    inimg = readfits (filename, hdu, /noscale)
    hdr   = readfits (filename, ext, /noscale, exten_no=1)

    imgtime = fxpar (hdu, 'TIME-OBS')
    itime   = STRMID (imgtime, 0, 5)
    scale   = fxpar (ext, 'DATASCL')		; data_scale
    offset  = fxpar (ext, 'DATAOFF')		; data_offset

    ;--- Apply photometric scaling.
    ;--- Accumulate image sum.

    img = (FLOAT (inimg) - offset) / scale 
    avgimg = avgimg + img
    rays [imgnum, *] = img [azplace, *] 

    imgnum = imgnum + 1

    ;--- Display image & average image.

    WSET, 0
    TV, BYTSCL (ALOG10 (img + 2.e-9), -9, -5)
    WSET, 1
    TV, BYTSCL (ALOG10 ((avgimg / imgnum) + 2.e-9), -9, -5)

    ;--- Add image time to extension header.

    inum1 = (imgnum - 1) / 3
    inum2 = (imgnum - 1) - inum1 * 3
    print, 'imgnum/inum1/inum2: ', imgnum, inum1, inum2, ' itime: ', itime

    CASE inum1 OF
    0: BEGIN						$
       avgtime1 = STRTRIM (avgtime1 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME1', avgtime1, cavgtime1
       END
    1: BEGIN						$
       avgtime2 = STRTRIM (avgtime2 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME2', avgtime2, cavgtime2
       END
    2: BEGIN						$
       avgtime3 = STRTRIM (avgtime3 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME3', avgtime3, cavgtime3
       END
    3: BEGIN						$
       avgtime4 = STRTRIM (avgtime4 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME4', avgtime4, cavgtime4
       END
    4: BEGIN						$
       avgtime5 = STRTRIM (avgtime5 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME5', avgtime5, cavgtime5
       END
    5: BEGIN						$
       avgtime6 = STRTRIM (avgtime6 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME6', avgtime6, cavgtime6
       END
    6: BEGIN						$
       avgtime7 = STRTRIM (avgtime7 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME7', avgtime7, cavgtime7
       END
    7: BEGIN						$
       avgtime8 = STRTRIM (avgtime8 + ' ' + itime, 2)
       fxaddpar, ext0, 'AVGTIME8', avgtime8, cavgtime8
       END

    ELSE: PRINT, 'imgnum: ', imgnum, ' greater than maximum allowed'
    ENDCASE

ENDWHILE

;--- Compute average image.

avgimg = avgimg / numimg 
TV, BYTSCL (ALOG10 (avgimg + 2.e-9), -9, -5)

;--- Write average image to disk in FITS format.

fxaddpar, ext0, 'PROCSTEP', is_cal_avg			; ProcSteps
outfile = STRMID (filename, 0, 6) + '.mk4.cpb.avg.fts'
writefits, outfile, FIX (avgimg * scale + offset), hdu0
writefits, outfile, hdr0, ext0, /append


WINDOW, 2, xs=512, ys=512

; What is mean and std dev for each height?
; An indicator of coronal structure.

avgheight      = FLTARR (numaz)
avgheightplus  = FLTARR (numaz)
avgheightminus = FLTARR (numaz)
avgheightmax   = FLTARR (numaz)
avgheightmin   = FLTARR (numaz)

FOR ii = 0, numpix-1 DO BEGIN
    res = MOMENT (avgimg [*, ii], sdev=derf)
    avgheight [ii]      = res [0]
    avgheightplus [ii]  = res [0] + derf
    avgheightminus [ii] = res [0] - derf
    avgheightmax [ii]   = max (avgimg [*, ii])
    avgheightmin [ii]   = min (avgimg [*, ii])
ENDFOR

writefits, 'corona_model', avgheight


;--- Create azimuthally binned data.
;    Again, this will work on average image of the day.

binsize = 20
binimg = FLTARR (numaz / binsize, numpix)

FOR jj = 0, numaz-1 DO BEGIN
	kk = jj / binsize
	binimg [kk, *] = binimg [kk, *] + img [jj, *]
ENDFOR

FOR kk = 0, numaz/binsize - 1 DO binimg [kk, *] = binimg [kk, *] / binsize
	

;--- Write plots to a Postscript file.

PRINT, "Set plot device"
SET_PLOT, 'PS'
DEVICE, filename = plotfile

;STOP

!p.title    = 'Average pB vs. Height' + '   Mk4 ' + date
plot_oo, height, avgheight
OPLOT, height, avgheightplus, psym=3
OPLOT, height, avgheightminus, psym=3

;OPLOT, height, avgheightmax, psym=3
;OPLOT, height, avgheightmin, psym=3

;STOP

;--- What is the stdev for each height over the day?
;
; This requires all the images going into the average
; image of the day to get RMS at each height
; over the day.

daystdev  = FLTARR (NUMPIX)
daymean   = FLTARR (NUMPIX)
polystdev = FLTARR (NUMPIX)
imgnum    = FINDGEN (numimg)

FOR jj = 0, numpix-1 DO BEGIN
	res  = MOMENT (rays [*, jj], sdev=derf)
	coef = POLY_FIT (imgnum, rays [*, jj], 2, yfit, yband, sigma)
	daystdev  [jj] = derf
	daymean   [jj] = res [0]
	polystdev [jj] = sigma
ENDFOR

;!p.title = 'Mean and Standard Deviation at One Azimuth'
!p.title = 'Mean & Std Deviation at One Azimuth' + '   Mk4 ' + date
!y.title = 'pB'

plot_oo, height, daymean
OPLOT, height, daystdev, psym=3
OPLOT, height, polystdev, psym=1
;STOP

;!p.title = 'pB:  10 degree binned azimuths'
!p.title = 'pB: 10 degree binned azimuths' + '   Mk4 ' + date

plot_oo, height, binimg [0, *], psym=3
FOR kk = 1, binsize - 1 DO BEGIN
	OPLOT, height, binimg [kk, *], psym=3
ENDFOR

DEVICE, /close
SET_PLOT, 'X'

END