FUNCTION scc_getbkgimg, inhdr, outhdr=outhdr, silent=silent, totalb=totalb, $ daily=daily, MATCH=match, interpolate=k_interp, ROLL=roll, $ NONLINEARITYCORRECTION=nonlinearitycorrection, _EXTRA=_extra ; ;+ ; $Id: scc_getbkgimg.pro,v 1.39 2008/07/02 19:51:25 nathan Exp $ ; ; Project : STEREO SECCHI ; ; Name : SCC_GETBKGIMG ; ; Purpose : Find and read SECCHI background image ; ; Explanation: This procedure returns the appropriate background image for a ; given SECCHI fits header ; ; Use : IDL> background = scc_getbkgimg(hdr,outhdr=outhdr) ; ; Inputs :hdr: secchi fits header (preferably as a structure) ; ; Outputs : background model image data and (optionally) header from the background file ; ; Keywords : /SILENT = Don't print out informational messages ; /TOTALB = Look for the total brightness background ; /INTERPOLATE = Interpolate between background images. ; /DAILY = Use daily background image instead of monthly. ; OUTHDR = returns the header of the background image ; /MATCH = Match state of bias, exptime, summing, size of output to input hdr ; /ROLL = if difference in CROTA is GT 1 deg, rotate bkg to match input ; /NONLINEARITYCORRECTION = apply COR2 double exposure non-lineariy correction (experimental) ; ; Calls from LASCO : ; ; Common : None. ; ; Restrictions: Need 'SECCHI_BKG' environment variable ; ; Side effects: None. ; ; Category : SECCHI, Calibration ; ; Prev. Hist. : None. ; ; Written : Karl Battams, NRL/I2, APR 2007 ; ; $Log: scc_getbkgimg.pro,v $ ; Revision 1.39 2008/07/02 19:51:25 nathan ; Updated outhdr CRPIX,CDELT if resizing; use CRPIX for rot if /ROLL (Bug 317) ; ; Revision 1.38 2008/06/17 18:29:44 nathan ; added some info messages ; ; Revision 1.37 2008/05/08 20:56:02 thompson ; Check against dates of major spacecraft repointings. ; Fixed bug extrapolating to recent dates. ; ; Revision 1.36 2008/04/15 21:07:01 nathan ; acknowledge /TOTALB in case polar value is incorrect ; ; Revision 1.35 2008/04/14 19:29:25 nathan ; forgot 1 case for closestisbefore var ; ; Revision 1.34 2008/04/14 14:53:06 nathan ; There were problems with the implementation of finding bracketing background ; images that was introduced in rev 1.19. So I reverted to the previous revision and then ; added the features added up to rev 1.33. I also changed some of the variable ; names for consistency through the different parts of the program. ; ; Revision 1.18 2007/10/26 22:46:16 nathan ; Made some changes for new SECCHI_BKG images: looks for files in ; $SECCHI_BKG/../newbkg (NRL only for beta testing--will be removed ; before SSW update); uses DATE-AVG (if available) instead of DATE-OBS ; for image (midpoint) time ; ; Revision 1.17 2007/10/09 20:35:33 thompson ; Removed 2 second wait ; ; Revision 1.16 2007/10/09 19:28:17 thompson ; Extended /SILENT to no-file-found message. ; ; Revision 1.15 2007/10/05 17:44:48 nathan ; Added /MATCH; outhdr.date_obs=inhdr.date_obs if /INTERP set; correct values ; in outhdr which are not currently set correctly in background image FITS ; headers ; ; Revision 1.14 2007/10/04 19:16:32 thompson ; Improved interpolation file selection. ; Changed CONGRID back to REBIN ; ; Revision 1.13 2007/10/03 21:18:01 thompson ; Added keyword /INTERPOLATE ; ; Revision 1.12 2007/09/27 21:11:30 thompson ; Correct NAXIS1/NAXIS2 in output header ; ; Revision 1.11 2007/09/06 15:08:24 nathan ; added _EXTRA; check for seb_prog=DOUBLE (still needs further tweaking to match TBr background to double image ; ; Revision 1.10 2007/09/05 19:00:13 thompson ; bug fix for /daily keyword ; ; Revision 1.9 2007/08/15 15:57:06 reduce ; Tweaked so it should work for mvi headers. Karl B. ; ; Revision 1.8 2007/07/25 21:00:32 thompson ; Added keyword /DAILY ; ; Revision 1.7 2007/07/25 15:49:42 thompson ; Changed REBIN to CONGRID ; ; Revision 1.6 2007/07/10 19:17:51 thompson ; Fixed small typo ; ; Revision 1.5 2007/07/09 20:31:33 thompson ; Rewrote to search by modified julian date ; ; Revision 1.4 2007/07/05 20:13:45 thompson ; Bug fixes. Don't use CD. Call function REDUCE. ; ; Revision 1.3 2007/07/05 19:02:53 reduce ; Took out wait command. Fixed bug with finding 0-deg files ; ; Revision 1.2 2007/07/03 20:27:33 colaninn ; added silent keyword ; ; Revision 1.1 2007/06/21 19:37:59 reduce ; Initial release. Karl B. ; ;- version='$Id: scc_getbkgimg.pro,v 1.39 2008/07/02 19:51:25 nathan Exp $' IF ~keyword_set(SILENT) THEN help,version ; Get the value of the INTERPOLATE keyword. This may be changed below. interp = keyword_set(k_interp) ; Do some verification, get camera, etc del=get_delim() ; get delimiter ;stop IF (datatype(inhdr) NE 'STC') then inhdr=fitshead2struct(inhdr,/DASH2UNDERSCORE) tel=trim(strupcase(inhdr.DETECTOR)) if strupcase(tel) eq 'EUVI' then cam = 'eu' else cam = $ strlowcase(strmid(tel,0,1) + strmid(tel,strlen(tel)-1,1)) ; Added the following stuff so that this routine should work for MVI headers too ; Karl B, 08/15/2007 tags_used=tag_names(inhdr) wOBS=where(tags_used EQ 'OBSRVTRY') IF (wOBS EQ -1) THEN BEGIN ; This must be an mvi header. ; Try to get sc identifier from filename flnm=inhdr.FILENAME sc=strcompress(strmid(filename,strlen(flnm)-5,1),/remove_all) ENDIF ELSE sc = strlowcase(strmid(inhdr.OBSRVTRY,7,1)) ; Make sure the spacecraft name is valid. IF (strupcase(sc) NE 'A') and (strupcase(sc) NE 'B') THEN BEGIN PRINT,'' PRINT,'ERROR!! Could not determine spacecraft from input header' PRINT,'' return, 1 ENDIF if keyword_set(daily) then $ rootdir = concat_dir(getenv('SECCHI_BKG'), sc + del + 'daily_med' + del) else $ rootdir = concat_dir(getenv('SECCHI_BKG'), sc + del + 'monthly_min' + del) ; Form the polar search string, based on the angle in the header, or on the ; keyword TOTALB. polstring='_pTBr_' polar=inhdr.polar IF polar LT 361 and polar GE 0 and ~keyword_set(TOTALB) THEN polstring='_p'+string(round(polar),format='(I3.3)')+'_' ; Get the date from the header, and look for a background file corresponding ; to this date. dt=inhdr.DATE_AVG cal = anytim2cal(dt, form=8, /date) sdir = strmid(cal,0,6) ;Monthly directory name sfil = strmid(cal,2,6) ;Date part of file name if keyword_set(daily) then fchar = 'd' else fchar = 'm' filesrch0 = concat_dir(rootdir, sdir + del + fchar + cam + strupcase(sc) + $ polstring + sfil + '.fts') files = file_search(filesrch0, count=count) delvarx, bkgfile0, bkgfile1 if count gt 0 then begin bkgfile = files[0] if interp then bkgfile0=bkgfile else goto, read_file ; bkgfile0 is same day is input endif ; Otherwise, start looking backwards and forwards in time until a background ; image is found. When /INTERPOLATE is set, look for three images, so that ; the two bracketing the date is set. utcin = anytim2utc(dt) taiin = anytim2tai(utcin) ; Set limits based on the major spacecraft repointings. case strupcase(sc) of 'A': repoint = ['2006-12-21T13:15', '2007-02-03T13:15'] 'B': repoint = ['2007-02-03T18:20', '2007-02-21T20:00'] endcase repoint = anytim2utc(repoint) tai_repoint = utc2tai(repoint) 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 ; Start searching nsearch = 0 while nsearch lt 30 do begin nsearch = nsearch + 1 ; ; Look before the header date. ; utc = utcin utc.mjd = utcin.mjd - nsearch if (utc.mjd gt mjdmin) and (utc.mjd lt mjdmax) then begin cal = anytim2cal(utc, form=8, /date) sdir = strmid(cal,0,6) ;Monthly directory name sfil = strmid(cal,2,6) ;Date part of file name files1 = file_search(concat_dir(rootdir, sdir + del + fchar + cam + $ strupcase(sc) + polstring + sfil + '.fts'), $ count=count1) end else count1 = 0 ; ; And after the header date. ; utc = utcin utc.mjd = utcin.mjd + nsearch if (utc.mjd gt mjdmin) and (utc.mjd lt mjdmax) then begin cal = anytim2cal(utc, form=8, /date) sdir = strmid(cal,0,6) ;Monthly directory name sfil = strmid(cal,2,6) ;Date part of file name files2 = file_search(concat_dir(rootdir, sdir + del + fchar + cam + $ strupcase(sc) + polstring + sfil + '.fts'), $ count=count2) end else count2 = 0 ; ; If files are found for both before and after, then pick the file closest in ; time. If the interpolate keyword is set, and the first background file ; hasn't been picked yet, then take both files. ; if (count1 gt 0) and (count2 gt 0) then begin if interp and (n_elements(bkgfile0) eq 0) then begin bkgfile0 = files1[0] bkgfile = files2[0] closestisbefore=1 goto, read_file endif dummy = sccreadfits(files1[0], h1, /nodata) dummy = sccreadfits(files2[0], h2, /nodata) IF h1.date_avg NE '' THEN h1t=h1.date_avg ELSE h1t=h1.date_obs IF h2.date_avg NE '' THEN h2t=h2.date_avg ELSE h2t=h2.date_obs dtai1 = taiin - anytim2tai(h1t) dtai2 = anytim2tai(h2t) - taiin if dtai2 lt dtai1 then begin bkgfile_found = files2[0] closestisbefore=0 goto, file_found end else begin bkgfile_found = files1[0] closestisbefore=1 goto, file_found endelse ; ; Otherwise, take the file that was found. ; end else if count1 gt 0 then begin bkgfile_found = files1[0] closestisbefore=0 goto, file_found end else if count2 gt 0 then begin bkgfile_found = files2[0] closestisbefore=1 goto, file_found endif ; goto, next_date ; ; Process the found file based on whether one is interpolating or not. ; file_found: bkgfile = bkgfile_found if interp then begin if n_elements(bkgfile0) eq 0 then begin bkgfile0 = bkgfile end else if n_elements(bkgfile1) eq 0 then begin bkgfile1 = bkgfile goto, read_file endif end else goto, read_file ; next_date: endwhile ; If /INTERPOLATE was set, and only one file was found, then proceed as if ; /INTERPOLATE was not set. if interp and (n_elements(bkgfile0) eq 1) then begin IF ~keyword_set(silent) THEN print, $ 'Only one background file found -- not interpolating' interp = 0 bkgfile = bkgfile0 goto, read_file endif ; If no files were found, return an error message. if ~keyword_set(silent) then PRINT, $ '%%SSC_GETBKGIMG: No models found near: ', filesrch0 RETURN, -1 ; Read in the background file(s). If /INTERPOLATE was set, and three files ; were found, then choose the two files which bracket the observation date, ; based on the filename. READ_FILE: if interp then begin IF ~keyword_set(silent) THEN print, 'reading closest to interp: ', bkgfile0 ; This is the closest file to input, could be before OR after bkgimg0 = SCCREADFITS(bkgfile0, outhdr0, silent=silent) IF outhdr0.date_avg NE '' THEN date_avg0=outhdr0.date_avg ELSE date_avg0=outhdr0.date_obs if n_elements(bkgfile1) eq 1 then begin bkgfile2 = bkgfile bkgfile = bkgfile1 break_file, bkgfile0, disk0, dir0, name0 d0 = '20' + strmid(name0,10,6) break_file, bkgfile1, disk1, dir1, name1 d1 = '20' + strmid(name1,10,6) break_file, bkgfile, disk2, dir2, name2 d2 = '20' + strmid(name2,10,6) d = anytim2cal(inhdr.date_obs,form=8,/date) if ((d gt d0) and (d gt d1) and (d lt d2)) or $ ((d lt d0) and (d lt d1) and (d gt d2)) then bkgfile = bkgfile2 endif endif ;IF ~keyword_set(silent) THEN help,bkgfile,bkgfile0,bkgfile1,bkgfile2 IF ~keyword_set(silent) THEN print, 'reading: ', bkgfile bkgimg = SCCREADFITS(bkgfile, outhdr, silent=silent) IF outhdr.date_avg NE '' THEN date_avg=outhdr.date_avg ELSE date_avg=outhdr.date_obs ; ++ Match size. For integer reduction factors, use the ; function REDUCE with /AVERAGE. Otherwise, use REBIN if (inhdr.naxis1 ne outhdr.naxis1) or (inhdr.naxis2 ne outhdr.naxis2) then begin i_reduce = float(outhdr.naxis1) / inhdr.naxis1 j_reduce = float(outhdr.naxis2) / inhdr.naxis2 if (i_reduce eq fix(i_reduce)) and (j_reduce eq fix(j_reduce)) then $ bkgimg = reduce(bkgimg, i_reduce, j_reduce, /average) else $ bkgimg = rebin(temporary(bkgimg),inhdr.naxis1,inhdr.naxis2) outhdr.naxis1 = inhdr.naxis1 outhdr.naxis2 = inhdr.naxis2 outhdr.summed = inhdr.summed ;outhdr.CRPIX1 = 0.5+(outhdr.crpix1-0.5)/i_reduce ;outhdr.CRPIX1A= 0.5+(outhdr.CRPIX1A-0.5)/i_reduce ;outhdr.CRPIX2 = 0.5+(outhdr.crpix2-0.5)/j_reduce ;outhdr.CRPIX2A= 0.5+(outhdr.CRPIX2A-0.5)/j_reduce ;--Because crpix in bkg headers up until May 2008 are incorrect, ; use value from inhdr outhdr.crpix1 = inhdr.crpix1 outhdr.crpix2 = inhdr.crpix2 outhdr.crpix1a = inhdr.crpix1a outhdr.crpix1a = inhdr.crpix1a outhdr.CDELT1 = outhdr.CDELT1*i_reduce outhdr.CDELT2 = outhdr.CDELT2*j_reduce outhdr.CDELT1A = outhdr.CDELT1A*i_reduce outhdr.CDELT2A = outhdr.CDELT2A*j_reduce endif ; If /INTERPOLATE was set, perform the same action on the other background ; image. if interp then begin if (inhdr.naxis1 ne outhdr0.naxis1) or (inhdr.naxis2 ne outhdr0.naxis2) then begin i_reduce = float(outhdr0.naxis1) / inhdr.naxis1 j_reduce = float(outhdr0.naxis2) / inhdr.naxis2 if (i_reduce eq fix(i_reduce)) and (j_reduce eq fix(j_reduce)) then $ bkgimg0 = reduce(bkgimg0, i_reduce, j_reduce, /average) else $ bkgimg0 = rebin(temporary(bkgimg0),inhdr.naxis1,inhdr.naxis2) outhdr0.naxis1 = inhdr.naxis1 outhdr0.naxis2 = inhdr.naxis2 endif ; Interpolate between the two images. if date_avg0 eq date_avg then bkgimg = bkgimg0 else begin tai0 = utc2tai((strsplit(date_avg0,' ',/extract))[0]) tai1 = utc2tai((strsplit(date_avg ,' ',/extract))[0]) factor=(taiin-tai0)/(tai1-tai0) IF ~keyword_set(SILENT) THEN help,date_avg0,date_avg,factor IF ~keyword_set(SILENT) and factor LT 0 THEN print,'Extrapolating...' bkgimg = bkgimg0 + (bkgimg-bkgimg0)*factor endelse outhdr1= outhdr outhdr = outhdr0 ;Header of closest background image outhdr.date_avg=inhdr.date_avg outhdr.obt_time=factor ; for testing IF (closestisbefore) THEN BEGIN outhdr.date_obs=outhdr0.date_obs outhdr.date_end=outhdr1.date_end ENDIF ELSE BEGIN outhdr.date_obs=outhdr1.date_obs outhdr.date_end=outhdr0.date_end ENDELSE endif goto, skiphdrcorrection ; ; -- Correct for bad (early) headers, using input values ; IF outhdr.biasmean LE 0 THEN outhdr.biasmean=inhdr.biasmean IF tel EQ 'COR2' and (outhdr.polar GE 0 and outhdr.polar LE 360) and outhdr.offsetcr LE 0 THEN BEGIN ; these bkgs do not have bias subtracted outhdr.offsetcr=0. ENDIF ELSE BEGIN outhdr.offsetcr=inhdr.biasmean ENDELSE IF tel EQ 'COR2' and outhdr.polar GT 1000 THEN outhdr.exptime=-1. ; TBr combo images used to make bkg IF (tel EQ 'HI1' or tel EQ 'HI2') and inhdr.offsetcr LE 0 THEN inhdr.offsetcr=inhdr.biasmean ; HI images have bias subtracted onboard before summing IF tel EQ 'HI1' THEN outhdr.exptime=1200. IF tel EQ 'HI2' THEN outhdr.exptime=4950. IF outhdr.summed LE 0 THEN BEGIN IF (tel EQ 'HI1' or tel EQ 'HI2') THEN BEGIN outhdr.summed=2. outhdr.ipsum=2 ENDIF IF tel EQ 'COR1' THEN BEGIN outhdr.summed=2. outhdr.ipsum=2 ENDIF IF tel EQ 'COR2' THEN BEGIN outhdr.summed=1. outhdr.ipsum=1 ENDIF ENDIF ; ; --Done header correction ; skiphdrcorrection: ;--Rotate IF difference greater than 1 degree rolldif=inhdr.crota - outhdr.crota IF ~keyword_set(SILENT) THEN help,rolldif IF abs(rolldif) GT 1. and keyword_set(ROLL) THEN BEGIN bc=scc_sun_center(outhdr) ;bkgimg = rot(temporary(bkgimg),rolldif,1.0,bc.xcen,bc.ycen,missing=0,/cubic,/pivot,_extra=ex) bkgimg = rot(temporary(bkgimg),rolldif,1.0,outhdr.crpix1-1,outhdr.crpix2-1,missing=0,/cubic,/pivot,_extra=ex) ;update header outhdr.crota = inhdr.crota outhdr.PC1_1 = inhdr.PC1_1 outhdr.PC1_2 = inhdr.PC1_2 outhdr.PC2_1 = inhdr.PC2_1 outhdr.PC2_2 = inhdr.PC2_2 IF ~keyword_set(SILENT) THEN BEGIN message,'Background rotated '+trim(rolldif)+' deg to match input header.',/info ENDIF ENDIF ELSE IF ~keyword_set(SILENT) THEN message,'Background NOT rotated to match input header.',/info IF keyword_set(MATCH) THEN BEGIN ; ;+++ Modify bkgimg and outhdr to match inhdr ; ; ++ Match binning (not size) sumdif=inhdr.ipsum - outhdr.ipsum IF sumdif NE 0 THEN BEGIN binfac=4^(sumdif) IF ~keyword_set(SILENT) THEN message,'Background corrected for binning * '+trim(binfac),/info bkgimg=bkgimg*binfac ENDIF ; ; Correct for double exposure non-linearity effect. ; if inhdr.seb_prog EQ 'DOUBLE' and keyword_set(NONLINEARITYCORRECTION) THEN begin a0 = 1.04418 a1 = -0.00645004 scl = (a0 + a1*alog(abs(inhdr.exptime))) + a1*alog(bkgimg>1) print,'Computed non-linearity factor (only >1 used):' maxmin,scl bkgimg=bkgimg*(scl>1.) endif ; ++ Match exptime bkgimg=bkgimg*abs(inhdr.exptime/outhdr.exptime) ; ++ Match bias state in OFFSETCR keyword ; secchi_prep removes bias BEFORE ipsum correction ; hi never has bias correction because it is done onboard IF (tel ne 'HI1' and tel ne 'HI2') THEN IF inhdr.offsetcr GT 0 and outhdr.offsetcr LE 0 THEN BEGIN ; subtract bias IF ~keyword_set(silent) THEN print,'subtracting bias...' bkgimg=bkgimg-outhdr.biasmean*4^sumdif outhdr.offsetcr=outhdr.biasmean*4^sumdif ENDIF IF (tel ne 'HI1' and tel ne 'HI2') THEN IF inhdr.offsetcr LE 0 and outhdr.offsetcr GT 0 THEN BEGIN ; add bias back in IF ~keyword_set(silent) THEN print,'adding bias...' bkgimg=bkgimg+outhdr.biasmean*4^sumdif outhdr.offsetcr=outhdr.offsetcr-outhdr.biasmean*4^sumdif ENDIF ;stop ; ++ Update bkg header outhdr.summed=inhdr.summed outhdr.exptime=abs(inhdr.exptime) ENDIF ; ; ++ Match subfield ; ; TBD ; Return the background image. IF ~keyword_set(silent) THEN maxmin,bkgimg return, bkgimg END