FUNCTION scc_getbkgimg, inhdr, outhdr=outhdr, silent=silent, totalb=totalb, $ daily=daily, MATCH=match, interpolate=k_interp, ROLL=roll, $ NONLINEARITYCORRECTION=nonlinearitycorrection, DOMEDIAN=domedian, $ calroll=calroll, raw_calroll=raw_calroll, NOSHIFT=noshift, $ double_totalb=double_totalb, _EXTRA=_extra ; ;+ ; $Id: scc_getbkgimg.pro,v 1.107 2024/12/16 19:37:34 thompson 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 ; /DOUBLE_TOTALB = When used together with /TOTALB, selects double ; exposure backgrounds ; /INTERPOLATE = Interpolate between background images. ; /DAILY = Use daily background image instead of monthly. ; /CALROLL= Use combination of monthly minimum and calibration roll ; backgrounds. These backgrounds may show some problems ; during the long gap between calibration rolls in mid to ; late 2007 ; /RAW_CALROLL = Use raw calibration roll backgrounds. This option ; is not recommended, and is included mainly to support ; generating the combined monthly/calroll backgrounds. ; 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) ; DOMEDIAN= Apply median filter of thix box size (>3) ; /NOSHIFT Do not use hi_align_image to match output with input (HI only) ; ; 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.107 2024/12/16 19:37:34 thompson ; Added COR1 dust events ; ; Revision 1.105 2023/08/15 21:29:16 thompson ; Added rectification change on 2023-08-12 ; ; Revision 1.104 2023/05/24 18:55:30 thompson ; Added 2021-05-07 event ; ; Revision 1.103 2021/05/10 20:25:54 thompson ; Added 2021-05-07 event ; ; Revision 1.102 2020/01/22 16:42:05 thompson ; *** empty log message *** ; ; Revision 1.101 2019/02/20 21:30:21 thompson ; Added COR1-A event on 2019-02-17 ; ; Revision 1.100 2019/02/19 22:03:16 thompson ; Added COR1-A event on 2019-02-05 ; ; Revision 1.99 2018/10/04 20:10:42 thompson ; Added COR1-A event on 2018-09-27 ; ; Revision 1.98 2017/08/10 21:49:49 nathan ; fix roll check for post-conjunction ; ; Revision 1.97 2016/09/07 15:59:38 thompson ; Added COR1 change on 2 Sep 2016 ; ; Revision 1.96 2016/06/20 15:45:55 thompson ; Add exposure time change back on June 8, 2016 ; ; Revision 1.94 2016/04/19 14:04:39 thompson ; Added COR1A event on 2016-04-13 ; ; Revision 1.93 2016/04/04 15:56:53 thompson ; Add change of COR1 exposure time to 0.7 seconds ; ; Revision 1.92 2015/11/23 22:07:54 thompson ; Added resumption of normal operations after 2015-Nov-16 as a break point for COR1 ; ; Revision 1.91 2014/09/03 21:55:55 thompson ; Changed so that 2014-08-19 repoint on Ahead doesn't apply to COR1 ; ; Revision 1.90 2014/09/03 17:37:09 secchia ; added a repoint day for SC_A for start of sidelobe 1 ; ; Revision 1.89 2014/09/02 15:41:15 thompson ; Added COR1 event of 2014-08-23 ; ; Revision 1.88 2014/07/17 20:13:29 thompson ; Include S/C testing July 6-11 for COR1-A ; ; Revision 1.87 2013/11/21 22:19:09 nathan ; fix level-1 test for not-hi case ; ; Revision 1.86 2013/10/28 18:48:41 nathan ; several adjustments to properly do hi_prep on bkg before returning ; ; Revision 1.85 2013/03/26 22:53:42 nathan ; re- remove doshift ; ; Revision 1.84 2013/03/22 18:23:44 nathan ; implement doshift and check tel for cor2_warp ; ; Revision 1.83 2012/05/24 15:29:07 nathan ; Undo default roll correction (pylon is a problem) ; ; Revision 1.82 2012/05/24 15:07:21 nathan ; By default roll bkg if difference gt 10 deg; set outhdr.jitrsdev based ; on maxshift; run fix_edges() if /DOSHIFT ; ; Revision 1.81 2012/05/23 15:20:22 secchib ; nr - change max roll for default behavior from 8 to 10 ; ; Revision 1.80 2012/03/06 16:14:43 thompson ; Added COR1-A event of 2012-02-19 ; ; Revision 1.79 2012-02-14 22:47:27 nathan ; make sure doshift is disabled ; ; Revision 1.78 2012/02/14 22:43:53 nathan ; update both input headers if /interp ; ; Revision 1.77 2012/02/02 20:13:18 nathan ; update dstop in outhdr ; ; Revision 1.76 2012/01/13 20:49:04 thompson ; Added keyword /DOUBLE_TOTALB ; ; Revision 1.75 2011-12-16 20:08:46 thompson ; Added COR1-A event of 2011-12-05. ; Only activate wait statements if DEBUGON is set. ; All these waits are really causing a problem. ; ; Revision 1.74 2011-12-15 17:54:54 nathan ; if /daily do not interpolate if input time not between found files ; ; Revision 1.73 2011/11/21 21:17:14 nathan ; Add /DOMEDIAN option; add hi_align_image code but disable; ; add code for HI interp special case but disable; do not ; check for rolled images in early mission ; ; Revision 1.72 2011/08/15 21:27:29 nathan ; expand postd (roll) search ; ; Revision 1.71 2011/08/10 18:21:30 nathan ; do not wait 2 sec if no double bkg found ; ; Revision 1.70 2011/08/10 18:16:41 nathan ; do not get special rolled bkg for COR1 ; ; Revision 1.69 2011/08/08 17:41:49 nathan ; never match roll for HI ; ; Revision 1.68 2011/08/04 22:36:27 nathan ; implement postd if sc_roll GT 8 ; ; Revision 1.67 2011/07/21 21:52:36 nathan ; for double images, go to TBr if dbTB not found ; ; Revision 1.66 2011/07/21 20:10:07 secchia ; do not go to TBr if dbl not available ; ; Revision 1.65 2011/07/13 19:32:00 thompson ; Remove restriction on use of 'tbr' with COR1 ; ; Revision 1.64 2011-05-09 23:10:34 secchib ; nr - re-commit new straylight pattern date for cor2-b ; ; Revision 1.63 2011/04/29 20:19:43 secchib ; nr - remove new stray light date until we have a new bkg ; ; Revision 1.62 2011/04/27 17:43:27 nathan ; add cor2-b event ; ; Revision 1.61 2011/03/21 17:45:36 thompson ; Added COR1-A event on 2011-03-08, and COR1-B event on 2011-03-11 ; ; Revision 1.60 2011/02/28 16:00:41 thompson ; Added COR1-A event on 2011-02-11 ; ; Revision 1.59 2011-02-15 21:50:28 nathan ; Handle all changes in background the same way, as dates in array pointing; ; utilize debugon flag for not-silent ; ; Revision 1.58 2011/02/04 20:29:22 nathan ; stray light change special case ; ; Revision 1.57 2011/02/02 16:57:21 secchib ; LM - set factor if bkg the same ; ; Revision 1.56 2011/01/28 16:30:44 thompson ; Added COR1-A event on 2011-01-12 ; ; Revision 1.55 2011/01/19 20:57:04 nathan ; treat HI islevelone differently ; ; Revision 1.54 2010/12/02 15:59:18 thompson ; Added COR1-A change on 2010-11-19 ; ; Revision 1.53 2010/11/19 14:47:29 mcnutt ; will look for double background if DOUBLE if not found will look for TB backgrounds ; ; Revision 1.52 2010/08/06 20:02:41 nathan ; fixed problem with different size bkg introduced with previous rev ; ; Revision 1.51 2010/05/03 18:46:01 nathan ; Check filename for level-1 and process accordingly; this necessitated moving ; a few things around. ; ; Revision 1.50 2010/04/09 16:27:51 thompson ; Add COR1-B event of 2010-03-24. ; ; Revision 1.49 2010/02/12 23:17:57 thompson ; Added COR1A event of 2010-01-27 ; ; Revision 1.48 2009/05/14 16:29:44 thompson ; Treat COR1 change from 1024x1024 to 512x512 as a repoint, except for raw calroll ; data. ; ; Revision 1.47 2009/03/26 16:23:34 nathan ; oops ; ; Revision 1.46 2009/03/26 15:41:00 nathan ; give a notice and exit early if cor1 tbr is requested ; ; Revision 1.45 2009/02/11 20:07:12 thompson ; Add COR1 change on 2009-01-30 ; ; Revision 1.44 2008-11-20 14:01:27 mcnutt ; if match checks distcorr and applies cor2_warp if different ; ; Revision 1.43 2008/09/26 19:29:18 nathan ; set ipsum in outhdr ; ; Revision 1.42 2008/09/18 21:50:05 nathan ; Fixed incorrect hdr when /interp; set nmissing=0 ; ; Revision 1.41 2008/08/28 16:28:32 nathan ; workaround for incorrect polar val in DOUBLE header (bug 222) ; ; Revision 1.40 2008/08/26 21:28:37 thompson ; Added keywords /calroll and /raw_calroll ; ; 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.107 2024/12/16 19:37:34 thompson Exp $' IF keyword_set(SILENT) THEN debugon=0 ELSE debugon=1 IF debugon 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 tel eq 'EUVI' then $ cam = 'eu' else $ cam = strlowcase(strmid(tel,0,1) + strmid(tel,strlen(tel)-1,1)) ishi = (cam EQ 'h1') or (cam EQ 'h2') doshift=0 ;IF ishi and not keyword_set(NOSHIFT) THEN doshift=1 ELSE doshift=0 ; 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') filename=inhdr.FILENAME IF (wOBS EQ -1) THEN BEGIN ; This must be an mvi header. ; Try to get sc identifier from filename sc=strlowcase(rstrmid(filename,4,1)) ENDIF ELSE sc = strlowcase(strmid(inhdr.OBSRVTRY,7,1)) ; Make sure the spacecraft name is valid. IF sc NE 'a' and sc NE 'b' THEN BEGIN PRINT,'' PRINT,'ERROR!! Could not determine spacecraft from input header' PRINT,'' return, 1 ENDIF ;--Check for processing level islevelone=0 levchar=strmid(filename,16,1) IF levchar EQ '1' or (levchar EQ '0' and ishi) THEN BEGIN islevelone=1 match=0 IF debugon THEN print,'Level-1 image detected.' ENDIF if keyword_set(raw_calroll) then begin ndays = 365 rootdir = concat_dir(getenv('SECCHI_BKG'), sc + del + 'roll_min' + del) fchar = 'r' end else if keyword_set(calroll) then begin ndays = 30 rootdir = concat_dir(getenv('SECCHI_BKG'), sc + del + 'monthly_roll' + del) fchar = 'mr' end else if keyword_set(daily) then begin ndays = 30 rootdir = concat_dir(getenv('SECCHI_BKG'), sc + del + 'daily_med' + del) fchar = 'd' end else begin ndays = 30 rootdir = concat_dir(getenv('SECCHI_BKG'), sc + del + 'monthly_min' + del) fchar = 'm' endelse ; Form the polar search string, based on the angle in the header, or on the ; keyword TOTALB. polstring='_pTBr_' if keyword_set(TOTALB) and keyword_set(DOUBLE_TOTALB) then polstring='_dbTB_' polar=inhdr.polar IF polar LT 361 and polar GE 0 and ~keyword_set(TOTALB) and inhdr.seb_prog NE 'DOUBLE' $ THEN polstring='_p'+string(round(polar),format='(I3.3)')+'_' IF polar GT 361 and ~keyword_set(TOTALB) and inhdr.seb_prog EQ 'DOUBLE' THEN polstring='_dbTB_' ;;if debugon and polstring EQ '_pTBr_' and cam EQ 'c1' then BEGIN ;; print ;; PRINT, '%%SSC_GETBKGIMG: There are no TB background images made for COR1.' ;; print, 'The default for SECCHI_PREP is to subtract a background for each' ;; print, 'polarization angle before combining into TB. Returning.' ;; print ;; RETURN, -1 ;;ENDIF ; Get the date from the header, and look for a background file corresponding ; to this date. dtin=inhdr.DATE_AVG cal = anytim2cal(dtin, form=8, /date) sdir = strmid(cal,0,6) ;Monthly directory name sfil = strmid(cal,2,6) ;Date part of file name utcin = anytim2utc(dtin) taiin = anytim2tai(utcin) ; Set limits based on the major spacecraft repointings, or other major events. ; ; WARNING!!! Do not input date of event until there is a background available after that date!!! ; case strupcase(sc) of 'A': begin repoint = ['2006-12-21T13:15', '2007-02-03T13:15','2015-05-19', '2023-08-12'] if strupcase(tel) eq 'COR1' then repoint = $ [repoint, '2010-01-27T16:49', '2010-11-19T16:00', $ '2011-01-12T12:23', '2011-02-11T04:23', '2011-03-08T17', $ '2011-12-05T12:03', '2012-02-19T02:33', '2012-03-16', '2014-07-08', $ '2014-08-23T17', '2015-11-16', '2016-03-15', '2016-04-13T01:08', $ '2016-05-25T05:40:05', '2016-06-01T12', '2016-06-08', $ '2016-09-02T16:19', '2018-09-27T16:44', '2019-02-05T21:07', $ '2019-02-17T05:57', '2020-01-14T23:52', '2021-05-07T23:53', '2023-05-15T23', $ '2023-08-31T19:59', '2024-11-29T09:55'] if strupcase(tel) ne 'COR1' then repoint = [repoint, '2014-08-19T00:00'] nomrollmjd=54160 ; Early mission roll: A=2007/03/01, B=2007/08/01 endcase 'B': begin repoint = ['2007-02-03T18:20', '2007-02-21T20:00'] if strupcase(tel) eq 'COR1' then repoint = $ [repoint, '2009-01-30T16:20', '2010-03-24T01:17', '2011-03-11'] if strupcase(tel) eq 'COR2' then repoint = $ [repoint, '2010-02-23T08:12', '2011-01-27T03:47','2011-04-25t18:30'] nomrollmjd=54313 endcase endcase ; Treat change to 512x512 as repoint, except for raw calibration rolls. This ; exception is made to allow calibration roll data to be used across the ; boundary. if (strupcase(tel) eq 'COR1') and not keyword_set(raw_calroll) then begin repoint = [repoint, '2009-04-18T23:59:59'] s = sort(repoint) repoint = repoint[s] endif 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 pntroll=inhdr.sc_roll IF (utcin.mjd GT 57160) and (utcin.mjd lt 60168) THEN pntroll=inhdr.crota ; post-conjunction IF pntroll LT 0 THEN titleangle = 360+pntroll ELSE titleangle = pntroll help,pntroll ; IF abs(pntroll) LT 10 or tel EQ 'COR1' or utcin.mjd LT nomrollmjd THEN postd='' $ ; COR1 bkg does not change much with roll ELSE BEGIN postd = 'r'+strmid(STRING(round(titleangle),format='(I3.3)'),0,1)+'*' ENDELSE IF ishi THEN roll=0 ; days when stray light pattern changed are a special case. ; bkg for day of change will be new bkg hhmm=strmid(filename,9,4) ;IF sfil EQ '110127' and sc EQ 'b' THEN IF hhmm LT '0353' THEN sfil='110126' ELSE interp=0 onemoretime: filesrch0 = concat_dir(rootdir, sdir + del + fchar + cam + strupcase(sc) + $ polstring + sfil +postd+ '.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. ; Start searching nsearch = 0 IF debugon THEN help,dtin while nsearch lt ndays do begin nsearch = nsearch + 1 fn1='--------------------' fn2='--------------------' ; ; Look before the header date. ; utc = utcin utc.mjd = utcin.mjd - nsearch ymd1=utc2yymmdd(utc) 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 fn1=fchar + cam + strupcase(sc) + polstring + sfil +postd + '.fts' files1 = file_search(concat_dir(rootdir, sdir + del + fn1), count=count1) end else count1 = 0 ; ; And after the header date. ; utc = utcin utc.mjd = utcin.mjd + nsearch ymd2=utc2yymmdd(utc) 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 fn2=fchar + cam + strupcase(sc) + polstring + sfil +postd+ '.fts' files2 = file_search(concat_dir(rootdir, sdir + del + fn2), 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 debugon THEN print,nsearch,' ',ymd1,' ',ymd2,byte(count1),byte(count2),' ',fn1,' ',fn2 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 postd NE '' THEN BEGIN message,'No '+postd+' models found within '+trim(nsearch)+' days of '+dtin,/info message,'Trying unrolled instead...',/info postd='' IF debugon THEN wait,2 goto, onemoretime endif if polstring EQ '_dbTB_' then begin ; look for cor2 pTBr if double is not found message,'No '+polstring+' '+postd+' models found within '+trim(nsearch)+' days of '+dtin,/info message,'Trying _pTBr_ instead...',/info polstring='_pTBr_' IF debugon THEN wait,2 goto, onemoretime endif ; If no files were found, return an error message. if ~keyword_set(silent) then PRINT, $ '%%SSC_GETBKGIMG: No models found near: ', filesrch0 if debugon then wait,10 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 keyword_set(DAILY) and interp THEN BEGIN ; Do not do interpolation with daily images if intime is not between the two found files x=sccreadfits(bkgfile0,h0,/nodata) x=sccreadfits(bkgfile,h1,/nodata) tai0=anytim2tai(h0.date_avg) tai1=anytim2tai(h1.date_avg) IF tai0 GT tai1 THEN BEGIN t0=tai1 & t1=tai0 ENDIF ELSE BEGIN t0=tai0 & t1=tai1 ENDELSE IF taiin LT t0 or taiin GT t1 THEN BEGIN IF ~keyword_set(silent) THEN BEGIN message, 'Input time not between 2 files -- not interpolating',/info if debugon then wait,2 ENDIF interp = 0 bkgfile = bkgfile0 ENDIF ENDIF if interp then begin IF ~keyword_set(silent) THEN print, 'reading closest to interp: ' IF ~keyword_set(silent) THEN print, bkgfile0 ; This is the closest file to input, could be before OR after IF (islevelone) and not ishi THEN $ ; HI is special case SECCHI_PREP,bkgfile0,outhdr0,bkgimg0,WARP_OFF=(inhdr.DISTCORR EQ 'F'), OUTSIZE=2048/(2^(inhdr.summed-1)) ELSE $ bkgimg0 = SCCREADFITS(bkgfile0, outhdr0, silent=silent) IF keyword_set(DOMEDIAN) THEN bkgimg0=median(temporary(bkgimg0),(domedian>3)) 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: ' IF ~keyword_set(silent) THEN print, bkgfile IF (islevelone) and not ishi THEN $ SECCHI_PREP,bkgfile,outhdr,bkgimg,WARP_OFF=(inhdr.DISTCORR EQ 'F'), OUTSIZE=2048/(2^(inhdr.summed-1)) ELSE $ bkgimg = SCCREADFITS(bkgfile, outhdr, silent=silent) IF keyword_set(DOMEDIAN) THEN bkgimg=median(temporary(bkgimg),(domedian>3)) 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 ; outhdr.summed set later i_reduce = float(outhdr.naxis1) / inhdr.naxis1 j_reduce = float(outhdr.naxis2) / inhdr.naxis2 if (inhdr.naxis1 ne outhdr.naxis1) or (inhdr.naxis2 ne outhdr.naxis2) then begin 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) endif IF ~(islevelone) THEN BEGIN ; Why is this ~islevelone? outhdr.dstop1 = inhdr.dstop1 outhdr.dstop2 = inhdr.dstop2 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.crpix2a = inhdr.crpix2a 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 maxshift=0. IF doshift THEN bkgimg=hi_align_image(outhdr,bkgimg,inhdr,/zeromiss,/verbos) maxshift=maxshift>outhdr.jitrsdev ; match bkgimg to inhdr AFTER matching size ; HI special case: if there are not enough medians in the bkg, then there are a lot of ; artifacts from stars. So rather than interp, use the bkg with most days in it. if 0 then begin ;IF ishi and interp THEN BEGIN n_images=outhdr.n_images n_images0=outhdr0.n_images usefile='' IF outhdr.n_images LT 20 and outhdr0.n_images GT 20 THEN BEGIN outhdr=outhdr0 bkgimg=bkgimg0 usefile=outhdr0.filename ENDIF IF outhdr.n_images GT 20 and outhdr0.n_images LT 20 THEN BEGIN outhdr0=outhdr bkgimg0=bkgimg usefile=outhdr.filename ENDIF IF usefile ne '' and ~keyword_set(SILENT) then BEGIN help,n_images,n_images0 message,'Not interpolating, only using '+usefile,/info ENDIF 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) endif IF ~(islevelone) THEN BEGIN ; necessary because wcs hdr used in hi_align_image below outhdr0.dstop1 = inhdr.dstop1 outhdr0.dstop2 = inhdr.dstop2 outhdr0.naxis1 = inhdr.naxis1 outhdr0.naxis2 = inhdr.naxis2 outhdr0.summed = inhdr.summed ;outhdr0.CRPIX1 = 0.5+(outhdr0.crpix1-0.5)/i_reduce ;outhdr0.CRPIX1A= 0.5+(outhdr0.CRPIX1A-0.5)/i_reduce ;outhdr0.CRPIX2 = 0.5+(outhdr0.crpix2-0.5)/j_reduce ;outhdr0.CRPIX2A= 0.5+(outhdr0.CRPIX2A-0.5)/j_reduce ;--Because crpix in bkg headers up until May 2008 are incorrect, ; use value from inhdr outhdr0.crpix1 = inhdr.crpix1 outhdr0.crpix2 = inhdr.crpix2 outhdr0.crpix1a = inhdr.crpix1a outhdr0.crpix2a = inhdr.crpix2a outhdr0.CDELT1 = outhdr0.CDELT1*i_reduce outhdr0.CDELT2 = outhdr0.CDELT2*j_reduce outhdr0.CDELT1A = outhdr0.CDELT1A*i_reduce outhdr0.CDELT2A = outhdr0.CDELT2A*j_reduce ENDIF IF doshift THEN bkgimg0=hi_align_image(outhdr0,bkgimg0,inhdr,/zeromiss,/verbos) maxshift=maxshift>(outhdr.jitrsdev) outhdr.jitrsdev=maxshift ; Interpolate between the two images. if date_avg0 eq date_avg then begin bkgimg = bkgimg0 factor=1.0 endif 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 IF doshift THEN BEGIN bi0=bkgimg bkgimg=fix_edges(bi0,ceil(maxshift)+2,0.1) ENDIF ;--All of this stuff gets set in SECCHI_PREP: ;if match check distcorr for bkgimg IF keyword_set(MATCH) and inhdr.DISTCORR eq 'T' and outhdr.DISTCORR eq 'F' and cam EQ 'c2' THEN $ bkgimg = cor2_warp(temporary(bkgimg),outhdr, INFO=histinfo, _EXTRA=ex) ;--------------------------------------------------------------- outhdr.nmissing=0 ; nr, 080918 - if 0, breaks hi_prep.pro 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 (ishi) and inhdr.offsetcr LE 0 THEN inhdr.offsetcr=inhdr.biasmean ; HI images from summing buffers 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 (ishi) 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,domedian 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 ; ; Correct for double exposure non-linearity effect. This is NOT done in SECCHI_PREP ; 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 IF keyword_set(MATCH) THEN BEGIN ;--All these are done in SECCHI_PREP (match set to zero above) ; ; ++ 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 outhdr.ipsum=inhdr.ipsum ENDIF ; ++ Match exptime bkgimg=bkgimg*abs(inhdr.exptime/outhdr.exptime) ; ++ Match bias state in OFFSETCR keyword ; HI has bias correction onboard if it went to HI summing buffer 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 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 ; ; ++ Do HI level-1 correction ; IF ishi and islevelone THEN BEGIN ; NOTE: this will not work if input header had /UPDATE_HDR_OFF unless: ; IF used /NOCALFAC, then islevelone is set. IF strlen(grep('Flat Field',inhdr.history)) GT 1 THEN nocalimg=0 ELSE nocalimg=1 IF strlen(grep('calibration factor',inhdr.history)) GT 1 THEN nocalfac=0 ELSE nocalfac=1 IF strlen(grep('exposure weighting',inhdr.history)) GT 1 THEN nodesmear=1 ELSE nodesmear=0 noexpcorr=1 IF strlen(grep('desmearing method',inhdr.history)) GT 1 $ or strlen(grep('Applied hi_desmear.pro',inhdr.history)) GT 1 $ or (nodesmear) THEN noexpcorr=0 corrforipsum=0 help,nocalimg,noexpcorr,nocalfac,nodesmear bk0=bkgimg inhdr1=inhdr IF stregex( inhdr.IP_00_19,' 37',/boolean) $ ; Use HI-1 summing buffer or stregex( inhdr.IP_00_19,' 38',/boolean) $ ; Use HI-2 summing buffer THEN BEGIN ; put values at same level as raw image for hi_prep inhdr1.ipsum=2. IF (nocalfac) THEN corrforipsum=1 ENDIF inhdr1.bunit='DN' bkgimg=bkgimg*inhdr.exptime*(inhdr1.ipsum)^2 hi_prep,inhdr1,bkgimg, saturation_limit=16384*4L, calimg_off=nocalimg, DESMEAR_OFF=nodesmear, $ EXPTIME_OFF=noexpcorr, SILENT=quiet, calfac_off=nocalfac IF (corrforipsum) THEN bkgimg=bkgimg/4. ENDIF ; ; ++ Match subfield ; ; TBD ; Return the background image. IF ~keyword_set(silent) THEN maxmin,bkgimg return, bkgimg END