pro scc_monthly_min,tel0,sc,dt,pol,NDAYS=ndays,FILES=files,OUTSIZE=outsize, OUTDIR=outdir, $
ROLL=roll, HISHIFT=HIshift, SMOOTHIT=smoothit, PEAKSMOOTH=peaksmooth, _EXTRA=_extra
;
;+
; $Id: scc_monthly_min.pro,v 1.32 2023/08/28 19:16:15 secchia Exp $
;
; Project : STEREO SECCHI
;
; Name : SCC_MONTHLY_MIN
;
; Purpose : This procedure generates the minimum of the daily median
; images for an interval (default is 27 days)
;
; Explanation:
;
; Use : IDL> scc_monthly_min,tel,sc,dt,pol
;
; Inputs :tel = cor1, cor2, hi1, hi2
; sc = 'A', 'B'
; dt = date to be processed
; pol= '0', '120', '240', 'tbr' (tbr = total brightness), 'dbl' (type DOUBLE), 'df' (type Deep Field)
;
; Outputs : FITS file in $SECCHI_BKG/[ab]/monthly_min/YYYY/
; Date in name is the MIDPOINT of time covered (DATE-AVG)
;
; Keywords : FILES: specify your own file list (this may not work properly yet!!)
; NDAYS= number of days to use for min (1/2 before and 1/2 after)
; OUTDIR= Location to put output files (defaults to $SECCHI_BKG/)
; OUTSIZE= Set output size LE daily medians; By default, output result is size of midpoint input
; ROLL= Make min from daily images of SC_ROLL this value
; /HISHIFT Apply hi_align_image (HI only)
; /SMOOTHIT Apply smooth(,17) to each daily image
; /PEAKSMOOTH Apply filter from sharpened image to daily image
;
;
; Calls from LASCO :
;
; Common :
;
; Restrictions: Need $SECCHI_BKG and appropriate permissions if not using savedir option
;
; Side effects:
;
; Category : DAILY
;
; Prev. Hist. : None.
;
; Written : Karl Battams, NRL/I2, Jan 2007
;
; $Log: scc_monthly_min.pro,v $
; Revision 1.32 2023/08/28 19:16:15 secchia
; committing merged revision
;
; Revision 1.31 2023/08/28 15:06:32 thompson
; Incorporate flip on 2023-08-12
;
; Revision 1.30 2023/06/05 14:16:17 secchia
; corrected DF file names when sending in file name
;
; Revision 1.29 2023/06/02 16:03:04 nathan
; add DF option to scc_monthly_min based on scc_monthly_min_ep.pro from Evangelos Paouris
;
; Revision 1.28 2017/08/10 21:53:04 nathan
; fix incorrect filters for cor2
;
; Revision 1.27 2014/09/03 17:27:45 secchia
; added a repoint day for SCA for start of sidelobe 1 and compares imsz for hsize of first image
;
; Revision 1.26 2012/03/07 21:01:29 nathan
; Add /SHOOTHIT, /PEAKSMOOTH keywords; use CRVAL median instead of mean in
; output header; put this info in COMMENT in header
;
; Revision 1.25 2011/12/27 19:45:17 nathan
; Change /NOSHIFT to /HISHIFT so hi_align_image is not default
;
; Revision 1.24 2011/12/27 19:43:27 nathan
; Allow gap of 5 days; implement hi_align_image default
;
; Revision 1.23 2011/08/15 21:54:17 nathan
; add cd,old for one case
;
; Revision 1.22 2011/08/08 17:41:03 nathan
; fix filelist conflict with FILES keyword
;
; Revision 1.21 2011/08/04 22:35:02 nathan
; Implement FILES= keyword so no other input is required, and generate filenames
; with rDEG indicating S/C roll (if GT 8 deg).
;
; Revision 1.20 2011/07/20 15:34:07 nathan
; documentation
;
; Revision 1.19 2011/05/18 15:29:25 secchib
; nr - increase max missing day interval from 2 to 3
;
; Revision 1.18 2011/05/02 20:35:03 secchib
; typo
;
; Revision 1.17 2011/04/27 17:43:10 nathan
; add cor2-b event
;
; Revision 1.16 2011/03/10 23:18:27 nathan
; Change way days are screened for inclusion to match how Bill is doing it for
; cor1, so that pointing/stray light changes are automatically screened out.
; Must edit this pro every time there is a stray light change.
;
; Revision 1.15 2011/03/08 13:06:42 secchia
; corrected premssion check when file does not exist
;
; Revision 1.14 2011/02/04 22:37:37 nathan
; check for write permission
;
; Revision 1.13 2011/01/19 18:54:13 nathan
; bunit should be dn/s (unchanged)
;
; Revision 1.12 2010/05/12 15:53:36 nathan
; Remove /NOREBIN and add OUTSIZE= keywords
;
; Revision 1.11 2009/04/09 12:03:48 mcnutt
; added cor2 pol=dbl backgrounds
;
; Revision 1.10 2008/12/03 18:08:01 nathan
; some corrections to days logic; mkdir so secchia can write to it
;
; Revision 1.9 2008/05/21 15:52:14 nathan
; added OUTDIR keyword and increased robustness
;
; Revision 1.8 2008/03/13 11:52:20 mcnutt
; change input tel to tel0 so not to change case of telj in secchi_reduce.pro
;
; Revision 1.7 2007/11/19 19:22:05 nathan
; use normal SECCHI_BKG directories
;
; Revision 1.6 2007/11/14 17:39:36 nathan
; Change mindays to 9 instead of 15
;
; Revision 1.5 2007/11/13 21:52:48 nathan
; added more info if not enough days
;
; Revision 1.4 2007/10/26 22:29:20 nathan
; DATE-OBS in header is start of interval; midpoint is DATE-AVG; saves output
; in $SECCHI_BKG/../newbkg; use secchi header structure for header and compute
; same averages as in scc_mk_daily_med.pro
;
; Revision 1.3 2007/10/11 16:37:12 reduce
; Trying to fix bugs where processing is abandoned in odd dirs. Karl.
;
; Revision 1.2 2007/07/05 19:07:28 reduce
; Fixed 0-deg polarizer filename error. Karl B.
;
; Revision 1.1 2007/02/27 21:03:03 reduce
; Initial Release -- Karl Battams, NRL/Interferometrics
;
;
;-
;+
version='$Id: scc_monthly_min.pro,v 1.32 2023/08/28 19:16:15 secchia Exp $'
len=strlen(version)
version=strmid(version,1,len-2)
dlm=get_delim()
good_files=0 ; just a counter
IF keyword_set(FILES) THEN BEGIN
filelist=files
jk=sccreadfits(filelist,allh,/nodata)
h0=allh[0]
tel0=h0.detector
sc=rstrmid(h0.obsrvtry,0,1)
srt=sort(allh.date_obs)
allh=allh[srt]
filelist=filelist[srt]
ndy=n_elements(filelist)
mindays=ndy
hdr0=allh[ndy/2]
; find middle date_obs, not necessarily the middle file
tai0=anytim2tai(h0.date_obs)
dt=utc2yymmdd(tai2utc(tai0+((anytim2tai(allh[ndy-1].date_end)-tai0)/2) ),/yyyy)
IF pol ne 'df' then begin ;--do not change pol or set roll value if pol='df'---added by LCH ---Paouris
IF tel0 EQ 'EUVI' THEN pol=trim(h0.wavelnth) ELSE $
IF tel0 EQ 'HI1' or tel0 EQ 'HI2' THEN pol='tbr' ELSE $
IF h0.seb_prog EQ 'DOUBLE' THEN pol='dbl' ELSE $
IF h0.polar GT 999 THEN pol='tbr' ELSE pol=trim(h0.polar)
IF h0.sc_roll LT 0 THEN addt=360 else addt=0
roll=round(avg(allh.sc_roll+addt))
ENDIF
polar=1
imsz=h0.naxis1
;---
ENDIF ELSE BEGIN
mindays=9 ; this is superceded by NDAYS
IF KEYWORD_SET(NDAYS) THEN BEGIN
IF ndays mod 2 EQ 0 THEN BEGIN
ndays=ndays+1
message,'NDAYS must be odd; using '+trim(ndays)+' days.',/info
ENDIF
ndy=ndays
mindays=ndays
ENDIF ELSE ndy=27
ENDELSE
; Do some checks...
sc=strlowcase(sc)
scuc=strupcase(sc)
tel=strlowcase(tel0)
IF strmid(tel,0,2) EQ 'co' THEN iscor=1 ELSE iscor=0
IF strmid(tel,0,2) EQ 'hi' THEN ishi =1 ELSE ishi =0
matchpos=0
IF keyword_set(HISHIFT) and ishi THEN matchpos=1
if ~ishi and ~iscor THEN message,'Input 1 must be "cor1","cor2","hi1",or "hi2"
IF (STRLEN(dt) LT 9) THEN BEGIN
; input is yymmdd format
udt = yymmdd2utc(dt)
ENDIF ELSE BEGIN
udt=anytim2utc(dt)
ENDELSE
dt = utc2yymmdd(udt)
dir_date=strmid(utc2yymmdd(udt,/yyyy),0,6)
date=utc2str(udt,/date_only)
; end checks
;
; Set directory for input and output; cd to daily_med directory.
;
dir = GETENV('SECCHI_BKG')+'/'+sc+dlm+'daily_med'+dlm ;+dir_date+dlm
CD,dir,curr=old
IF keyword_set(OUTDIR) THEN BEGIN
outdir1=outdir
IF strmid(outdir,0,1) EQ '.' THEN BEGIN
len=strlen(outdir)
outdir1=concat_dir(old,strmid(outdir,2,len-2))
ENDIF
ENDIF ELSE $
outdir1= GETENV('SECCHI_BKG')+dlm+sc+dlm+'monthly_min'+dlm+dir_date
if not file_exist(outdir1) THEN BEGIN
cmd='mkdir -p -m 775 '+outdir1
spawn,cmd,/sh
endif
; Get some constants...
IF ishi THEN pol='tbr'
IF ishi THEN units='deg' ELSE units='arcsec'
if strlowcase(pol) EQ 'tbr' THEN polstr='pTBr' ELSE BEGIN
pp=float(pol)
if (pp EQ 0) then polstr='p000' ELSE polstr='p'+trim(pol)
ENDELSE
IF strlowcase(pol) EQ 'dbl' then polstr='dbTB'
IF strlowcase(pol) EQ 'df' then polstr='dbDF' ;--- This is the new 'dbDF' for filename -Paouris 05/30/2023
tel_str=strmid(tel,0,1)+strmid(tel,strlen(tel)-1,1) ; EG hi1 becomes h1, cor2 becomes c2
srchstring=tel_str+scuc+'_'+polstr+'_' ; e.g. hi1_pTBr_
srchstring_part1='d'+srchstring
IF keyword_set(ROLL) THEN BEGIN
IF roll LT 0 THEN roll=roll+360
postd = 'r'+STRING(round(roll),format='(I3.3)')
ENDIF ELSE postd=''
; end constants
; Setup some initial values
mjd = STR2UTC (date +' 12:00:00')
;mjd.mjd = mjd.mjd-ndy
mjd0 = mjd
mjd0.mjd = mjd.mjd - ndy/2
first = 1
totim = 0
totdays=0
tutc = mjd0
endutc = mjd0
ilast=0
diff1=0
pwd
mname = 'm'+srchstring+dt+postd+'.fts'
outname=concat_dir(outdir1,mname)
file_acc, outname, exi, rea, wri, exe, typ
IF exi and not (wri) THEN BEGIN
print,'Permission error: ',outname
message,'Not writing '+mname,/info
stop
ENDIF
; end intial value setup
; do the rest...
if datatype(filelist) EQ 'UND' THEN filelist=strarr((ndy+1)/2) ELSE goto, readfiles
; Set limits based on the major spacecraft repointings, or other major events.
;
; These dates should be identical to values in scc_getbkgimg.pro.
; COR1 is not included here--see cor1_monthly_min.pro.
;
case scuc of
'A': begin
repoint = ['2006-12-21T13:15', '2007-02-03T13:15','2014-08-19T00:00', $
'2023-08-12']
endcase
'B': begin
repoint = ['2007-02-03T18:20', '2007-02-21T20:00']
if strupcase(tel) eq 'COR2' then repoint = $
[repoint, '2010-02-23T08:12', '2011-01-27T03:47','2011-04-25t18:30']
endcase
endcase
repoint = anytim2utc(repoint)
tai_repoint = utc2tai(repoint)
taiin = utc2tai(date)
i1 = max(where(tai_repoint lt taiin, tcount))
if tcount gt 0 then mjdmin = repoint[i1].mjd else mjdmin = 0
i2 = min(where(tai_repoint gt taiin, tcount))
if tcount gt 0 then mjdmax = repoint[i2].mjd else mjdmax = 99999
;
; ** This first loop tests available days to determine the correct ndy and stdate
; There must be same number of days before and after mjd.
;
maxgap=5 ; days
hdrset=0
FOR i=0,ndy-1 DO BEGIN
; check for pointing/stray light change
IF (tutc.mjd gt mjdmin) and (tutc.mjd lt mjdmax) then begin
day = UTC2yymmdd (tutc,/yyyy)
monthdir = strmid(day,0,4) + strmid(day,4,2)
name = concat_dir(monthdir,srchstring_part1+strmid(day,2,6)+postd+'.fts')
print,name,file_exist(name),i
IF file_exist(name) THEN BEGIN
; make sure same number of days before as after
idiff=i-ilast
;print,name,idiff
;printf,lunm,'Testing ',name,n
;IF n GT 0 or (idiff GT 6) THEN BEGIN
IF NOT keyword_set(NDAYS) THEN BEGIN
; max 3 day gap
; ** Second half of interval
IF idiff GT maxgap+1 AND tutc.mjd GT mjd.mjd THEN GOTO, loopexit
; ** First half of interval
IF idiff GT maxgap+1 AND tutc.mjd LT mjd.mjd THEN first = 1
ENDIF
IF first THEN BEGIN
diff1 = mjd.mjd - tutc.mjd
;IF NOT keyword_set(NDAYS) AND diff1 LT 7 OR $
; (keyword_set(NDAYS) AND (tutc.mjd - mjd0.mjd) GT 1) THEN BEGIN
IF diff1 LT ((mindays/2)<((ndy/2)-1)) THEN BEGIN
message, 'Not enough first half days to make model for '+mname,/info
help,day,diff1,name,mindays
wait,4
cd,old
return
ENDIF
ff0 = name
first = 0
idiff=0
ENDIF ELSE BEGIN
ff0=[ff0,name]
;endutc = tutc
ENDELSE
ilast = i
endutc=tutc
;ENDIF
IF tutc.mjd GE mjd.mjd and NOT hdrset THEN BEGIN
; Use header of file close to midpoint.
print,'Using ',name,' for default header.'
x0=sccreadfits(name,hdr0)
imsz=hdr0.naxis1
IF keyword_set(OUTSIZE) THEN imsz=outsize
IF (hdr0.naxis1 NE imsz) THEN x=scc_putin_array(x0,hdr0,imsz,/new)
hdrset=1
ENDIF
ENDIF ELSE wait,2
ENDIF
tutc.mjd = tutc.mjd + 1
ENDFOR
loopexit:
; ** Second half of interval
diff2 = endutc.mjd - mjd.mjd
;IF NOT keyword_set(NDAYS) AND diff2 LT 7 OR $
; (keyword_set(NDAYS) AND (mjd0.mjd + ndy - endutc.mjd) GT 1) THEN BEGIN
IF diff2 LT ((mindays/2)<((ndy/2)-1)) THEN BEGIN
message, 'Not enough second half days to make model for '+tel+sc+'_p'+pol+'_'+dt+postd,/info
help,day,diff2,name,mindays
;printf,lunm,tt,' ',imgtype,' ',date,ndy,' Diff2: ', diff2
;printf,lunm,'Not enough days to make model for '+td
;close,lunm
;free_lun,lunm
wait,4
cd,old
return
ENDIF
mindiff = diff1 < diff2
IF NOT keyword_set(NDAYS) THEN BEGIN
mjd0.mjd = mjd.mjd - mindiff
ndy = mindiff*2+1
ENDIF
help,diff1,diff2,mindiff,ndy
; printf,lunm,tt,' ',imgtype,' ',date,ndy,diff1,diff2,mindiff
wait,3
; +++++++++++++++++++++++ END OF MAIN FILE CHECKING SECTION +++++++++++++++++++++++
readfiles:
tmjd = mjd0
; first day of series
first=1
notallcorrected=0
maxshift=0.
;
; ** This part computes the result based on ndy and stdate now defined
;
daysused=strarr(ndy)
FOR i=0,ndy-1 DO BEGIN
day = UTC2yymmdd (tmjd,/yyyy)
monthdir = strmid(day,0,4) + strmid(day,4,2)
IF keyword_set(FILES) THEN name=filelist[i] ELSE name = concat_dir(monthdir,srchstring_part1+strmid(day,2,6)+postd+'.fts')
ex=file_exist(name)
;stop
IF i EQ 0 THEN IF ex EQ 0 THEN timeobs = '23:59:59' ELSE timeobs = '12:00:00'
; --> first day in sequence missing
IF (ex) THEN BEGIN
;FOR f=0, n-1 DO BEGIN ;** loop over years, look at images for same days but different years
;IF KEYWORD_SET(ALL_YEARS) THEN name = ff(f)
PRINT,'Processing file '+name
a = sccREADFITS(name,h)
;-- Compile header info
if (first) THEN hall=h ELSE hall=[hall,h]
IF (first) THEN crotas=h.crota ELSE crotas=[crotas,h.crota]
IF (first) THEN biases=h.biasmean ELSE biases=[biases,h.biasmean]
IF (first) THEN exptimes=h.exptime ELSE exptimes=[exptimes,h.exptime]
; should always be 1.0 sec
IF (first) THEN offsetcrs=h.offsetcr ELSE offsetcrs=[offsetcrs,h.offsetcr]
IF (first) THEN cadences=h.cadence ELSE cadences=[cadences,h.cadence]
IF (first) THEN readtimes=h.readtime ELSE readtimes=[readtimes,h.readtime]
IF (first) THEN cleartimes=h.cleartim ELSE cleartimes=[cleartimes,h.cleartim]
IF (first) THEN pitches=h.crval1 ELSE pitches=[pitches,h.crval1]
IF (first) THEN yaws=h.crval2 ELSE yaws=[yaws,h.crval2]
IF (first) THEN nmissings=h.nmissing ELSE nmissings=nmissings+h.nmissing
IF (first) THEN nsats=h.datasat ELSE nsats=[nsats,h.datasat]
IF (first) THEN cosmicss=h.cosmics ELSE cosmicss=[cosmicss,h.cosmics]
daysused[i]='Used '+name
;--Cannot Assume input images are already 1024x1024, SUMMED=2
IF first and h.naxis1 LT imsz THEN imsz=h.naxis1
IF h.naxis1 LT imsz THEN BEGIN
inp=''
print,'Input daily median size LT chosen OUTSIZE.'
read, 'y to continue or n to exit and re-choose OUTSIZE',inp
IF inp EQ 'n' THEN BEGIN
cd,old
return
ENDIF
ENDIF
IF keyword_set(SMOOTHIT) THEN BEGIN
IF smoothit NE 1 THEN fac=smoothit ELSE fac=17
message,'smooth(,'+trim(fac)+')',/info
a=smooth(a,fac)
ENDIF
IF keyword_set(PEAKSMOOTH) THEN BEGIN
IF peaksmooth NE 1 THEN fac=peaksmooth ELSE fac=17
message,'smoothpeaks (,'+trim(fac)+', 0.3)',/info
a=smoothpeaks(a,fac,0.3,/dark,/bright)
ENDIF
a=rebin(a,imsz,imsz)
; see scc_putin_array() earlier for header update
CASE imsz OF
2048: h.summed=1
1024: h.summed=2
512 : h.summed=3
256 : h.summed=4
ELSE: message,'Invalid OUTSIZE'
ENDCASE
a0 = a/(h.EXPTIME) ; exptime should be 1.0; IPSUM should already be corrected
maxmin,a0
IF h.exptime NE 1 THEN BEGIN
print,'Exptime=',h.exptime
wait,10
ENDIF
maxshift=maxshift>h.jitrsdev ; shift in daily could be greater than output
IF (matchpos) THEN BEGIN
IF h.ravg EQ -881 THEN notallcorrected=1
a=hi_align_image(h,a0,hdr0,/zerom,/verbose)
maxshift=maxshift>h.jitrsdev
ENDIF ELSE a=a0
IF (first) THEN BEGIN
alldays=fltarr(ndy,imsz,imsz)
alldays[0,*,*]=a
img = a
first = 0
stdate =h.date_obs
ENDIF ELSE BEGIN
alldays[i,*,*]=a
wi = where (img le 0,ni)
wa = where (a le 0,na)
img1=img<a
IF (ni gt 0) THEN img1(wi) = a(wi)>img(wi)
IF (na gt 0) THEN img1(wa) = a(wa)>img(wa)
endate=h.date_end
img = img1
ENDELSE
totim = totim+1
;ENDFOR
totdays=totdays+1
mjdn=tmjd
;printf,lunm,daysused
ENDIF
tmjd.mjd = tmjd.mjd+1
ENDFOR
exptime_stdv =stdev(exptimes,exptime_mean)
bias_stdv =stdev(biases,bias_mean)
offsetcr_stdv =stdev(offsetcrs,offsetcr_mean)
crota_stdv =stdev(crotas,crota_mean)
nsat =min(nsats)
help,exptime_mean
dsatval =float(hdr0.dsatval/exptime_mean) ; normalize this also to 1 second
dateobstai =anytim2tai(stdate)
datendtai =anytim2tai(h.date_end)
date_avg =utc2str(tai2utc(dateobstai + (datendtai-dateobstai)/2))
cadence_stdv =stdev(cadences,cadence_mean)
readtime_stdv =stdev(readtimes,readtime_mean)
cleartime_stdv =stdev(cleartimes, cleartime_mean)
crval1_stdv =stdev(pitches)
crval1_mean=median(pitches)
crval2_stdv =stdev(yaws)
crval2_mean=median(yaws)
cosmics_stdv =stdev(cosmicss, cosmics_mean)
PRINT,'Minimum CROTA value: ',strcompress(string(min(crotas)))
PRINT,'Maximum CROTA value: ',strcompress(string(max(crotas)))
PRINT,'Average CROTA value: ',trim(crota_mean)
;stop
tmjd.mjd = tmjd.mjd-1
;endate=utc2str(mjdn,/date_only)
IF not (ex) THEN timeobs='00:00:00' ; --> last day in sequence missing
print
IF datatype(endutc) NE 'UND' THEN print,utc2str(endutc,/date_only)
help,endate,stdate,ex
print,date,' ',timeobs
sz = SIZE (img)
IF (sz(0) NE 0) THEN BEGIN
;stop
;mname = 'm'+strcompress(STRMID(ff2[0],1,9),/remove_all)+dt+'.fts'
expt = median(exptimes) ; again, should be 1.0
help,expt
img2 = float(expt*img>0)
;stop
; For HI-2 we have to smooth the median array, otherwise is screws up the star brightnesses
if (hdr0.DETECTOR EQ 'HI2') then BEGIN
PRINT,'HI2 image detected.'
PRINT,'Will apply median filter to the minimum before writing to disk...'
img2=median(img2,12)
ENDIF
outhdr0=scc_update_hdr(img2>0,hdr0,verbose=verbose,satmax=dsatval)
outhdr1 = rem_tag2(outhdr0,'TIME_OBS')
outhdr = struct2fitshead(outhdr1,/allow_crota,/dateunderscore)
;IF nsat LE 0 THEN it=where(median_array_out GE dsatval, nsat)
;fxaddpar,outhdr,'DSATVAL',dsatval
;fxaddpar,outhdr,'DATASAT',nsat
help,dsatval,nsat,bscale
get_utc,dte,/ecs
fxaddpar,outhdr,'DATE-OBS',stdate,' Start of first exposure in N_IMAGES (N_DAYS)'
fxaddpar,outhdr,'DATE-END',h.date_end,' End of last exposure in N_IMAGES (N_DAYS)'
fxaddpar,outhdr,'DATE-AVG',date_avg,' Midpoint between OBS and END'
fxaddpar,outhdr,'CADENCE', cadence_mean,' sec; stdev='+trim(cadence_stdv)
fxaddpar,outhdr,'READTIME',readtime_mean,' sec; stdev='+trim(readtime_stdv)
fxaddpar,outhdr,'CLEARTIM',cleartime_mean,' sec; stdev='+trim(cleartime_stdv)
IF matchpos THEN BEGIN
fxaddpar,outhdr,'JITRSDEV',maxshift,'Max hi_align_image pixel shift in x or y'
fxaddpar,outhdr,'HISTORY','Used hi_align_image.pro'
IF notallcorrected THEN fxaddpar,outhdr,'COMMENT','Some daily medians were not aligned.'
ENDIF ELSE BEGIN
fxaddpar,outhdr,'CRVAL1',crval1_mean,' '+units+' median; stdev='+string(crval1_stdv,'(f4.2)')+'; was '+trim(hdr0.crval1)
fxaddpar,outhdr,'CRVAL2',crval2_mean,' '+units+' median; stdev='+string(crval2_stdv,'(f4.2)')+'; was '+trim(hdr0.crval2)
fxaddpar,outhdr,'CROTA',crota_mean,'deg; avg; stdev='+string(crota_stdv,'(f4.2)')+'; was '+trim(hdr0.crota)
print,'CROTA diff from PC matrix:',crota_mean-hdr0.crota
ENDELSE
fxaddpar,outhdr,'COSMICS',cosmics_mean,' stdev='+trim(cosmics_stdv)
fxaddpar,outhdr,'DATE',dte
;fxaddpar,outhdr,'BUNIT','DN'
fxaddpar,outhdr,'BSCALE',1.0
fxaddpar,outhdr,'HISTORY','Data values normalized to 1 sec'
fxaddpar,outhdr,'EXPTIME', 1.0,'effective exptime; actual avg='+trim(exptime_mean)+'; stdev='+trim(exptime_stdv)
fxaddpar,outhdr,'BIASMEAN',bias_mean,'=avg over N_IMAGES'
fxaddpar,outhdr,'BIASSDEV',bias_stdv,' Over N_IMAGES'
fxaddpar,outhdr,'OFFSETCR',offsetcr_mean,'=avg over N_IMAGES; stdev='+trim(offsetcr_stdv)
fxaddpar,outhdr,'N_IMAGES',totim,' number of daily median images used'
; stop
FXADDPAR,outhdr,'FILENAME',mname
FXADDPAR,outhdr,'MID_DATE',mjd.mjd
FXADDPAR,outhdr,'MID_TIME',mjd.time
FXADDPAR,outhdr,'N_DAYS',totdays,' Daily median images used to compute this image'
FXADDPAR,outhdr,'HISTORY',version
fxaddpar,outhdr,'COMMENT','%SCC_MONTHLY_MIN: header values from middle day unless noted'
wkwd=[keyword_set(SMOOTHIT), keyword_set(PEAKSMOOTH), keyword_set(HISHIFT), keyword_set(ROLL), keyword_set(FILES)]
kwds=['/SMOOTHIT','/PEAKSMOOTH','/HISHIFT','/ROLL','/FILES']
wkw=where(wkwd)
IF wkw[0] GE 0 THEN fxaddpar,outhdr,'COMMENT','Keywords used '+arr2str(kwds[wkw])
PRINT,'Writing monthly fits file: ',outname,' using',totim,' out of',ndy,' days.'
wait,2
WRITEFITS,outname,img2,outhdr
ENDIF
help,old
CD,old
RETURN
END