pro scc_mk_daily_med,tel,sc,Date,polar=polar, DOREBIN=dorebin, SAVEDIR=savedir,FILES=files, VERBOSE=verbose ; ;+ ; $Id: scc_mk_daily_med.pro,v 1.19 2008/12/11 19:18:51 nathan Exp $ ; ; Project : STEREO SECCHI ; ; Name : SCC_MK_DAILY_MED ; ; Purpose : This procedure generates an image by taking all the files of a given ; type (up to 40) for one day and finding the median value for each pixel. ; ; Explanation: ; ; Use : IDL> scc_mk_daily_med,'CAM','SC','YYYYMMDD',polar='POL' ; ; Inputs :CAM = cor1, cor2, hi1, hi2 ; SC = 'A', 'B' ; YYYYMMDD = date to be processed ; polar= '0', '120', '240', 'tbr' (tbr = total brightness) ; ; Outputs : FITS file in $SECCHI_BKG/[ab]/daily_med/YYYYMM/ ; ; Keywords : SAVEDIR: specify your own directory to save file to ; FILES: specify your own file list (this may not work properly yet!!) ; /DOREBIN output 512x512 instead of 1024x1024 ; ; ; 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_mk_daily_med.pro,v $ ; Revision 1.19 2008/12/11 19:18:51 nathan ; do not include zero values in median ; ; Revision 1.18 2008/08/06 19:16:58 nathan ; correct cor1 inp directory; fix bug if savedir=. ; ; Revision 1.17 2008/02/12 14:55:18 secchia ; nr - added /new to scc_putin_array call ; ; Revision 1.16 2008/02/11 15:31:06 secchia ; nr - fixed case where ng=0 ; ; Revision 1.15 2007/11/19 19:22:05 nathan ; use normal SECCHI_BKG directories ; ; Revision 1.14 2007/11/01 21:49:14 nathan ; added image roll check and dstart/stop correction ; ; Revision 1.13 2007/10/26 22:24:12 nathan ; Added /DOREBIN; added cadence, readtime, cleartim, crval, nmissing, cosmics ; to output header; tried to get dsatval right; use start of first image for ; DATE-OBS and midpoint is DATE-AVG; scale COR2 images to type INT before ; saving; save HI as type FLOAT; saves results in $SECCHI_BKG/../newbkg ; ; Revision 1.12 2007/09/25 22:55:52 nathan ; mostly done--still need to update scc_monthly_min.pro ; ; Revision 1.11 2007/09/25 16:00:11 nathan ; This has been re-written to utilize secchi_prep and the full SECCHI FITS header ; and to resolve Bug 224 ; ; Revision 1.10 2007/08/16 15:08:12 nathan ; made sure default output dir is group writable ; ; Revision 1.9 2007/07/25 14:15:05 reduce ; needed uppercase S/C identifier. Karl B ; ; Revision 1.8 2007/07/05 19:03:36 reduce ; Fixed 0-deg polarizer filename error. Karl B. ; ; Revision 1.7 2007/07/02 19:28:26 reduce ; Added a 'cd,orig' cmd. Karl B. ; ; Revision 1.6 2007/06/27 14:01:15 reduce ; Exits gracefully if too few files found. Karl B. ; ; Revision 1.5 2007/06/22 20:46:25 nathan ; cd,orig before return in Error case ; ; Revision 1.4 2007/06/21 19:20:00 reduce ; No longer crashes on missing images. Karl B. ; ; Revision 1.3 2007/06/18 20:48:57 reduce ; Couple of bug fixes. Should be ready for release now. Karl B. ; ; Revision 1.2 2007/02/01 20:43:52 reduce ; More mods and bug fixes. KB ; ; Revision 1.1 2007/01/30 20:29:21 reduce ; Initial Release by KB ; ;- ; ; PROCEDURE: ; For each image that satifies the selection conditions, (default naxis1=1024, ; filter and polarizer as requested), the median image is computed of ; the median value of all the images for a single day after being ; normalized to the median exposure time. ; ; If the number of images is less than 7, there is a second pass. ; In the second pass, images within +/- 2 days of the given day are used, ; up to 15 per day. ; ; ; Convert the telescope number into lower case ; And select the standard size parameters ; IF (STRLOWCASE(tel) EQ 'cor1') THEN BEGIN PRINT,' PRINT,'**********************************************************' PRINT,' Creation of median images for COR1 is not yet implemented PRINT,' as per the wishes of the COR1 instrument team. PRINT,'**********************************************************' PRINT,' return ENDIF IF (n_params() LT 3) then BEGIN PRINT,'Incorrect number of input parameters.' PRINT,' PRINT,'Usage: PRINT,' scc_mk_daily_med,CAM,S/C,DATE,POLAR PRINT,'EXAMPLE: PRINT,' IDL> scc_mk_daily_med,"hi1","a","20070115",polar="tbr"' PRINT,' --- OR --- PRINT,' IDL> scc_mk_daily_med,"cor2","b","20070222",polar="120"' PRINT,' ENDIF version='$Id: scc_mk_daily_med.pro,v 1.19 2008/12/11 19:18:51 nathan Exp $' len=strlen(version) version=strmid(version,1,len-2) n_median_min = 5 ; number of images needed to generate median cam = STRLOWCASE(tel) xaxis=1 yaxis=1 cor_flag=0 ; set to 1 if CORs are chosen if not keyword_set(polar) then BEGIN PRINT,'Keyword "polar" needs to be defined. PRINT,'Options are as follows:' PRINT,' "tbr" -- total brightness (COR or HI) PRINT,' "0", "120" or "240" -- polarized (COR only) ENDIF pol=strlowcase(polar) IF keyword_set(VERBOSE) THEN strn = ' Using:' ELSE strn = '' if (strlen(date) ne 8) THEN BEGIN PRINT,'Date must be of the form YYYYMMDD' RETURN ENDIF CASE cam OF 'cor1': BEGIN if pol EQ 'tbr' THEN BEGIN s=GETENV('SECCHI_P0')+'/'+strcompress(strlowcase(sc),/remove_all)+'/cor1/' pol_str='pTBr' cor_flag=1 ENDIF ELSE BEGIN s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/seq/cor1/' pp=strcompress(pol,/remove_all) IF (float(pp) EQ 0) THEN pol_str='p000' ELSE pol_str='p'+pp ENDELSE pref='c1' ; filename prefix IF not keyword_set(anysize) THEN BEGIN xaxis=2048 yaxis=2048 ENDIF END 'cor2': BEGIN stdexptime=4. ; am I using this? if pol EQ 'tbr' THEN BEGIN s=GETENV('SECCHI_P0')+'/'+strcompress(strlowcase(sc),/remove_all)+'/cor2/' pol_str='pTBr' cor_flag=1 ENDIF ELSE BEGIN s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/seq/cor2/' pp=strcompress(pol,/remove_all) IF (float(pp) EQ 0) THEN pol_str='p000' ELSE pol_str='p'+pp ENDELSE pref='c2' ; filename prefix IF not keyword_set(anysize) THEN BEGIN xaxis=2048 yaxis=2048 ENDIF END 'hi1': BEGIN IF KEYWORD_SET(polar) THEN pol=strupcase(polar) ELSE pol='' s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/hi_1/' pref='h1' pol='tbr' pol_str='pTBr' IF not keyword_set(anysize) THEN BEGIN xaxis=1024 yaxis=1024 ENDIF END 'hi2': BEGIN pol='' s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/hi_2/' pref='h2' pol='tbr' pol_str='pTBr' IF not keyword_set(anysize) THEN BEGIN xaxis=1024 yaxis=1024 ENDIF END ELSE: BEGIN PRINT,'Unrecognized telescope code: '+cam RETURN END ENDCASE img_type='' IF (datatype(date) NE 'STR') THEN BEGIN PRINT,'Input date must be a string E.G. 20061223' RETURN ENDIF dt2=strmid(date,0,4)+'-'+strmid(date,4,2)+'-'+strmid(date,6,2) ; stop IF keyword_set(SAVEDIR) THEN BEGIN ; in case savedir= '.' cd,savedir spawn,'pwd',outp,/sh savedir=outp[0] ENDIF ; Check if data exists yet... ; s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/'+imgdir ; got here ************** CD,s,curr=orig f=file_search(date+'/*fts') sz=size(f) IF (sz(0) EQ 0) THEN BEGIN PRINT,'No directory for '+cam+' telescope on '+date CD,orig RETURN ENDIF spawn,'pwd',/SH ; Find appropriate files... IF keyword_set(FILES) THEN BEGIN fnames = files ng=n_elements(files) good = indgen(ng) print,'Using ',trim(ng),' files for median.' ENDIF ELSE BEGIN if (not cor_flag) then summary=scc_read_summary(DATE=dt2,SPACECRAFT=sc,TELESCOPE=strupcase(tel)) $ else summary=scc_read_summary(DATE=dt2,/totalb,SPACECRAFT=sc) ;stop PRINT,'Found a total of ',strcompress(n_elements(summary)),' files for the day in the summary file.' IF (n_elements(summary) LE 5) THEN BEGIN PRINT,'Based on the summary file, there do not appear to be enough files for making a median for today on this telescope.' CD,orig return ENDIF IF datatype(summary) EQ 'INT' THEN BEGIN CD,orig return ENDIF ; Either an individual COR polarizer angle is picked, or a HI image (no polarizer) or total brightness is picked. if (pol NE 'tbr') then good = where(summary.VALUE EQ pol and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark' and summary.PROG NE 'Doub', ng) $ else if not cor_flag then good = where(summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.DEST EQ 'SSR1'and summary.PROG NE 'Dark', ng) $ else good = where(summary.VALUE EQ 1001 and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark', ng) print,'Found ',trim(ng),' appropriate images for ',dt2,' from summary file.' IF ng EQ 0 THEN BEGIN print,'ERROR: There are not enough acceptable images to make a good daily median image. Returning.' wait,5 CD,orig return ENDIF ELSE files=summary[good].FILENAME ENDELSE path=s+date+'/' CD,path,curr=old print,'Screening file headers for acceptability,,,' ; - Figure out what roll should be ut=anytim2utc(dt2) uta=replicate(ut,5) uta.mjd=uta.mjd + indgen(5)-2 ra=-1*get_stereo_roll(uta,sc,/verbose) minr=min(ra,max=maxr) help, minr,maxr ; i=1 f=0 ;stop WHILE i LE 40 and f LT ng DO BEGIN ; There are more files than we need so lets take the first and last 20 of the day. ;tmp1=files[0:19] ;tmp2=files[n-20:n-1] ;files=[tmp1,tmp2] ;print,'reading ',files[f] jnk=sccreadfits(files[f],h,/nodata) IF h.sc_roll GT maxr or h.sc_roll LT minr THEN $ print,'Not using ',files[f],' because SC_ROLL=',h.sc_roll $ ELSE BEGIN IF i EQ 1 THEN gfiles=files[f] ELSE gfiles=[gfiles,files[f]] i=i+1 ENDELSE f=f+1 ENDWHILE files=gfiles n=n_elements(files) if n LT n_median_min then begin print,'ERROR: There are not enough acceptable images to make a good daily median image. Returning.' wait,5 CD,orig return endif ELSE PRINT,'Using ',trim(n),' files for median.' datacube= fltarr(1024,1024,n) crotas= fltarr(n) biases= fltarr(n) exptimes= dblarr(n) offsetcrs= fltarr(n) cadences= fltarr(n-1) readtimes= dblarr(n) cleartimes= dblarr(n) pitches= fltarr(n) yaws= fltarr(n) nmissings= lonarr(n) nsats= lonarr(n) cosmicss= lonarr(n) midpoint=n/2 n_used=n ; used later ; This is a temporary bug fix for the COR total-B summary files last_char=strmid(files[0],strlen(files[0])-1,1) if (last_char) NE 's' then files=files+'s' ; end bug fix IF keyword_set(VERBOSE) THEN quiet=0 ELSE quiet=1 FOR m=0,n-1 DO BEGIN if file_exist(files[m]) THEN fn=files[m] ELSE $ if file_exist(files[m-1]) THEN BEGIN fn=files[m-1] PRINT,'Could not find file: ',files[m] PRINT,'Will use previous file instead...' ENDIF ELSE BEGIN PRINT,'OK, I can not find the previous file either. I give up.' CD,orig return ENDELSE print,fn IF (cor_flag) THEN BEGIN im = sccreadfits(fn,hdr) ; im is type float IF hdr.naxis1 EQ 2048 THEN BEGIN hdr.dstart1=1 hdr.dstop1 =2048 ; some early TBr images are incorrect ENDIF im = SCC_PUTIN_ARRAY(im,hdr,1024,/new) ss='' k=0 WHILE ss NE 'EXPTIMEs' DO BEGIN k=k+1 ss=strmid(hdr.comment[k],0,8) ENDWHILE print,hdr.comment[k] parts=strsplit(hdr.comment[k],' ,',/extract) recexp=avg(float(parts[1:3])) effexp=1.0 ; TBr from SECCHI_P0 are effectively 1 second dsatval=30000. ; approximate for 6sec exposure, from 1 image example (20071005_025230_0B4c2B.fts) ENDIF ELSE BEGIN secchi_prep,fn, hdr,im, /NOCALFAC, /calimg_off, /update_hdr_off, /EXPTIME_OFF, OUTSIZE=1024, SILENT=quiet ; image is corrected for IP including IPSUM ; Do not correct for EXPTIME in secchi_prep because of weighting done in hi_correction.pro recexp=hdr.exptime ;help,recexp ;wait,2 im=im/recexp effexp=1.0 dsatval=hdr.dsatval IF dsatval LE 0 THEN BEGIN IF strmid(pref,0,1) EQ 'h' THEN ccdsat=14000. ELSE ccdsat=15000. dsatval = long(ccdsat * hdr.n_images * (2^(hdr.IPSUM - 1))^2) IF strmid(pref,0,1) NE 'h' THEN dsatval = dsatval<60000 ENDIF ENDELSE thistai=anytim2tai(hdr.date_obs) IF (m GT 0) THEN cadences[m-1]= thistai-lasttai lasttai=thistai readtimes[m]= hdr.readtime cleartimes[m]= hdr.cleartim pitches[m]= hdr.crval1 yaws[m]= hdr.crval2 nmissings[m]= hdr.nmissing cosmicss[m]= hdr.cosmics crotas[m]=hdr.crota biases[m]=hdr.biasmean exptimes[m]=recexp offsetcrs[m]=hdr.offsetcr ; bias that has been subtracted nsats[m]=hdr.datasat datacube[*,*,m]=im ;--Datacube is normalized to 1 second ; ; stop if (m EQ midpoint) THEN outhdr0 = hdr if (m EQ 0) THEN dateobstai=thistai if (m EQ 0) THEN date_obs=hdr.date_obs ENDFOR ; This is experimental but what I'm trying to do here is remove from the datacube any images ; that would not make a good median image. So I start by removing images that are rolled more ; than 5-sigma of the median roll value for that day. Later I will remove pixels that are ; freakishly bright compared to the others. Note that 5-sigma might not be tight enough... we'll see... ;crota_stdv=stddev(hdrs.crota) ;crota_med=median(hdrs.crota) ;good=where( abs(hdrs.crota - crota_med) LT (5*crota_stdv) ) <<<<<<< scc_mk_daily_med.pro IF datatype(summary) EQ 'INT' THEN BEGIN CD,orig return ENDIF ; Either an individual COR polarizer angle is picked, or a HI image (no polarizer) or total brightness is picked. if (pol NE 'tbr') then good = where(summary.VALUE EQ pol and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark' and summary.PROG NE 'Doub') $ else if not cor_flag then good = where(summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.DEST EQ 'SSR1'and summary.PROG NE 'Dark') $ else good = where(summary.VALUE EQ 1001 and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark') ; Before we continue, just make sure we haven't screwed up and dropped too many images... if n_elements(good) LT n_median_min then begin print,'ERROR: There are not enough images to make a good daily median image.' CD,orig return endif else files=summary[good].FILENAME ENDELSE n=n_elements(files) PRINT,'Found a total of ',strcompress(n),' appropriate files for the day.' path=s+date+'/' CD,path,curr=old if (n gt 40) THEN BEGIN ; There are more files than we need so lets take the first and last 20 of the day. tmp1=files[0:19] tmp2=files[n-20:n-1] files=[tmp1,tmp2] n=n_elements(files) PRINT,'Trimming down number of files to ',strcompress(n) ENDIF datacube=lonarr(xaxis,yaxis,n) crotas=fltarr(n) midpoint=n/2 n_used=n ; used later ; This is a temporary bug fix for the COR total-B summary files last_char=strmid(files[0],strlen(files[0])-1,1) if (last_char) NE 's' then files=files+'s' ; end bug fix FOR m=0,n-1 DO BEGIN next: if file_exist(files[m]) THEN im=readfits(files[m],hdr) $ ELSE BEGIN ;maillist='karl.battams@nrl.navy.mil sungrazer@nrl.navy.mil' ;spawn,'mailx -s SCC_MK_DAILY_MED: Missing image found '+maillist+' < '+files[m],/sh PRINT,'Could not find file: ',files[m] PRINT,'Will use previous file instead...' if file_exist(files[m-1]) THEN im=readfits(files[m-1],hdr) $ ELSE BEGIN PRINT,'OK, I can not find the previous file either. I give up.' CD,orig return ENDELSE ENDELSE crotas[m]=float(sxpar(hdr,'CROTA')) im = scc_img_trim(im, hdr) datacube[*,*,m]=im ; stop if (m EQ midpoint) THEN BEGIN ;stop if datatype(hdr) NE 'STC' THEN h=fitshdr2struct(hdr) else h = hdr mid_date=hdr.DATE_OBS mid_time=strmid(mid_date,11,12) exptime=abs(hdr.EXPTIME) filter=hdr.FILTER polar=hdr.P1COL p1col=hdr.P1COL p1row=hdr.P1ROW p2col=hdr.P2COL p2row=hdr.P2ROW ENDIF ENDFOR ; This is experimental but what I'm trying to do here is remove from the datacube any images ; that would not make a good median image. So I start by removing images that are rolled more ; than 5-sigma of the median roll value for that day. Later I will remove pixels that are ; freakishly bright compared to the others. Note that 5-sigma might not be tight enough... we'll see... ;crota_stdv=stddev(hdrs.crota) ;crota_med=median(hdrs.crota) ;good=where( abs(hdrs.crota - crota_med) LT (5*crota_stdv) ) crota_stdv=stddev(crotas) crota_med=median(crotas) good=where( abs(crotas - crota_med) LT (5*crota_stdv) ) ; stop if (n_elements(good) NE n) THEN BEGIN if (n_elements(good) GT 0) THEN datacube=datacube[*,*,good] ELSE BEGIN PRINT,'ERROR! Something went wrong when trying to eliminate rolled images from the datacube...' stop ENDELSE ENDIF ; stop median_array=lonarr(xaxis,yaxis) for j=0,xaxis-1 DO BEGIN for k=0,yaxis-1 DO BEGIN pixel_row = datacube[j,k,*] ; Get rid of zero-value pixels ; stop w = where(pixel_row[0,0,*] GT 0) if (n_elements(w) GT 1 ) THEN good_pix = pixel_row[0,0,*] ELSE good_pix = pixel_row ; Now try to eliminate extreme-value pixels. ; Note that 10-sigma might be too high... ; pix_med = median(pixel_row) ; pix_sd = stddev(pixel_row) ; ok=where( abs(pixel_row - pix_med) LT (10*pix_sd) ) median_array[j,k]=median(good_pix) endfor endfor ======= crota_stdv=stddev(crotas) crota_med=median(crotas) good=where( abs(crotas - crota_med) LT (5*crota_stdv) ) help,good >>>>>>> 1.17 if (n_elements(good) NE n) THEN BEGIN if (n_elements(good) GT 0) THEN datacube=datacube[*,*,good] ELSE BEGIN message,'ERROR! Something went wrong with range of CROTA in selected images...' ENDELSE ENDIF ; stop exptime_stdv =stdev(exptimes[good],exptime_mean) bias_stdv =stdev(biases[good],bias_mean) offsetcr_stdv =stdev(offsetcrs[good],offsetcr_mean) crota_stdv =stdev(crotas[good],crota_mean) nsat =min(nsats) help,exptime_mean dsatval =float(dsatval/exptime_mean) ; normalize this also to 1 second datendtai =anytim2tai(hdr.date_end) date_avg =utc2str(tai2utc(dateobstai + (datendtai-dateobstai)/2)) cadence_stdv =stdev(cadences[good],cadence_mean) readtime_stdv =stdev(readtimes[good],readtime_mean) cleartime_stdv =stdev(cleartimes[good], cleartime_mean) crval1_stdv =stdev(pitches[good], crval1_mean) crval2_stdv =stdev(yaws[good], crval2_mean) cosmics_stdv =stdev(cosmicss[good], cosmics_mean) outhdr0.nmissing=total(nmissings[good]) median_array=fltarr(1024,1024) for j=0,1024-1 DO BEGIN for k=0,1024-1 DO BEGIN pixel_row = datacube[j,k,*] ; Get rid of zero-value pixels ; stop w = where(pixel_row[0,0,*] GT 0) if (n_elements(w) GT 1 ) THEN median_array[j,k]=median(pixel_row[0,0,w]) ; Now try to eliminate extreme-value pixels. ; Note that 10-sigma might be too high... ; pix_med = median(pixel_row) ; pix_sd = stddev(pixel_row) ; ok=where( abs(pixel_row - pix_med) LT (10*pix_sd) ) endfor endfor ; stop window,xsize=512,ysize=512,retain=2 ;if (xaxis GT 1024) THEN median_array=rebin(median_array,1024,1024) maxmin,median_array maxmedian=max(median_array) IF dsatval LT 32000/exptime_mean THEN BEGIN ; COR bscale= 1./exptime_mean median_array_out= fix(round( ((median_array>000,outhdr0,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',nsatload help,dsatval,nsat,bscale get_utc,dte,/ecs fxaddpar,outhdr,'DATE-OBS',date_obs,' Start of first exposure in N_IMAGES' fxaddpar,outhdr,'DATE-END',hdr.date_end,' End of last exposure in N_IMAGES' fxaddpar,outhdr,'DATE-AVG',date_avg,' Midpoint between OBS and END' fxaddpar,outhdr,'CADENCE', median(cadences),' sec (median); stdev='+trim(cadence_stdv) fxaddpar,outhdr,'READTIME',readtime_mean,' sec; stdev='+trim(readtime_stdv) fxaddpar,outhdr,'CLEARTIM',cleartime_mean,' sec; stdev='+trim(cleartime_stdv) fxaddpar,outhdr,'CRVAL1',crval1_mean,' arcsec; stdev='+trim(crval1_stdv) fxaddpar,outhdr,'CRVAL2',crval2_mean,' arcsec; stdev='+trim(crval2_stdv) fxaddpar,outhdr,'COSMICS',cosmics_mean,' stdev='+trim(cosmics_stdv) fxaddpar,outhdr,'DATE',dte fxaddpar,outhdr,'BUNIT','DN/sec' fxaddpar,outhdr,'BSCALE',bscale 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,'CROTA',crota_mean,'deg; =avg over N_IMAGES; stdev='+trim(crota_stdv) fxaddpar,outhdr,'N_IMAGES',n_used,' scc_mk_daily_med.pro' ;stop sdir1 = GETENV('SECCHI_BKG')+'/'+strlowcase(sc)+'/daily_med/' sdir2 = sdir1+strmid(date,0,6) if not file_exist(sdir2) THEN BEGIN PRINT,'Directory ',sdir2, ' does not exist.' cmd='mkdir -p -m 775 '+sdir2 PRINT,cmd spawn,cmd,/sh ENDIF IF keyword_set(SAVEDIR) THEN sdir = savedir ELSE sdir = sdir2 CD,sdir fname0 = 'd'+pref+strupcase(sc)+'_'+pol_str+'_'+strmid(date,2,6)+'.fts' PRINT,'Writing daily fits file: '+sdir+'/'+fname0 print,' using '+string(n_used)+' images' fxaddpar,outhdr,'FILENAME',fname0 fxaddpar,outhdr,'HISTORY',version fxaddpar,outhdr,'COMMENT','%SCC_MK_DAILY_MED: header values from middle image unless noted' WRITEFITS,fname0,median_array_out,outhdr ;stop CD,orig ;IF floor(r) NE last_r THEN GOTO, beginning RETURN END