;+ ; NAME: ; smei_star_remove ; PURPOSE: ; Removes bright stars from SMEI orbital sky maps ; CATEGORY: ; ucsd/camera/idl/star ; CALLING SEQUENCE: PRO smei_star_remove, wanted_map, $ destination = destination , $ to_smeidb = to_smeidb , $ checkversion = checkversion , $ overwrite = overwrite , $ force = force , $ silent = silent , $ keepbkgnd = keepbkgnd , $ fit_with_nobkgnd = fit_with_nobkgnd , $ cleanedge = cleanedge , $ version = version , $ get_version = get_version , $ reverse_loop = reverse_loop , $ use_bkgnd_version = use_bkgnd_version , $ use_zld_version = use_zld_version , $ _extra = _extra ; INPUTS: ; wanted_map array[n]; type: string, time structure or integer ; passed to href=smei_getfile= ; Selects skymap files to be processed ; OPTIONAL INPUT PARAMETERS: ; ; /to_smeidb if set, write to directory structure as used ; for SMEI data base (separate subdirs for mode ; and camera) with 'destination' as the base dir. ; destination scalar; type: string; default: $TUB ; array[2]; type: string ; destination directory for skymaps and time series file ; if /to_smeidb NOT set then maps are written ; directly to 'destination'. ; The default in this case is $TUB ; if /to_smeidb SET then files are written to ; a directory structure as used for the ; SMEI data base with 'destination' as the ; base directory. ; The default in this case is 'SMEIDB?' (the ; SMEI database). ; If two directories are specified the first is used ; for the skymap and the second for the time series ; If destination includes a "user@host:" identifier then ; the output is created locally in $TUB ; and is then transferred to the remote host using scp; ; the local files are deleted after the transfer; ; note that this will only work if the remote host has ; been configured for ssh using private/public keys ; allowing login without passwords. ; /overwrite forces overwriting of existing skymaps and timeseries files ; (by default no files are overwritten) ; /checkversion if present, files will be overwritten only if the current ; star subtraction version number and the current sidereal ; background version is higher then the version in the Fits ; header of the skymap. ; (if /keepbkgnd is set then the background version is not ; checked) ; ; silent=silent if set, suppresses informational messages ; upto_mags=upto_mags ; scalar; type: float; default: 6.0 ; maximum SMEI magnitude; only stars brighter than this ; are removed. (passed to href=smei_star_fit) ; /keepbkgnd if set, the sidereal background (Milky Way, faint stars, ; etc.) is not subtracted from the skymaps after ; the bright stars are subtracted. ; ; /cleanedge if set, a correction is made for straylight from sky ; just outside the fov into the outer areas (direction ; cosine angles larger than 20 deg) of the camera fov. ; /fit_with_nobkgnd if set, stars are fit with the sidereal background removed ; ; _extra=_extra The _extra keyword is passed to href=smei_star_fit=, ; and href=smei_getfile=, so check those for ; additional input keywords. ; ; /get_version retrieve software version number ; (returned in keyword "version") ; use_bkgnd_version=use_bkgnd_version ; scalar; type: integer 1,2 or 3; default: 3 ; selects the sidereal background map ; use_zld_version=use_zld_version ; scalar; type: integer 1,2,3,4; default: 4 ; selects the zodiacal brightness model ; (passed to smei_zld_model) ; OUTPUTS: ; Files in destination directory ; OPTIONAL OUTPUTS: ; version=version Only defined if /get_version is set ; array[2]; type: float ; version[0]: version of bright star subtraction ; version[1]: version of sidereal background map ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, CheckDir, IsType, destroyvar ; smei_filename, smei_filepath, smei_getfile, smei_star_fit ; smei_camera_gain, smei_star_writepnt, smei_sky_cleanedge_fov ; GetFileSpec, TimeSystem, TimeUnit, gzip_file ; hide_env ; EXAMPLE: ; smei_star_remove, 1500 ; Subtract stars for orbit 1500 (specified as integer) for camera 1 ; smei_star_remove, smei_coriolis([1500,1501]), camera=1 ; Subtract stars for orbits 1500 and 1501 (specified as time ; structure) for camera 1 ; smei_star_remove, '2003/04/16 18:25:57' ; Subtract stars for orbit 1500 specified as time string ; smei_star_remove, '*.fts.gz' ; Use dialog_pickfile to select single skymap ; smei_star_remove, 'c1sky_2006_045_021031.fts.gz' ; Subtract stars for specified map ; PROCEDURE: ; For each skymap a new skymap (Fits file) with the star-subtracted ; maps is created, and an ASCII file with information about the ; fitted stars. Both files are created in the directories ; specified in keyword 'destination'. ; ; The current star subtraction is build around the following ; set of keywords: ; ; smei_star_remove, f, /cleanedge, /use_weights, /auto_wing ; ; This implies a wind_radius value of 1.4 degrees. ; A good way to redo a section of maps in a time range 'trange' in the SMEI ; skymap data base on $SMEISKY0 is the command: ; ; smei_star_remove, trange,cam=3,/cleanedge,/use_weights,/auto_wing, $ ; /checkversion,silent=3,/to_smeidb ; ; MODIFICATION HISTORY: ; NOV-2005, Paul Hick (UCSD/CASS) ; JUN-2006, Paul Hick (UCSD/CASS) ; Modified to accept arrays as input in wanted_map argument. ; JUN-2006, Paul Hick (UCSD/CASS) ; Added Fits keyword VERSION to sky map ; Added NAME, CREATED and VERSION keywords to time series file ; Added keyword /overwrite ; JUL-2006, Jordan Vaughan, Paul Hick (UCSD/CASS) V2.0 ; Added subtraction of sidereal background map (can be disabled ; by keyword /keepbkgnd) ; ; Added two extensions (copied from original skymap): ; -extension 6: fraction of orbit (extension 7 from original skymap) ; -extension 7: time in sec since first frame used ; (extension 14 from original skymap) ; ; Keyword VERSION is now called BSTARVER (to distinguish it from ; keyword BKGNDVER (the version of the sidereal background map) ; Keyword MILKYWAY is now called BKGND ; If the sidereal backgound is subtracted keyword BKGND is set ; to 1, and two more keywords are added: ; BKGNDVER is the version of the sidereal background map ; BKGAIN is the gain correction (multiplicative factor) applied ; to the sidereal background map before subtracting it. ; ; Fixed problem with header for extension 4 (map of removed stars ; in equatorial map. ; AUG-2006, Paul Hick (UCSD/CASS) V3.0 ; Modifications related to separate fitting of background and PSF. ; Added couple of entries to smei_star_fit structure. ; Time series file is now written by procedure smei_star_writepnt ; AUG-2006, Paul Hick (UCSD/CASS) V3.0 ; Padded records in time series file with 5 spaces to make them ; exactly the same length as the titles record. ; SEP-2006, Paul Hick (UCSD/CASS) V3.1 ; Added /use_filled_psf keyword ; Added star_info keyword to smei_star_fit and smei_star_writepnt ; NOV-2006, Paul Hick (UCSD/CASS) V4.0 ; Rearranged layout of time series file ; Added extension (lowres map copied from original skymap) ; -extension 8: number of pixels per skybin ; (extension 16 from original skymap) ; MAR-2007, Paul Hick (UCSD/CASS) V4.1 ; Added extension (lowres map copied from original skymap) ; -extension 8: direction cosine angle ; (extension 15 from original skymap) ; -extension 9: old extension 8 (pixels per skybin) ; APR-2007, Paul Hick (UCSD/CASS) V4.2 ; Added cleanedge keyword to allow removal of straylight ; from sky just outside the camera fov. ; MAY-2007, Paul Hick (UCSD/CASS) V4.3 ; Modified fovangle field and added dfovangle field ; to star_fit structure ; MAY-2007, Paul Hick (UCSD/CASS) V4.31 ; If keywords /cleanedge is set and /keepbkgnd NOT set, ; then the sidereal background is added again AFTER ; the star subtraction, and PRIOR to the removal of ; the scattered light along the edge of the FOV. After ; the scattered light correction the sidereal background ; is taken out again. ; Added keyword /to_smeidb to write directly into SMEI data base ; Added keyword remote to do transfers to remote host ; JUN-2007, Paul Hick (UCSD/CASS) ; Added keyword /from_smeidb ; JUL-2007, Paul Hick (UCSD/CASS) ; Bug fix. removebkgnd was written into Fits file as byte ; instead of short int. ; JUL-2007, Paul Hick (UCSD/CASS), V5.01 ; Added three extensions to output Fits file. ; These contain 3 skymaps (equatorial,north and south pole) ; with the stars subtracted and the core PSF replaced by ; the background model. ; Modified use of area defined with wing_radius. Bins in this ; area that are part of the core PSF of a brighter neighbour ; star (which has been subtracted already) are excluded from ; being used in the lsq fit. They are still used when the ; star model is subtracted. ; Removed keywords wing_radius and use_filled_psf ; (now passed to smei_star_fit using _extra keyword). ; NOV-2007, Paul Hick (UCSD/CASS), V5.02 ; Fixes bug in calculation of camera gain for cameras 1 and 3 ; smei_camera_gain. ; DEC-2007, Paul Hick (UCSD/CASS), V5.03 ; Added /inside_range keyword to be able to process ; files in a specified time range (specified as a 2-element ; time array in 'wanted_map') ; Changed keyword CREATED in time series file to DATE. ; Shortened header for equatorial maps in ext 3 and 6 ; (TX1MOON, TX2MOON, TX3MOON removed). ; DEC-2007, Paul Hick (UCSD/CASS), V5.04 ; Added two lines to header for timeseries files ; (clean edge and remove bkgnd keywords) ; DEC-2007, Paul Hick (UCSD/CASS), V5.05 ; Added some more keyword to timeseries header: ; BSTAR,BSTARVER,BKGNDVER,BKGAIN,CLNEDGE ; (these are also present in the Fits file with the ; star-subtracted skymaps). ; DEC-2007, Paul Hick (UCSD/CASS), V5.06 ; Added keywords ONEORBIT and ORBIT to timeseries header. ; JAN-2008, Paul Hick (UCSD/CASS), V5.07 ; Fixed some format issues in pnt file. ; Star name was truncated at 11 char, should be 12. ; Planet names were print right-adjusted, should be ; left-adjusted ; FEB-2008, Paul Hick (UCSD/CASS) ; Location of files now done using smei_getfile. ; Removed keyword remote (functionally is now provided by ; keyword destination). ; FEB-2008, Paul Hick (UCSD/CASS), V5.10 ; Added star_fit.remove entry. ; Bug fix in calculation of psfwidth (sometimes gave NaN). ; FEB-2008, Paul Hick (UCSD/CASS), V5.11 ; Added USNO asteroids and 'stars_not_subtracted' catalogue to ; list of point sources (these are fitted, but not subtracted). ; MAR-2008, Paul Hick (UCSD/CASS), V5.12 ; Added keyword SMEI_SKY, SMEI_CAL and SMEI_ORB to pnt header ; MAR-2008, Paul Hick (UCSD/CASS), V5.13 ; Fixed problem with two DATE keywords in header of output ; Fits file. ; MAR-2008, Paul Hick (UCSD/CASS), V5.2 ; Bug fix in edge-clean procedure. The bug may have been ; introduced with V5.01 (JUL-2007) with the inclusion of the ; 'flat' star-subtracted maps. ; APR-2008, Paul Hick (UCSD/CASS), V5.3 ; Both sidereal background and zodiacal dust brightness ; are now subtracted prior to fitting stars. ; The weighting in href=smei_star_lsqfit= was changed to ; from 1/b to 1/(50+b), where b is the residual brightness ; (after background and zodiac subtraction). ; Added keyword removezld. ; Added header keywords ZLD, ZLDVER to output files. ; Renamed header keyword gain from BKGAIN to SKYGAIN ; APR-2008, Paul Hick (UCSD/CASS), V5.31 ; Replaced sidereal background map with new version. ; JUN-2008, Paul Hick (UCSD/CASS) ; Added PSF angle to output Fits file. ; JUL-2008, Paul Hick (UCSD/CASS), V5.32 ; Changed gain decrease for camera 1 from 1%/yr to 0.5%/yr ; (in smei_camera_gain) ; OCT-2008, Paul Hick (UCSD/CASS), V5.33 ; Fixed bug in calculation of observation time (used stime as ; time origin instead of extracted the origin from the header ; for the time map (keyword TORIGIN). ; DEC-2008, Paul Hick (UCSD/CASS), V6.00 ; New sidereal background map (version 3) ; New zodiacal dust cloud model (model 4.01) ; Gain change for camera 1 set back to -1% per year. ; AUG-2009, Paul Hick (UCSD/CASS) ; Added keyword use_bkgnd_version to be able to change the ; selected sidereal background map. ; SEP-2009, Paul Hick (UCSD/CASS) ; Added keyword use_zld_version to be able to change the ; selected zld model. ; OCT-2009, John Clover (UCSD/CASS; jclover@ucsd.edu), V6.01 ; Changed the default behavior of star fitting, the sidereal ; background is no longer removed prior to the star fit ; if background is removed, it is done so after star fit. ; Also changed default background to version 4. ; DEC-2009, Paul Hick (UCSD/CASS), V6.02 ; The star_sky array returned from smei_star_fit is now ; contains "flat" stars, i.e. the fitted stars with the ; fitted background model subtracted ; MAY-2010, John Clover (UCSD/CASS), V7.01 ; Updated to include ZLD removed FTS headers. ; 1, 2, 3 contain star removed "flat" maps (eq, np, sp) ; with ZLD removed ; 4, 5, 6 contain star removed "flat" maps (eq, np, sp) ; with ZLD not removed (unchanged from previous version) ; Removed keyword removezld (zld now always removed) ; Added keyword fit_with_nobkgnd (default off) ; JUN-2010, Paul Hick (UCSD/CASS), V7.02 ; After bugfix in construction of flat_sky in smei_star_fit ; JUL-2010, John Clover (UCSD/CASS) ; Incrimented use_bkgnd_version to 5, new background maps ; have been from engineering mode maps with the bugfix ; applied ; AUG-2012, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V7.03 ; Added two columns (for ccd_mode and area of (exteneded) ; PSF. Changed format of column Istd. Allowed value 2 ; in last column to indicate stars not fitted because there ; are not enough good bins. ;- bstar_version = 7.03 ; Version nr for bright star subtraction InitVar, to_smeidb , /key ; Write to $SMEISKY0 InitVar, destination, ([getenv('TUB'),'SMEIDB?'])[to_smeidb] InitVar, silent , 0 InitVar, keepbkgnd , /key InitVar, cleanedge , /key InitVar, get_version , /key InitVar, reverse_loop, /key InitVar, checkversion, /key InitVar, overwrite , /key InitVar, force , /key InitVar, use_bkgnd_version, 6 InitVar, fit_with_nobkgnd , /key overwrite OR= force ; Make sure the sidereal background FITS file exists ; If it exists, read the equatorial, north and south polar background maps tmp = filepath(root=smei_filepath(mode='sid'),subdir='v'+strcompress(use_bkgnd_version,/rem),'c0sid*.fts.gz') bkgnd_map = (file_search(tmp))[0] have_bkgnd = bkgnd_map NE '' ; Background maps exist IF get_version THEN BEGIN version = [bstar_version,-1.00] IF have_bkgnd THEN $ ; Sidereal background version version[1] = fxpar(headfits(bkgnd_map,/silent),'BKGNDVER') RETURN ENDIF from_mode = 'sky' to_mode = [ 'equ', 'pnt'] to_type = ['.fts','.txt'] ; Need two destinations: ; for skymap and time series file InitVar, destination , count=1, [destination,destination] to_remote = strpos(destination,':') NE -1 ; Set up local and remote destination ndestination = n_elements(destination) FOR i=0,ndestination-1 DO BEGIN CASE to_remote[i] OF 0: BEGIN boost, local_destination , destination[i] boost, remote_destination, '' END 1: BEGIN boost, local_destination , getenv('TUB') boost, remote_destination, destination[i] END ENDCASE ENDFOR given_map = smei_getfile(wanted_map, $ mode = from_mode , $ count = count , $ _extra = _extra ) IF count EQ 0 THEN $ RETURN CASE have_bkgnd OF 0: message, /info, 'no sidereal background, '+hide_env(tmp) 1: BEGIN message, /info, 'sidereal background, '+hide_env(bkgnd_map) eq_bkgnd0 = readfits(bkgnd_map, hdr, exten_no=0, /silent); Equatorial map bkgnd_version = fxpar(hdr, 'BKGNDVER') ; Version number of sidereal background map eq_cra = fxpar(hdr,'CRPIX1')-1 eq_cdec = fxpar(hdr,'CRPIX2')-1 eq_bdeg = 1.0/fxpar(hdr,'CDELT2') np_bkgnd0 = readfits(bkgnd_map, hdr, exten_no=1, /silent); North pole map np_cra = fxpar(hdr,'CRPIX1')-1 np_cdec = fxpar(hdr,'CRPIX2')-1 np_bdeg = 1.0/fxpar(hdr,'CDELT2') sp_bkgnd0 = readfits(bkgnd_map, hdr, exten_no=2, /silent); South pole map sp_cra = fxpar(hdr,'CRPIX1')-1 sp_cdec = fxpar(hdr,'CRPIX2')-1 sp_bdeg = 1.0/fxpar(hdr,'CDELT2') END ENDCASE removebkgnd = have_bkgnd AND (1-keepbkgnd) IF NOT removebkgnd THEN $ IF silent LE 2 THEN $ message, /info,'not subtracting sidereal background' stop_time = TimeSystem(/silent) IF silent LE 3 THEN message, /info, strcompress(count,/rem)+' '+from_mode+' map(s)' imap_low = 0L *(1-reverse_loop)+(count-1)*reverse_loop imap_high = (count-1)*(1-reverse_loop)+ 0L *reverse_loop imap_step = 1L *(1-reverse_loop)+( -1L)*reverse_loop FOR imap=imap_low,imap_high,imap_step DO BEGIN start_time = stop_time skyfile = given_map[imap] skyhide = hide_env(skyfile) skytime = smei_filename(skyfile,camera=camera,postfix=postfix) destroyvar, workname, workfile bad_dir = 0 FOR i=0,ndestination-1 DO BEGIN boost, workname, smei_filename(skytime,camera=camera,mode=to_mode[i],type=to_type[i]) CASE 1 OF to_remote[i]: tmp = local_destination[i] ; Local destination is $TUB to_smeidb : tmp = smei_filepath(source=local_destination[i],mode=to_mode[i],camera=camera,postfix=postfix) ELSE : tmp = local_destination[i] ENDCASE IF NOT CheckDir(tmp) THEN bad_dir += 1 boost, workfile, filepath(root=tmp,workname[i]) ENDFOR IF bad_dir NE 0 THEN continue IF checkversion THEN BEGIN tmp = (file_search(workfile[0]+'*'))[0] IF tmp NE '' THEN BEGIN ; Skymap exists hdr = headfits(tmp) map_bstar_version = fxpar(hdr,'BSTARVER') has_bstar_version = !err EQ 0 ; CONTINUE not allowed inside CASE block IF have_bkgnd THEN BEGIN map_bkgnd_version = fxpar(hdr,'BKGNDVER') has_bkgnd_version = !err EQ 0 IF (has_bstar_version AND map_bstar_version GE bstar_version) AND $ (has_bkgnd_version AND map_bkgnd_version GE bkgnd_version) THEN BEGIN IF silent LE 2 THEN message, /info, '>= version, '+hide_env(workfile[0]) CONTINUE ENDIF ENDIF ELSE BEGIN IF has_bstar_version AND map_bstar_version GE bstar_version THEN BEGIN IF silent LE 2 THEN message, /info, '>= version, '+hide_env(workfile[0]) CONTINUE ENDIF ENDELSE ENDIF ENDIF ELSE IF NOT overwrite THEN BEGIN tmp = (file_search(workfile[0]+'*'))[0] IF tmp NE '' THEN BEGIN ; Skymap exists IF silent LE 2 THEN message, /info, hide_env(tmp)+' exists' CONTINUE ENDIF ENDIF IF silent LE 2 THEN message, /info, skyhide ; Read maps and headers; then remove stars smei_sky,skyfile,/exists,/equatraw,sky=eq_raw,hdr=eq_struct,/noplot,silent=silent+1,nbin=1, exten_no=eq_exten eq_exists = size(eq_raw,/n_dim) NE 0 eq_full = n_elements(eq_raw) EQ 3600L*1200L IF eq_exists THEN eq_good = where(finite(eq_raw)) ELSE eq_good = -1 ccd_mode = eq_struct.ccd_mode smei_sky,skyfile,/exists,/northraw,sky=np_raw,hdr=np_struct,/noplot,silent=silent+1,nbin=1, exten_no=np_exten np_exists = size(np_raw,/n_dim) NE 0 np_full = n_elements(np_raw) EQ 800L*800L IF np_exists THEN np_good = where(finite(np_raw)) ELSE np_good = -1 smei_sky,skyfile,/exists,/southraw,sky=sp_raw,hdr=sp_struct,/noplot,silent=silent+1,nbin=1, exten_no=sp_exten sp_exists = size(sp_raw,/n_dim) NE 0 sp_full = n_elements(sp_raw) EQ 800L*800L IF sp_exists THEN sp_good = where(finite(sp_raw)) ELSE sp_good = -1 skygain = smei_camera_gain(skytime, camera=camera) ;view, in=rebin(eq_raw < 500,900,300),/stretch ; Calculate the sidereal background and the zodiacal light ; distribution for all good bins IF NOT (eq_full AND np_full AND sp_full) THEN message, 'not yet implemented' ; Set the weights to be used in the linear fit in smei_star_lsqfit ; from the skymaps prior to subtracting sidereal background and ; zodiacal light distribution. This makes the weights (and hence the ; result of the fit) largely independent of the background. ; The weights will only be used (in smei_star_lsqfit) if /use_weights is set. ; (For now we are depending on the default in smei_star_lsqfit of 1/(50+brightness) ;eq_weights = 1.0/eq_raw ;np_weights = 1.0/np_raw ;sp_weights = 1.0/sp_raw IF eq_good[0] NE -1 THEN BEGIN IF have_bkgnd THEN eq_bkgnd = eq_bkgnd0[eq_good]*skygain ELSE eq_bkgnd = 0 tmp = ArrayLocation(eq_good,dim=size(eq_raw,/dim)) tmp = cvsmei(from_map=tmp,mode='eqmap',/to_equatorial,/degrees,/silent) eq_zld = skygain*smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent, $ version=zld_version,use_zld_version=use_zld_version, camera=camera) ENDIF IF np_good[0] NE -1 THEN BEGIN IF have_bkgnd THEN np_bkgnd = np_bkgnd0[np_good]*skygain ELSE np_bkgnd = 0 tmp = ArrayLocation(np_good,dim=size(np_raw,/dim)) tmp = cvsmei(from_map=tmp,mode='npmap',/to_equatorial,/degrees,/silent) np_zld = skygain*smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent, $ version=zld_version,use_zld_version=use_zld_version, camera=camera) ENDIF IF sp_good[0] NE -1 THEN BEGIN IF have_bkgnd THEN sp_bkgnd = sp_bkgnd0[sp_good]*skygain ELSE sp_bkgnd = 0 tmp = ArrayLocation(sp_good,dim=size(sp_raw,/dim)) tmp = cvsmei(from_map=tmp,mode='spmap',/to_equatorial,/degrees,/silent) sp_zld = skygain*smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent, $ version=zld_version,use_zld_version=use_zld_version, camera=camera) ENDIF ; Remove zodiacal light prior to star subtraction ; If fit_with_nobkgnd is set then The sidereal background is removed ; before subtracting stars. Presumably this makes the background near the galactic plane ; flatter (removing first and higher spatial derivatives), and thus ; improving the star fitting. IF silent LE 2 THEN $ message, /info, 'remove zodiacal light prior to star subtraction' IF eq_good[0] NE -1 THEN eq_raw[eq_good] -= eq_zld+fit_with_nobkgnd*eq_bkgnd ; Equator IF np_good[0] NE -1 THEN np_raw[np_good] -= np_zld+fit_with_nobkgnd*np_bkgnd ; North pole IF sp_good[0] NE -1 THEN sp_raw[sp_good] -= sp_zld+fit_with_nobkgnd*sp_bkgnd ; South pole ;view, in=rebin(eq_raw < 500,900,300),/new,/stretch ; Fit and subtract the stars: *_raw -> *_sky smei_star_fit, hdr=eq_struct, sky=eq_raw, final_sky=eq_sky , $ flat_sky=eq_flat, star_sky=eq_stars, star_fit=eq_fit , $ count=eq_count, /degrees, badfit=0.0, silent=silent , $ _extra=_extra, star_info=star_info, weights=eq_weights , $ ccd_mode=ccd_mode print print smei_star_fit, hdr=np_struct, sky=np_raw, final_sky=np_sky , $ flat_sky=np_flat, star_sky=np_stars, star_fit=np_fit , $ count=np_count, /degrees, badfit=0.0, silent=silent , $ _extra=_extra, /northpole, weights=np_weights , $ ccd_mode=ccd_mode print print smei_star_fit, hdr=sp_struct, sky=sp_raw, final_sky=sp_sky , $ flat_sky=sp_flat, star_sky=sp_stars, star_fit=sp_fit , $ count=sp_count, /degrees, badfit=0.0, silent=silent , $ _extra=_extra, /southpole, weights=sp_weights , $ ccd_mode=ccd_mode print print IF silent LE 2 THEN BEGIN tmp = [eq_count,np_count,sp_count] message, /info, strjoin(strcompress(tmp,/rem),'+')+'=' + $ strcompress(round(total(tmp)),/rem)+' stars, '+skyhide ENDIF IF cleanedge AND NOT (eq_full AND np_full AND sp_full) THEN BEGIN cleanedge = 0 message, /info, '/cleanedge ignored: not a full-skymap, '+skyhide ENDIF ;view, in=rebin(eq_sky > (-25) < 25,900,300),/new,/stretch ;view, in=rebin(eq_flat > (-25) < 25,900,300),/new,/stretch ;return ; At this point the zodiacal light and all the stars have been removed. ; If fit_with_nobkgnd is set then the sidereal background is also gone. ; We now have two forms of star-subtracted maps: ; eq_sky, np_sky, sp_sky: in these maps the fitted PSF has been subtracted. ; eq_flat, np_flat, sp_flat: in these maps the planar background has ; been substituded for the star. ; In previous versions (< 7.0) both these maps where stored in the ; the output file: *_sky maps in Fits extensions 0,1,2; *_flat in ; extensions 3,4,5. ; Since version 7.0 the *_sky maps are not longer output. Instead ; ony the *_flat maps are used. ; eq_sky (1st three maps): sidereal, stars and zld removed ; eq_flat (2nd three maps): sidereal, stars removed; zld still in ; The sidereal background removal is optional (depending on removebkgnd setting) ; So eq_sky is eq_flat with the zld removed CASE cleanedge OF ; Edge cleaning skipped 0: BEGIN ; Subtract the sidereal map from the "flat" skymaps ; if removebkgnd is set. ; Copy the "flat" map into the *_sky maps. ; Then add the zld model back into the "flat" maps IF eq_good[0] NE -1 THEN BEGIN eq_flat[eq_good] += (fit_with_nobkgnd-removebkgnd)*eq_bkgnd eq_sky = eq_flat eq_flat[eq_good] += eq_zld ; Equator ENDIF IF np_good[0] NE -1 THEN BEGIN np_flat[np_good] += (fit_with_nobkgnd-removebkgnd)*np_bkgnd np_sky = np_flat np_flat[np_good] += np_zld ; North pole ENDIF IF sp_good[0] NE -1 THEN BEGIN sp_flat[sp_good] += (fit_with_nobkgnd-removebkgnd)*sp_bkgnd sp_sky = sp_flat sp_flat[sp_good] += sp_zld ; South pole ENDIF END 1: BEGIN IF silent LE 2 THEN $ message, /info, 'scattered light at edge of FOV' ; Add the zodiacal light and the sidereal background (if ; fit_with_nobkgnd is set) back in ; (these contribute to the scattered light removed by smei_sky_cleanedge). IF eq_good[0] NE -1 THEN BEGIN tmp = eq_zld+fit_with_nobkgnd*eq_bkgnd ; should this have eq_zld+tmp or just tmp like the others? eq_sky [eq_good] += tmp ; Equator eq_flat[eq_good] += tmp ; Equator ENDIF IF np_good[0] NE -1 THEN BEGIN tmp = np_zld+fit_with_nobkgnd*np_bkgnd np_sky [np_good] += tmp ; North pole np_flat[np_good] += tmp ; North pole ENDIF IF sp_good[0] NE -1 THEN BEGIN tmp = sp_zld+fit_with_nobkgnd*sp_bkgnd sp_sky [sp_good] += tmp ; South pole sp_flat[sp_good] += tmp ; South pole ENDIF ; Get the low res data needed for the scattered light removal smei_sky, skyfile, /exists, /fovx , sky=fovx_map, /noplot, silent=silent+1, /degrees smei_sky, skyfile, /exists, /orbsecs, sky=time_map, /noplot, silent=silent+1 ; Remove scattered light from just outside fov ; Concatenate eq_sky and eq_flat. ; eq_sky is used to calculate the scattered light ; The same is subtracted from both eq_sky and eq_flat ; (same for np and sp maps) eq_sky = reform([[eq_sky],[eq_flat]], [size(eq_sky,/dim),2]) np_sky = reform([[np_sky],[np_flat]], [size(np_sky,/dim),2]) sp_sky = reform([[sp_sky],[sp_flat]], [size(sp_sky,/dim),2]) smei_sky_cleanedge_fov , $ eq_sky = eq_sky , $ np_sky = np_sky , $ sp_sky = sp_sky , $ time_map= time_map , $ fovx_map= fovx_map , $ /degrees , $ silent = silent eq_flat = eq_sky[*,*,1] eq_sky = eq_sky[*,*,0] np_flat = np_sky[*,*,1] np_sky = np_sky[*,*,0] sp_flat = sp_sky[*,*,1] sp_sky = sp_sky[*,*,0] ; If necessary, remove sidereal background and/or ; zodiacal dust brightness. ; If one of the *_good arrays is -1, then the corresponding ; *_flat and *_sky arrays should be entirely filled with NaNs. IF eq_good[0] NE -1 THEN BEGIN eq_flat[eq_good] -= removebkgnd*eq_bkgnd ; Equator eq_sky = eq_flat eq_sky [eq_good] -= eq_zld ; Equator ENDIF IF np_good[0] NE -1 THEN BEGIN np_flat[np_good] -= removebkgnd*np_bkgnd ; North pole np_sky = np_flat np_sky [np_good] -= np_zld ; North pole ENDIF IF sp_good[0] NE -1 THEN BEGIN sp_flat[sp_good] -= removebkgnd*sp_bkgnd ; South pole sp_sky = sp_flat sp_sky [sp_good] -= sp_zld ; South pole ENDIF END ENDCASE ;view, in=rebin(eq_sky < 500,900,300),/new,/stretch ;return ; =================================================== ; Start writing the Fits file for the subtracted maps. ; hdr is the primary header from the original skymap ; Mark bad data with BAD_DATA value from Fits header hdr = headfits(skyfile, silent=silent, exten=eq_exten) bad_data = fxpar(hdr,'BAD_DATA') sky_version = fxpar(hdr,'SMEI_SKY') IF !err LT 0 THEN sky_version = 0.00 cal_version = fxpar(hdr,'SMEI_CAL') IF !err LT 0 THEN cal_version = 0.00 orb_version = fxpar(hdr,'SMEI_ORB') IF !err LT 0 THEN orb_version = 0.00 IF eq_exists THEN BEGIN bad = where(1-finite(eq_sky)) IF bad[0] NE -1 THEN eq_sky [bad] = bad_data bad = where(1-finite(eq_flat)) IF bad[0] NE -1 THEN eq_flat [bad] = bad_data bad = where(1-finite(eq_stars)) IF bad[0] NE -1 THEN eq_stars[bad] = bad_data ENDIF IF np_exists THEN BEGIN bad = where(1-finite(np_sky)) IF bad[0] NE -1 THEN np_sky [bad] = bad_data bad = where(1-finite(np_flat)) IF bad[0] NE -1 THEN np_flat [bad] = bad_data bad = where(1-finite(np_stars)) IF bad[0] NE -1 THEN np_stars[bad] = bad_data ENDIF IF sp_exists THEN BEGIN bad = where(1-finite(sp_sky)) IF bad[0] NE -1 THEN sp_sky [bad] = bad_data bad = where(1-finite(sp_flat)) IF bad[0] NE -1 THEN sp_flat [bad] = bad_data bad = where(1-finite(sp_stars)) IF bad[0] NE -1 THEN sp_stars[bad] = bad_data ENDIF ; Set up Fits header template ; Make main header. Copy most entries from the main header ; of the original skymap. Fits keyword SMEI_SKY (was SMEI_HTM) ; comes right after the mandatory keywords. mkhdr, main_hdr, IsType(eq_sky), size(eq_sky,/dim), /extend tmp = (where(strpos(hdr,'SMEI_SKY') EQ 0))[0] IF tmp EQ -1 THEN tmp = (where(strpos(hdr,'SMEI_HTM') EQ 0))[0] ; This cuts of the IDL Fits header after the EXTEND keyword, ; dropping the DATE keyword and a long COMMENT line. ; (there is another DATE keyword copied from the primary header ; of the input skymap, which is updated below). main_hdr = [main_hdr[0:(where(strpos(main_hdr,'DATE =') EQ 0))[0]-1],hdr[tmp:*]] ; Main header updates fxaddpar, main_hdr, 'NAME' , GetFileSpec(workfile[0],part='name') fxaddpar, main_hdr, 'IMG_TYPE', 'skymap_'+to_mode[0] fxaddpar, main_hdr, 'DATE' , TimeGet(TimeSystem(/silent),format='YYYY-MN-DD hh:mm',/scalar) ; Main header additions. Add new keywords in reverse order as needed in Fits file. fxaddpar, main_hdr, 'SKYGAIN' , skygain , ' Gain factor', after='BAD_DATA' fxaddpar, main_hdr, 'ZLDVER' , zld_version , ' Zodiacal light model' , after='BAD_DATA' ; Don't write byte into Fits file! IF removebkgnd THEN $ fxaddpar, main_hdr, 'BKGNDVER', bkgnd_version, ' Sidereal background version' , after='BAD_DATA', format='(F5.2)' ; Don't write byte into Fits file! fxaddpar, main_hdr, 'BKGND' , fix(removebkgnd) , ' Sidereal background removed (0=N,1=Y)' , after='BAD_DATA' fxaddpar, main_hdr, 'BSTARVER' , bstar_version , ' Star subtraction version' , after='BAD_DATA', format='(F5.2)' fxaddpar, main_hdr, 'BSTAR' , 1 , ' Bright stars subtracted (0=N,1=Y)', after='BAD_DATA' ; Basic header for equatorial maps mkhdr, eq_hdr, IsType(eq_sky), size(eq_sky,/dim), /image eq_hdr = [eq_hdr[0:(where(strpos(eq_hdr,'END ') EQ 0))[0]-1],main_hdr[(where(strpos(main_hdr,'MAP') EQ 0))[0]:(where(strpos(main_hdr,'BAD_DATA') EQ 0))[0]]] ; Basic header for maps of north pole hdr = headfits(skyfile, silent=silent, exten=np_exten) mkhdr, np_hdr, IsType(np_sky), size(np_sky,/dim), /image np_hdr = [np_hdr[0:(where(strpos(np_hdr,'END ') EQ 0))[0]-1],hdr[(where(strpos(hdr,'MAP') EQ 0))[0]:*]] ; Basic header for maps of south pole hdr = headfits(skyfile, silent=silent, exten=sp_exten) mkhdr, sp_hdr, IsType(sp_sky), size(sp_sky,/dim), /image sp_hdr = [sp_hdr[0:(where(strpos(sp_hdr,'END ') EQ 0))[0]-1],hdr[(where(strpos(hdr,'MAP') EQ 0))[0]:*]] ; =========== ; Start writing file ; Extension 0, equatorial map with stars and zld removed hdr = main_hdr fxaddpar, hdr, 'MAP', 'Equatorial map RA=[0,360],DEC=[-60,+60] (no zld)', ' Map description' fxaddpar, hdr, 'CLNEDGE', (['F','T'])[cleanedge], ' Clean edge' writefits, workfile[0], eq_sky, hdr ; Extension 1, map of north pole with stars and zld removed hdr = np_hdr fxaddpar, hdr, 'MAP', 'Map of north pole DEC=[+50,+90] (no zld)', ' Map description' writefits, workfile[0], np_sky, hdr, /append ; Extension 2, map of south pole with stars and zld removed hdr = sp_hdr fxaddpar, hdr, 'MAP', 'Map of south pole DEC=[-50,-90] (no zld)', ' Map description' writefits, workfile[0], sp_sky, hdr, /append ; Extension 3, equatorial map with stars removed, but zld still in hdr = eq_hdr fxaddpar, hdr, 'MAP', 'Equatorial map RA=[0,360],DEC=[-60,+60] (with zld)', ' Map description' writefits, workfile[0], eq_flat, hdr, /append ; Extension 4, map of north pole with stars removed, but zld still in hdr = np_hdr fxaddpar, hdr, 'MAP', 'Map of north pole DEC=[+50,+90] (with zld)', ' Map description' writefits, workfile[0], np_flat, hdr, /append ; Extension 5, map of south pole with stars removed, but zld still in hdr = sp_hdr fxaddpar, hdr, 'MAP', 'Map of south pole DEC=[-50,-90] (with zld)', ' Map description' writefits, workfile[0], sp_flat, hdr, /append ; Extension 6, equatorial map of removed stars hdr = eq_hdr fxaddpar, hdr, 'MAP' , 'Stars removed from equator RA=[0,360],DEC=[-60,+60]', ' Map description' fxaddpar, hdr, 'NSTARS', round(eq_count), ' Number of stars removed' writefits, workfile[0], eq_stars, hdr, /append ; Extension 7, map of stars removed from north pole hdr = np_hdr fxaddpar, hdr, 'MAP' , 'Stars removed from north pole DEC=[+50,+90]', ' Map description' fxaddpar, hdr, 'NSTARS', round(np_count), ' Number of stars removed' writefits, workfile[0], np_stars, hdr, /append ; Extension 8, map of stars removed from south pole hdr = sp_hdr fxaddpar, hdr, 'MAP' , 'Stars removed from south pole DEC=[-50,-90]', ' Map description' fxaddpar, hdr, 'NSTARS', round(sp_count), ' Number of stars removed' writefits, workfile[0], sp_stars, hdr, /append ; Copy two time extensions from original skymap ; Extension 9, time (fraction of orbit) RA=[0,360],DEC=[-90,+90] tmp = readfits(skyfile, hdr, exten_no=5+2, /silent) writefits, workfile[0], tmp, hdr, /append ; Extension 10, time (sec since TORIGIN) RA=[0,360],DEC=[-90,+90] tmp = readfits(skyfile, hdr, exten_no=5+9, /silent) writefits, workfile[0], tmp, hdr, /append ; Extension 11, Direction cosine angle-x RA=[0,360],DEC=[-90,+90] tmp = readfits(skyfile, hdr, exten_no=5+10, /silent) writefits, workfile[0], tmp, hdr, /append ; Extension 12, contributing pixel count RA=[0,360],DEC=[-90,+90] tmp = readfits(skyfile, hdr, exten_no=5+11, /silent) writefits, workfile[0], tmp, hdr, /append ; Extension 13, PSF position angle from north RA=[0,360],DEC=[-90,+90] tmp = readfits(skyfile, hdr, exten_no=5+3, /silent) writefits, workfile[0], tmp, hdr, /append ; File complete ; ============= ; /compress on writefits is ignored with /append present, so gzip here. tmp = gzip_file(workfile[0], /force) workfile[0] += '.gz' IF silent LE 3 THEN message, /info, hide_env(workfile[0])+' is '+to_mode[0]+' map' dec_cutoff = 55.0 ; Same units as RA,dec in fit_* structure (degrees) ; Write the .pnt text file. ; First write the header only IF IsType(eq_fit,/defined) OR IsType(np_fit,/defined) OR IsType(sp_fit,/defined) THEN BEGIN IF silent LE 3 THEN $ message, /info, hide_env(workfile[1])+' is '+to_mode[1]+' file' boost, star_info, $ ['SMEI_SKY: '+string(sky_version ,format='(F5.2)') , $ 'SMEI_CAL: '+string(cal_version ,format='(F5.2)') , $ 'SMEI_ORB: '+string(orb_version ,format='(F5.2)') , $ 'BSTAR : yes' , $ 'BSTARVER: '+string(bstar_version,format='(F5.2)') ] boost, star_info, 'BKGND : '+(['no','yes'])[removebkgnd] IF removebkgnd THEN $ boost, star_info, 'BKGNDVER: '+string(bkgnd_version,format='(F5.2)') boost, star_info, 'ZLDVER : '+string(zld_version,format='(F5.2)') boost, star_info, 'SKYGAIN : '+strcompress(skygain,/remove_all) boost, star_info, 'CLNEDGE : '+(['no','yes'])[cleanedge] oneorbit = fxpar(main_hdr, 'ONEORBIT') boost, star_info, 'ONEORBIT: '+(['no','yes'])[oneorbit] IF oneorbit THEN BEGIN orbitnum = fxpar(main_hdr,'ORBIT') boost, star_info, 'ORBIT : '+strcompress(orbitnum, /rem) ENDIF smei_star_writepnt, workfile[1] , $ dec_cutoff = dec_cutoff , $ camera = camera , $ star_info = star_info , $ /degrees ENDIF ELSE $ IF silent LE 3 THEN message, /info, 'no time series file created' ; Write the structures (nothing happens if the structure doesn't exist) smei_star_writepnt, workfile[1], eq_fit, dec_cutoff=dec_cutoff, /append, /degrees smei_star_writepnt, workfile[1], np_fit, dec_cutoff=dec_cutoff, /append, /degrees smei_star_writepnt, workfile[1], sp_fit, dec_cutoff=dec_cutoff, /append, /degrees FOR i=0,ndestination-1 DO BEGIN ; If destination is a remote directory then scp the files ; and delete the local copies IF to_remote[i] THEN BEGIN tmp = strpos(remote_destination[i],':')+1 CASE to_smeidb OF 0: tmp = remote_destination[i] 1: tmp = strmid(remote_destination[i],0,tmp) + $ smei_filepath(source=strmid(remote_destination[i],tmp),mode=to_mode[i],camera=camera,postfix=postfix) ENDCASE tmp = 'scp '+workfile[i]+' '+filepath(root=tmp,workname[i]) print, tmp spawn, tmp tmp = do_file(workfile[i], /delete) ENDIF ENDFOR stop_time = TimeSystem(/silent) IF silent LE 2 THEN BEGIN tmp = TimeOp(/subtract,stop_time,start_time,TimeUnit(/min)) tmp = string(tmp,format='(F10.2)') message, /info, skyhide+' '+strtrim(tmp,2)+' min' print print ENDIF ENDFOR RETURN & END