;+ ; NAME: ; smei_star_fit ; PURPOSE: ; Fit bright stars in sky map ; CATEGORY: ; ucsd/camera/idl/stars ; CALLING SEQUENCE: PRO smei_star_fit, file , $ hdr = hdr , $ sky = sky , $ camera = camera , $ ccd_mode = ccd_mode , $ origin = origin , $ scale = scale , $ southpole = southpole , $ northpole = northpole , $ torigin = torigin , $ time_map = time_map , $ psf_map = psf_map , $ fovx_map = fovx_map , $ nofovx = nofovx , $ stime = stime , $ etime = etime , $ final_sky = final_sky , $ flat_sky = flat_sky , $ star_sky = star_sky , $ star_fit = star_fit , $ star_info = star_info , $ count = count , $ upto_mags = upto_mags , $ bkgnd_origin= bkgnd_origin , $ bkgnd_radius= bkgnd_radius , $ bkgnd_count = bkgnd_count , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ auto_wing = auto_wing , $ star_list = star_list , $ star_remove = star_remove , $ star_index = star_index , $ catalogues = catalogues , $ degrees = degrees , $ kill_one = kill_one , $ fix_centroid= fix_centroid , $ fix_fovangle= fix_fovangle , $ fix_psfangle= fix_psfangle , $ fix_pattern = fix_pattern , $ noplane = noplane , $ noyzero = noyzero , $ use_filled_psf=use_filled_psf,$ use_weights = use_weights , $ weights = weights , $ skip_dist = skip_dist , $ incl_dist = incl_dist , $ max_stars = max_stars , $ badfit = badfit , $ show = show , $ area_show = area_show , $ max_show = max_show , $ min_show = min_show , $ silent = silent , $ sigma = sigma , $ magnify = magnify , $ ;no_stars_mask = no_stars_mask , $ no_corepsf_mask = no_corepsf_mask, $ rmbkgnd = rmbkgnd , $ rmzld = rmzld , $ known_stars = known_stars , $ known_names = known_names , $ known_mags = known_mags ; INPUTS: ; file scalar; type; string ; fully-qualified file name of sky map ; OPTIONAL INPUT PARAMETERS: ; /degrees if set then all angles are in degrees (default: radians) ; ; hdr=hdr array[1]; type: structure ; sky map header ; sky=sky array[n,m]; type: float ; 2D sky map ; ; camera=camera ; ccd_mode=ccd_mode ; origin=origin ; scale=scale ; ; /southpole if set, then fit stars in north polar map ; /northpole if set, then fit stars in south polar map ; By default stars are fit in the equatorial map ; ; star_list=star_list array; type: SMEI_STAR_LIST structure ; list of stars to be fitted ; These can be extracted from the SMEI ; star catalogues with href=smei_star_list=. ; array: type: string ; alternatively a list of star names are specified ; A call to href=smei_star_list= is made to convert ; the list to a SMEI_STAR_LIST structure ; catalogues=catalogues scalar or array; type: string; ; list of star catalogues to be used ; Passed to href=smei_star_list=. ; NOT used if keyword star_list is set. ; upto_mags=upto_mags scalar; type: float; default: 6.0 ; maximum SMEI magnitude; only stars ; brighter than this are removed. ; All other objects are fitted, but NOT removed ; from the skymaps ; bkgnd_origin=bkgnd_origin ; array[2]; type: float ; eventually processed by href=smei_star_standard=. ; bkgnd_radius=bkgnd_radius ; scalar or array[2]; type: float ; eventually processed by href=smei_star_standard=. ; bkgnd_count=bkgnd_count ; scalar; type: float (<= 1.0) or integer (> 1) ; passed to href=smei_star_lsqbins=. ; ; wing_origin=wing_origin array[2]; type: float ; eventually processed by href=smei_star_standard=. ; wing_radius=wing_radius array[2]; type: float ; eventually processed by href=smei_star_standard=. ; /auto_wing if set, then wing_radius varies with the SMEI magnitude as ; wing_radius+0.1*(3-m_SMEI) degrees ; i.e. the wing radius is smaller for fainter stars by ; one high-res skybin (0.1 degrees) per magnitude. ; wing_radius is set to 1.4 degrees if it is not defined. ; (note that href=smei_star_standard= will prevent the ; wing_radius to drop below a minimum value though). ; This keyword was added for use in the SMEI pipeline only, ; to define a larger wing radius for very bright stars. ; /fix_centroid if set, then the star centroid is included in the fit. ; If multiple stars are fitted simultaneously then this will ; adjust the centroid of the group of stars while keeping the ; relative locations of the stars fixed (as implied by the ; catalogue positions). ; This keyword is primarily useful for correcting for pointing ; errors, which affect all stars in a group in the same way. ; /fix_pattern /fix_pattern implies /fix_centroid. ; /fix_pattern has no effect for stars fitted individually. ; if set, then for stars fitted simultaneously the centroids ; of all stars in the group are adjusting individually ; (changing the relative locations of the stars) to improve the fit. ; /fix_psfangle if set, then the PSF angle is included in the fit ; /fix_fovangle if set, then the direction cosine angles are included ; in the fit. ; (this slows down the star fitting considerably). ; ; /noplane if set, suppress fit of planar background, and fit a ; constant background instead (default: fit ; planar background). ; /noyzero passed to href=smei_star_lsqfit= ; /use_filled_psf passed to href=smei_star_standard= ; /use_weights passed to href=smei_star_lsqfit= ; ; weights=weights array[n,m]; type: float ; Only used if /use_weights SET. ; Weights to be used for skybins ; ; /no_corepsf_mask by default, when a star is fitted, all skybins in the core PSF ; of stars brighter then m_smei+2 are excluded from being ; using in the fitting procedure. ; Set this keyword to bypass this (bad idea!) ; ; /kill_one by default, when multiple (nearby) stars are fitted simultaneously ; all stars in the group are subtracted. If /kill_one is set the ; group is fitted as a whole, but only the brightest star is subtracted ; and the other stars will be fitted separately later. ; ; skip_dist=skip_dist scalar; type: float; default: 0.25 degrees ; nearby fainter stars at angular distances less ; then 'skip_dist' away are not fit ; incl_dist=incl_dist scalar; type: float; default: 0.75 degrees ; nearby fainter stars at angular distances ; less than 0.75 degrees (and more than /skip_dist) ; are included when fitting a star. ; max_stars=max_stars scalar; type: integer; default: 4 ; maximum number of stars fitted simultaneously. ; If more than 'max_stars' stars pass the ; 'incl_dist' criterion, than only the ; brightest 'max_stars-1' are included in the fit ; (for a fit of max_stars stars). ; sigma=sigma scalar; type: float ; ? ; show=show scalar; type: star ; name of star for which results are displayed ; This should be one of the names from the ; SMEI star catalogue. ; Set show='*' to display results for all stars ; Five views of the star are displayed: ; - the area around the star in the original skymap ; - the same area after the fit is subtracted ; - the same area with the PSF replaced by the fitted background ; - the same area showing the bins used in the fit ; - a scatter plot of model vs. data ; ; min_show scalar or array[2]; type: numeric; default: none ; max_show scalar or array[2]; type: numeric; default: none ; when keyword 'show' is set the display is auto-scaled ; by default. These two keywords can be used to ; set the dynamic range. If a scalar is set then the ; value is used for all three diplays of the area around the ; star. If a 2-element array is specified then the ; 1st value is used for the first display (the original ; skymap), and the 2nd is used for the two subtracted ; displays. ; ; Sometimes it is convenient to be able to specify characteristics of certain 'known' stars, ; rather than fitting them. There are two mechanisms in place to do this. ; ; known_stars = known_stars ; array[3,*]; type: float ; specifies RA, dec and SMEI magnitude of a list of 'known' stars. ; These stars are removed from the map prior to fitting other ; stars by subtracting the standard star multiplied by the SMEI ; brightness (as determined by smei_i2m from the SMEI magnitude) ; at the specified location in the skymap. ; Note that it is the responsibility of the user to find the ; correct location and SMEI magnitude of the star. Note also ; that these star do not need to be in one of the SMEI ; star catalogues. ; ; known_names array; type: string ; known_mags array, same size as 'known_names'; type: numeric ; list of stars (catalogue names and SMEI magnitudes) ; of stars for which the brightness is considered ; known. The brightness for these stars is not fitted. ; Typically used when a non-variable star (with known ; brightness), intrudes on the location of a variable ; star of interest. Both stars can be fit simultaneously, ; while keeping the brightness of the non-variable star ; fixed, i.e. only the background and brightness of the ; variable star is fit. ; OPTIONAL OUTPUT PARAMETERS: ; count=count scalar; type: integer ; number of stars fitted (but not necessarily ; subtracted) ; count=0 if the skymap was not found or if ; no stars with valid intensity or psf angle ; were found. ; If count=0 then star_fit will not exist. ; star_fit=star_fit array[count]; SMEI_STAR_FIT structure ; contains the fit results for all stars ; star_index=star_index array[count]; integer ; position of fitted star in original ; star_list (useful only if star_list is ; is one of the input arguments) ; If count=0 then star_index = -1 ; final_sky=final_sky array[3600,1200] or array[800,800] ; skymap with stars subtracted ; flat_sky=flat_sky array[3600,1200] or array[800,800] ; skymap with stars subtracted, but also with ; the core PSF for each fitted star explicitly ; replaced by the background model ; star_sky=star_sky array[3600,1200] or array[800,800] ; skymap containing the subtracted stars ; ; If there is a problem reading the skymap Fits files ; these variables will contain a scalar -1 ; If count=0 then final_sky and sky will contain the ; unmodified input skymap, and star_sky will contain ; zeroes. ; ; area_show=area_show array[50,50,4]; type: float ; If keyword 'show' is used, this returns the content ; of the the four displays showing the area around ; the star (after applying the min_show and max_show ; keywords). Note that only the last display is saved, ; so this keyword only makes sense if a single star ; is specified in the 'show' keyword. ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, IsType, ToRadians, ToDegrees, who_am_i, smei_star_list ; smei_star_info, smei_sky, AngleRange, AngleUnits, destroyvar ; MagnifyArray, smei_star_stdmaps, smei_star_lsqfit, gridgen, TimeSet ; TimeUnit, hide_env, IsTime, TimeOp, smei_coriolis, jpl_body, big_eph ; TimeLimits, BadValue, TimeGet, sphere_distance, view, twin ; LocalExtrema, smei_star_standard, smei_star_formatpnt ; smei_star_lsqbins, smei_star_enough_bins ; EXTERNAL: ; smei_star__define ; PROCEDURE: ; The information needed for the star removal is: ; ; sky 2D-sky map ; map type type of sky map (south pole, north pole, equatorial map) ; origin array indices for RA=dec=0 (eq map) or indices of pole (polar map) ; resolution number of degrees/radians per bin ; camera SMEI camera (each camera has its own standard star map) ; ; psf_map equatorial 2D-sky map with rotation angles for the stars ; ; time_map equatorial 2D-sky map with times (secs since start of orbit) ; ; This information can be supplied in a number of different ways: ; ; 1. Set 'file' to the fully-qualified file name of a sky map ; All information is extracted from the file by href=smei_sky=. ; Information in other keywords is ignored. ; Note that by default the equatorial map is processed, unless ; /southpole or /northpole is set. ; ; 2. Set keyword 'hdr' to the header of a sky map previously read by ; href=smei_sky=. 'hdr' provides all the information needed, ; include the map type (equatorial, south- or northpole). ; The skymap itself can be supplied in keyword 'sky' (it probably ; is available from the same smei_sky call that set 'hdr'. If absent ; then the original skymap is read (the file name is in 'hdr'). ; All other keywords are ignored. The psf and time maps are always ; read from file by smei_sky. ; ; 3. Explicitly provide all information in the appropriate keywords. ; Defaults are set to the values currently used for the sky maps ; produced by href=smei_skyd=. ; ; ; Stars are fitted and subtracted starting with the brightest stars ; and working down to the fainter stars down to the cutoff magnitude ; 'upto_mags'. ; ; Groups of upto 'max_stars' stars are fitted simultaneously. This group ; is selected as follows: ; - the first star in the group is the brightest star in the sky ; that has not been fit yet. ; - all stars within 'incl_dist' are added to the group ; - the same is done iteratively for stars added to the group in this ; way: for each a neighbourhood of 'incl_dist' is searched, and if ; stars are found these are added to the group ; - stop if no more stars are found within 'incl_dist' of any star ; in the group, or if the group size exceeds 'max_crowd' members ; - pass through all group members in order of decreasing brightness ; and remove fainter stars that lie within 'skip_dist'. These ; stars are marked as "done" immediately. ; - only keep at most 'max_star' of the remaining stars are not ; fitted on the assumption that its psf is indistinguishable from ; the psf of the brighter star. ; ; The results of the fit are by default subtracted for all the stars ; in the group simultaneously. If /kill_single is set then only the ; fit for the bright star is subtracted. The fainter stars will be ; fitted again separately (after the bright star has been removed). ; ; By default a planar (sloped) background is fitted. If keyword ; /noplane is set then a constant background is fitted. ; ; By default the star locations in the SMEI star catalogue are used ; as fixed centroids. If keyword /fix_centroid is set then a brute ; force iteration scheme is used to improve the fit by fitting the ; centroid location also. This is mainly useful for detecting 'bad ; quaternion' episodes in the SMEI data. Setting /fix_centroids ; slows the program down considerably. ; ; Fitting of the plane background and the star brightnesses is done ; using linear least squares fits. href=lsqLinearFit= implements fits ; upto four dimension as analytic expressions, allowing for fitting ; of two stars with a planar background of three stars with a constant ; background. For fits in more dimension a matrix inversion is used. ; ; The planets Mercury,Venus,Mars,Jupiter,Saturn and Uranus are ; subtracted using their apparent visual magnitudes to control ; the order of subtraction. ; ; MODIFICATION HISTORY: ; NOV-2005, Paul Hick (UCSD/CASS) ; FEB-2006, Paul Hick (UCSD/CASS) ; Replaced keyword /plane by keyword /noplane. By default ; a sloped background is fitted unless suppressed by /noplane. ; AUG-2006, Paul Hick (UCSD/CASS) ; Keyword star_list can now also be a list of star names from ; the SMEI catalogues. ; If no stars are left to be fitted the procedure will now always ; return with count=0 and final_sky and star_sky existing. ; (it used to abort in some situations). ; SEP-2006, Paul Hick (UCSD/CASS) ; Added /use_filled_psf keyword ; Added star_info keyword ; NOV-2006, Jordan Vaughan (UCSD; jtvaugha@ucsd.edu) ; Added gain field to star_fit ; Set star_fit.close to '-' by default ; FEB-2007, Jordan Vaughan, Paul Hick (UCSD/CASS) ; Added stars_mask array, passed to smei_star_stdmaps ; to exclude skybins from previously-subtracted star ; from the background skybins. ; APR-2007, Paul Hick (UCSD/CASS) ; The lowres time map is now used as is (i.e. we do not ; take the absolute value anymore). ; MAY-2007, Paul Hick (UCSD/CASS) ; Added keyword fix_psfangle fix_fovangle ; In the iterative fitting a point far away from the ; current best fit is now used in an attempt to speed up ; the iteration when the starting point is far away from ; the correct solution (seems to be most useful if ; keyword /fix_fovangle is set). ; JUN-2007, Paul Hick (UCSD/CASS) ; Bug fix. psfwidth calculation for case with only one ; skybin insided psf would crash. ; JUL-2007, Paul Hick (UCSD/CASS) ; Substantial rewrite of section that selects the ; group of stars to be fitted simultaneously. ; JAN-2008, Paul Hick (UCSD/CASS) ; Modified storage of planet names in star_list ; structure (left-adjusted, instead of right-adjusted) ; FEB-2008, Paul Hick (UCSD/CASS) ; Added star_fit.remove entry. ; Added USNO asteroids and the 'stars_not_subtracted' ; catalogue to list of point sources. These are fitted but ; not subtracted from the sky map. ; Bug fix in calculation of psfwidth (sometimes gave NaN). ; MAR-2008, Paul Hick (UCSD/CASS) ; Bugfix. The return arrays final_sky, flat_sky and ; star_sky would not exist if count=0 (no stars subtracted). ; Now final_sky and flat_sky contain the input map, ; and star_sky contains zeroes. ; APR-2008, Paul Hick (UCSD/CASS) ; Made sure that dradec is really zero when no centroid ; fitting is done (used to be a really small number ; sometimes instead of identical to zero). ; Added keyword weights. Added keyword original_weights ; to all smei_star_lsqfit calls. ; JUL-2008, Paul Hick (UCSD/CASS) ; Added keyword /fix_pattern to allow adjusting of centroids ; of individual stars when multiple stars are fitted ; simulataneously. ; FEB-2010, Paul Hick (UCSD/CASS) ; Added keyword area_show. Added documentation. ; Fixed bug in processing of keyword star_list (would crash ; if none of the stars on the list would be in one of ; the selected SMEI star catalogues). ; JUN-2010, Paul Hick (UCSD/CASS) ; Fixed bug in the construction of flat_sky array (pieces ; of final_sky were erroneously copied into flat_sky) ; JUN-2010, John Clover (UCSD/CASS) ; Bumped incl_dist to 1.0 degrees ; APR-2012, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added more checks in the iteration procedures to pick up ; cases where there are not enough bins left for an lsq fit ; (also prevents crashes due to an 'array index out of ; range' error on an empty index array when no bins are ; left at all) ;- InitVar, show , '' IF IsType(show,/generic_int) THEN show = '*' InitVar, silent , 0 InitVar, upto_mags , 6.0 InitVar, fix_centroid , /key InitVar, fix_psfangle , /key InitVar, fix_fovangle , /key InitVar, fix_pattern , /key fix_centroid OR= fix_pattern fix_something = fix_centroid OR fix_psfangle OR fix_fovangle InitVar, nofovx , /key InitVar, kill_one , /key InitVar, degrees , /key InitVar, noplane , /key InitVar, noyzero , 0 InitVar, use_filled_psf , /key InitVar, auto_wing , /key InitVar, magnify , 0.05 ;InitVar, no_stars_mask , /key InitVar, no_corepsf_mask, /key rpm = ToRadians(degrees=degrees) dpm = ToDegrees(degrees=degrees) usec = TimeUnit(/sec) InitVar, skip_dist , 0.25/dpm InitVar, incl_dist , 1.00/dpm InitVar, max_stars , 4 InitVar, ccd_mode , -1 IF ccd_mode LT 0 OR ccd_mode GT 2 THEN ccd_mode = 9 ; Keyword auto_wing only does something if ; wing_radius is set. IF auto_wing THEN $ InitVar, wing_radius, 1.4/dpm has_sigma = IsType(sigma ,/defined) fit_bkgnd = IsType(bkgnd_radius ,/defined) fit_wing = IsType(wing_radius ,/defined) format = (['(F7.5)','(F5.3)'])[degrees] destroyvar, area_show ; Store fit parameters in return keyword star_info as string array star_info = [ $ 'm_SMEI threshold: '+strcompress(upto_mags,/rem) , $ 'fix centroid : '+(['no','yes'])[fix_centroid] , $ 'fix psfangle : '+(['no','yes'])[fix_psfangle] , $ 'fix fovangle : '+(['no','yes'])[fix_fovangle] , $ 'fix pattern : '+(['no','yes'])[fix_pattern ] , $ 'kill one : '+(['no','yes'])[kill_one ] , $ 'filled PSF : '+(['no','yes'])[use_filled_psf] , $ 'angular units : '+(['radians','degrees'])[degrees] , $ 'skip distance : '+string(skip_dist,format=format) , $ 'incl distance : '+string(incl_dist,format=format) , $ 'max crowd : '+strcompress(max_stars,/rem) , $ 'fovx correction : '+(['no','yes'])[1-nofovx] , $ 'bkgnd plane : '+(['no','yes'])[1-noplane] , $ 'bkgnd separate : '+(['no','yes'])[IsType(bkgnd_radius,/defined)], $ 'extended wing : '+(['no','yes'])[IsType(wing_radius ,/defined)] $ ] IF IsType(bkgnd_radius,/defined) THEN BEGIN star_info = [star_info, 'noyzero : '+(['no','yes'])[noyzero] ] star_info = [star_info, 'bkgnd radius : '+strjoin(string(bkgnd_radius,format=format),' - ') ] IF IsType(bkgnd_origin,/defined) THEN $ star_info = [star_info, 'bkgnd origin : '+strjoin(string(bkgnd_origin,format=format),' - ') ] IF IsType(bkgnd_count,/defined) THEN $ star_info = [star_info, 'bkgnd count : '+strcompress(bkgnd_count,/rem) ] ENDIF IF IsType(wing_radius,/defined) THEN BEGIN star_info = [star_info, 'wing radius : '+string(wing_radius,format=format) ] IF IsType(wing_origin,/defined) THEN $ star_info = [star_info, 'wing origin : '+strjoin(string(wing_origin,format=format),' - ') ] star_info = [star_info, 'wing Ms-adjusted: '+(['no','yes'])[auto_wing]] ENDIF IF has_sigma THEN $ star_info = [star_info, 'sigma threshold : '+strcompress(sigma,/rem) ] IF silent LE 1 THEN FOR i=0,n_elements(star_info)-1 DO print, star_info[i] CASE 1 OF IsType(file,/string): BEGIN ; Name of sky map provided IF (file_search(file))[0] EQ '' THEN BEGIN count = 0 error = hide_env(file)+' not found' IF silent LE 1 THEN message, /info, error RETURN ENDIF InitVar, southpole, /key InitVar, northpole, /key equator = 1-southpole AND 1-northpole maptype = -1*southpole+1*northpole mapname = (['south pole','equatorial','north pole'])[maptype+1]+' map' smei_sky, file, southraw=southpole, northraw=northpole, equatraw=equator, $ hdr=hdr, sky=sky, nbin=1, /noplot, silent=silent+1, $ rmbkgnd=rmbkgnd, rmzld=rmzld IF size(sky,/n_dim) EQ 0 THEN BEGIN count = 0 final_sky = -1 flat_sky = -1 star_sky = -1 error = hide_env(file)+', '+mapname+' not present' IF silent LE 1 THEN message, /info, error RETURN ENDIF ndim = size(sky,/dim) stime = hdr.stime etime = hdr.etime camera = hdr.camera origin = [hdr.cx, hdr.cy] scale = 1.0/(hdr.bd*dpm) ; Read the PSF rotation, fovx and time map. IF IsType(psf_map,/undefined) THEN $ smei_sky, file, /psfn, sky=psf_map, /noplot, silent=silent+1, degrees=degrees IF IsType(fovx_map,/undefined) THEN $ smei_sky, file, /fovx, sky=fovx_map, /noplot, silent=silent+1, degrees=degrees IF IsType(time_map,/undefined) THEN BEGIN smei_sky, file, /orbsecs, sky=time_map, /noplot, silent=silent+1, exten_no=tmp torigin = TimeSet( fxpar(headfits(file,exten=tmp),'TORIGIN') ) ENDIF END IsType(hdr,/structure): BEGIN southpole = strpos(hdr.map,'Map of south pole') NE -1 northpole = strpos(hdr.map,'Map of north pole') NE -1 equator = strpos(hdr.map,'Equatorial map' ) NE -1 IF southpole+northpole+equator NE 1 THEN $ message, 'no hdr, hdr.map="'+hdr.map+'"' maptype = -1*southpole+0*equator+1*northpole mapname = (['south pole','equatorial','north pole'])[maptype+1]+' map' file = hdr.file IF IsType(sky,/undefined) THEN $ smei_sky, file, southraw=southpole, northraw=northpole, equatraw=equator, $ hdr=hdr, sky=sky, nbin=1, /noplot, silent=silent+1, $ rmbkgnd=rmbkgnd, rmzld=rmzld CASE size(sky,/n_dim) OF 0: BEGIN count = 0 final_sky = -1 flat_sky = -1 star_sky = -1 error = hide_env(file)+', '+mapname+' not present' IF silent LE 1 THEN message, /info, error RETURN END 2: ELSE: BEGIN count = 0 error = hide_env(file)+', skymap must be a 2D array' IF silent LE 1 THEN message, /info, error RETURN END ENDCASE ndim = size(sky,/dim) stime = hdr.stime etime = hdr.etime camera = hdr.camera origin = [hdr.cx, hdr.cy] scale = 1.0/(hdr.bd*dpm) ; Read the rotation, time and fovx map. ; Extract the scaling factors for the equatorial mapping from the header smei_sky, file, /psfn , sky=psf_map, /noplot, degrees=degrees, silent=silent+1 smei_sky, file, /fovx , sky=fovx_map, /noplot, degrees=degrees, silent=silent+1 smei_sky, file, /orbsecs, sky=time_map, /noplot, silent=silent+1, exten_no=tmp torigin = TimeSet( fxpar(headfits(file,exten=tmp),'TORIGIN') ) END ELSE: BEGIN InitVar, southpole, /key InitVar, northpole, /key equator = 1-southpole AND 1-northpole maptype = -1*southpole+1*northpole mapname = (['south polar','equatorial','north polar'])[maptype+1]+' map' IF size(sky, /n_dim) NE 2 THEN BEGIN count = 0 error = hide_env(file)+', skymap must be a 2-dimensional' IF silent LE 1 THEN message, /info, error RETURN ENDIF InitVar, camera, 1 ndim = size(sky,/dim) CASE equator OF 0: BEGIN InitVar, origin, 0.5*(ndim-1) InitVar, scale , (80.0/dpm)/ndim[0] END 1: BEGIN InitVar, origin, [-0.5,0.5*(ndim[1]-1)] InitVar, scale , (360.0/dpm)/ndim[0] END ENDCASE IF size(psf_map , /n_dim) NE 2 THEN message, 'psf map must be 2D array' IF size(fovx_map, /n_dim) NE 2 THEN message, 'fovx map must be 2D array' IF NOT IsTime(stime) THEN message, 'no time origin for time array specified' IF NOT IsTime(etime) THEN etime = TimeOp(/add,stime,smei_coriolis(stime,/orbital_period)) IF size(time_map, /n_dim) NE 2 THEN message, 'time map must be 2D array' IF NOT IsTime(torigin) THEN message, 'no time origin for time map specified' END ENDCASE mapcode = (['spmap','eqmap','npmap'])[maptype+1] IF silent LE 0 THEN $ message, /info, mapname+' for camera'+strcompress(camera) CASE IsType(star_list,/defined) OF 0: BEGIN ; No stars specified; built a list ; Read stars from SMEI catalogues ; rmstar=1 for stars on the bright star and duh-stars atalogue ; rmstar=0 for all other stars (nova_list, stars_not_subtracted) star_list = smei_star_list(silent=silent+1,catalogues=catalogues,rmstar=rmstar,count=count) ; Only stars brighter than the cutoff magnitude are removed ; from the skymaps. Fainter stars are fitted, but not removed. rmstar AND= smei_star_info(star_list,/get_mags) LE upto_mags ;================================= ; Add planets and asteroids to the list of stars ; Time halfway through orbit ; Only the brighter planets are subtracted. ; Asteroids are not subtracted. ut = TimeLimits([stime,etime],/mid) ssbody = ['Mercury','Venus','Mars','Jupiter','Saturn','Uranus','Neptune','Pluto'] nrbody = jpl_body(ssbody) mvbody = jpl_mag(ut,nrbody) rmbody = [ 1 , 1 , 1 , 1 , 1 , 1 , 0 , 0 ] ssastr = usno_body(number=nrastr,count=i) mvastr = replicate(upto_mags,i) ; Puts asteroids at end after sorting on magnitude rmastr = replicate(0 ,i) ; Fit, but do not subtract i = jpl_body(/earth) tmp = [[jpl_eph (ut,nrbody,center=i,/to_sphere,degrees=degrees,/location)], $ [usno_eph(ut,ssastr,center=i,/to_sphere,degrees=degrees,/location,/silent)] ] boost, ssbody, ssastr boost, nrbody, nrastr boost, mvbody, mvastr boost, rmbody, rmastr destroyvar, ssastr, nrastr, mvastr, rmastr ; Kinda clumsy: the star structure stores RA as hours,min,sec ; and decl as deg,min,sec, while AngleUnits returns ; msec in addition. tmp = AngleUnits(from_angle=tmp[0:1,*],/to_almanac,degrees=degrees,/singlesign) tmp = float(tmp) tmp[2,*,*] += tmp[3,*,*]/1000 ; Add msec to sec tmp = tmp[0:2,*,*] ; Drop msec body_list = replicate(star_list[0],n_elements(ssbody)) body_list.nr = -nrbody ; ID number (negative of JPL id) ; Make sure the name has the right number of chars ; (12) with trailing spaces. i = strlen(star_list[0].name) j = strjoin(replicate(' ',i)) body_list.name = strmid(ssbody+j,0,i) ; Planet name body_list.obafg = '' body_list.ra = reform(tmp[*,0,*]) ; RA body_list.dec = reform(tmp[*,1,*]) ; Planet ; Use visual magnitudes throughout. ; Probably not too good for SMEI mags. body_list.magb = mvbody body_list.magv = mvbody body_list.mags = mvbody boost, star_list, body_list boost, rmstar , rmbody destroyvar, ssbody, nrbody, mvbody, rmbody, body_list ;======================================== ; Sort on magnitude. ; At this point the bright planets will come floating to the top. i = sort(smei_star_info(star_list, /get_mags)) star_list = star_list[i] rmstar = rmstar [i] count = n_elements(star_list) ;tmp = where(smei_star_info(star_list,/get_mags) LE upto_mags, count) ;CASE count EQ 0 OF ;0: star_list = star_list[tmp] ;1: error = 'no stars brighter than SMEI magnitude'+strcompress(upto_mags) ;ENDCASE END 1: BEGIN count = n_elements(star_list) IF IsType(star_list, /string) THEN BEGIN star_list = smei_star_list(name=star_list, count=count, /silent, catalogues=catalogues) IF count EQ 0 THEN error = 'star name(s) not found in catalogues' ENDIF IF count NE 0 THEN InitVar, star_remove, replicate(1,count), set=rmstar END ENDCASE all_count = count ; count : # stars on star_list ; star_list [count] : list of SMEI_STAR_LIST structures ; rmstar [count] : 0/1: do not/do remove from skymap ; From here on 'good stars' are tracked by the index array 'valid' ; (with number of elements 'count'). This is an index array into ; 'star_list' (which is not modified anymore). ; Other arrays (rmstar, rastar, xystar, star_fit) will only have ; 'count' entries. IF count NE 0 THEN BEGIN ; Check for stars in/outside map ; Exclude stars that are not inside the map ; Pick up RA,dec for all 'count' stars rastar = smei_star_info(star_list, /radec, degrees=degrees) one = replicate(1.0,count) ; Translate RA,dec into array indices in the sky map xystar = cvsmei(from_equatorial=rastar, degrees=degrees,/to_map, mode=mapcode,/silent) ; Only keep stars inside map (valid array indices for centroid after roundoff) xx = round( reform(xystar[0,*]) ) yy = round( reform(xystar[1,*]) ) old_count = count valid = where( $ ; Subset of stars inside map 0 LE xx AND xx LT ndim[0] AND $ ; Valid horizontal index 0 LE yy AND yy LT ndim[1], $ ; Valid vertical index count ) CASE count NE 0 OF 0: error = 'no stars inside '+mapname 1: IF count LT old_count THEN $ IF silent LE 1 THEN $ message, /info, strcompress(old_count-count,/rem)+'/'+ $ strcompress(all_count,/rem)+' stars outside map' ENDCASE all_count = count ENDIF IF count NE 0 THEN BEGIN ; Don't fit stars for which there is no valid brightness and/or psf ; angle available. ; (not sure whether the test for brightness is such a good idea. ; Because it's done on the hires skymap, it throws out stars in ; relatively small 'holes' in the skymap for which a PSF is available ; in the lores psf map. ; Pick up intensity and psf angle for all stars inside map tmp = round(xystar[*,valid]) adu = reform(sky[tmp[0,*],tmp[1,*]]) ; Intensity at catalogue position ; Array indices on low-res map tmp = cvsmei(from_equatorial=rastar[*,valid],degrees=degrees,/to_map,mode='lores',/silent) tmp = round(tmp) psf = reform(psf_map[tmp[0,*],tmp[1,*]]); PSF angles ; Fit only stars with a finite intensity and a finite psf angle ; 'valid' is an array of indices to be used to access ; arrays star_fit, star_list[valid], rastar, xystar old_count = count tmp = where(finite(adu) AND finite(psf),count) CASE count NE 0 OF 0: error = 'no stars with valid intensity and psfangle at catalogue location' 1: BEGIN valid = valid[tmp] IF silent LE 1 THEN $ message, /info, strcompress(old_count-count,/rem)+'/'+ $ strcompress(all_count,/rem)+' without PSF values' END ENDCASE destroyvar, adu, psf ENDIF ; At this point, 'sky' contains a valid skymap. ; Initialize return array here to make sure they exist ; if no stars end up being fitted (count=0 returned) ; final_sky and flat_sky contain versions of the subtracted maps. ; star_sky is used to store the subtracted stars final_sky = sky star_sky = 0*sky IF count NE 0 THEN BEGIN timefmt = smei_star_formatpnt(/get_timefmt) ; 'valid' now refers to all stars in 'star_list' located inside map ; and with valid intensity and psf angle rmstar = rmstar[ valid] rastar = rastar[*,valid] xystar = xystar[*,valid] ; Create single structure for storing star fitting info ; and fill in as many fields as possible with default values ; prior to replicating it to an array for all stars InitVar, badfit, BadValue(0.0), set=tmp star_fit = {smei_star} star_fit.maptype = strmid(mapcode,0,2) star_fit.camera = camera star_fit.ccd_mode = ccd_mode star_fit.time = TimeGet(!TheTime,format=timefmt,/scalar) ; Time of observation star_fit.total_adus = tmp star_fit.bright_raw = tmp star_fit.bright_mdl = tmp star_fit.gain = 1.0 star_fit.name_close = '-' star_fit.bright = tmp ; Fitted brightness star_fit.back_const = tmp ; Fitted background star_fit.back_slope = tmp star_fit.n_bkgnd = 0 star_fit.n_fitpsf = 0 ; # Pixels used in fit star_fit.stdev_bkgnd= tmp star_fit.stdev_psf = tmp star_fit.cvfit = tmp star_fit.dradec = tmp star_fit.dpsfangle = tmp star_fit.dpsfstretch= tmp star_fit.psfwidth = tmp star_fit.remove = 1 star_fit.done = 0 ; Replicate structure to form array for all stars inside the map. star_fit = replicate(star_fit, count) star_fit.name = smei_star_info(star_list[valid], /get_name) star_fit.remove = rmstar star_fit.radec = rastar ; Pick up psf rotation angles for all stars inside map ; Pick up fovx direction cosine angles for all stars inside map ; If keyword nofovx is set then set the angle to zero. tmp = cvsmei(from_equatorial=rastar,degrees=degrees,/to_map,mode='lores',/silent) tmp = round(tmp) star_fit.psfangle = reform(psf_map[tmp[0,*],tmp[1,*]]) star_fit.fovangle = (1-nofovx)*[fovx_map[tmp[0,*],tmp[1,*]],dblarr(1,count)] star_fit.psfstretch = cos(star_fit.fovangle*rpm) ; Pick up times for all stars inside map with valid psf angles ; (TimeSet acts funny if a badvalue is input). time = reform(time_map[tmp[0,*],tmp[1,*]]) ; Times at star locations ; Fill in time of observation and gain. tmp = where(finite(time),complement=i) ; Stars with finite psfangle but no time ; This is not supposed to happen, but it does. ; If the psf angle is finite, then the time should be finite too. ; This probably needs to be fixed in the indexing program IF i[0] NE -1 THEN $ message, /info, 'stars with valid PSF angle but no time present' IF tmp[0] NE -1 THEN BEGIN i = TimeOp(/add,torigin,TimeSet(/difference,time[tmp],usec)) star_fit[tmp].time = TimeGet(i,format=timefmt,upto=usec,/roundt) star_fit[tmp].gain = smei_camera_gain(i, camera=star_fit[tmp].camera) ENDIF ; Freeze the brightnesses of stars in known_names to values in known_mags fixed_bright = replicate(BadValue(star_fit.bright),count) IF IsType(known_names,/defined) THEN BEGIN tmp = where_common(strtrim(star_fit.name),known_names) IF tmp[0] NE -1 THEN BEGIN fixed_bright[tmp] = smei_i2m(from_ms=known_mags) FOR i=0,n_elements(tmp)-1 DO $ message, /info, star_fit[tmp[i]].name + $ ' set to brightness'+strcompress(fixed_bright[tmp[i]]) ENDIF ENDIF ; std_ndim represents a 50x50 box. This matches the 200x200 standard star at ; one quarter the resolution ; (0.025 deg/bin in the standard star; 0.1 in the equatorial skymap). std_ndim = smei_star_standard(camera, $ ; Returns [5,5] (degrees) use_filled_psf = use_filled_psf, $ wing_radius = wing_radius , $ /get_fov , $ degrees=degrees) xtmp = round(std_ndim/scale) ; Makes [50,50] (skybins) std_ndim = cvsmei(from_equatorial=[[0,0],[std_ndim]],/to_map,mode='eqmap',degrees=degrees,/silent) std_ndim = round(std_ndim[*,1]-std_ndim[*,0]) ; Makes [50,50] (skybins) IF xtmp[0] NE std_ndim[0] OR xtmp[1] NE std_ndim[1] THEN message, 'stop std_ndim' CASE maptype EQ 0 OF 0: fullmap = n_elements(sky) EQ 800L* 800L ; North/south-polar map 1: fullmap = n_elements(sky) EQ 3600L*1200L ; Equatorial map ENDCASE view_star = fltarr(std_ndim[0],std_ndim[1]) corepsf_sky = smei_star_corepsf( $ fullmap , $ camera , $ maptype , $ origin , $ scale , $ size(sky,/dim) , $ std_ndim , $ ndim , $ count , $ lindgen(count) , $ star_list[valid] , $ star_fit , $ xystar , $ rastar , $ degrees = degrees , $ use_filled_psf = use_filled_psf, $ wing_radius = wing_radius , $ magnify = magnify*fix_something) IF show NE '' THEN $ IF NOT fullmap THEN $ view, in=(sky > (-100)) < 100,/stretch, /bc ; Estimated SMEI magnitude names = star_fit.name m_smei = smei_star_info(star_list[valid],/get_mags) ; Set up some strings for informational messages. cstar = names + $ strjoin(AngleUnits(from_angle=rastar,/to_almanac,/string,degrees=degrees),' ')+$ ' psf='+string(star_fit.psfangle ,format='(F7.2)') + $ ' tx=' +string(star_fit.fovangle[0],format='(F6.2)') + $ ' Ms=' +string(m_smei,format='(F5.2)') ;================================== ; Subtract known stars from star map ; know_stars[3,*]: ra, dec, smei-magnitude FOR ibright=0,n_elements(known_stars)/3-1 DO BEGIN m_known = known_stars[2,ibright] b_known = smei_i2m(from_ms=m_known) print, m_known, b_known raknown = known_stars[0:1,ibright] xyknown = cvsmei(from_equatorial=raknown, $ degrees=degrees,/to_map, mode=mapcode,/silent) ; Pick up coordinates for part of skymap around star as array ; indices into original sky map. ipix = smei_star_box( $ fullmap , $ maptype , $ std_ndim , $ ndim , $ xyknown , $ ibeg , $ iend , $ jbeg , $ jend , $ rr , $ ibeg0 = ibeg0 , $ jbeg0 = jbeg0 ) ; Shift origin to location of brightest star rr0 = rr-xyknown#replicate(1.0,(iend-ibeg+1)*(jend-jbeg+1)) ; The box is now (iend-ibeg+1) x (jend-jbeg+1) pixels ; Extract the box around the star to be fitted from the skymap original_star = final_sky[ipix,jbeg:jend] IF IsType(weights,/defined) THEN $ original_weights = weights[ipix,jbeg:jend] ; Extract the same box from the map storing the subtracted stars. ; This is passed to smei_star_stdmaps to exclude skybins that ; are part of an already-subtracted star to become part of the ; background in_bkgnd. ;IF NOT no_stars_mask THEN $ ; stars_mask = star_sky[ipix,jbeg:jend] ; corepsf_mask: ; 0: bins in core PSF of stars brighter than m_smei+2 ; 1: bins in background and core PSF of stars fainter than m_smei+2 IF NOT no_corepsf_mask THEN $ corepsf_mask = corepsf_sky[ipix,jbeg:jend] GT (m_known+2) ; Get the brightness at the locations in the standard star ; equivalent to locations rr. Note that if multiple stars are ; fitted then array standard_stars will have a trailing ; dimension equal to the number of stars fitted simultaneously. ; The PSF and background area are defined in in_fitpsf and ; in_bkgnd, based on the location of bins relative to the ; standard star. tmp = cvsmei(from_equatorial=raknown,degrees=degrees,/to_map,mode='lores',/silent) tmp = round(tmp) psfangle = reform(psf_map[tmp[0],tmp[1]]) psfstretch = cos((1-nofovx)*[fovx_map[tmp[0],tmp[1]],0.0d0]*rpm) IF auto_wing THEN $ wing_delta = scale*(3-m_known) standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin , $ scale , $ raknown , $ psfangle , $ psfstretch , $ degrees = degrees , $ in_fullpsf = in_fullpsf , $ in_fitpsf = in_fitpsf , $ in_corepsf = in_corepsf , $ in_bkgnd = in_bkgnd , $ std_rr = std_rr , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea , $ stdadus = stdadus ) view, in=MagnifyArray(original_star < 150,5), /stretch, /bl in_fullpsf = where(in_fullpsf) IF in_fullpsf[0] NE -1 THEN original_star[in_fullpsf] -= b_known*standard_stars[in_fullpsf] view, in=MagnifyArray(original_star < 150,5), /stretch, /br final_sky[ipix,jbeg:jend] = original_star ENDFOR ;================================== ENDIF flat_sky = final_sky cskip_dist = '<'+string(skip_dist,format='(F4.2)') cincl_dist = '<'+string(incl_dist,format='(F4.2)') ; The loop over all stars processes stars from the brightest ; stars down, i.e. in order of increasing magnitude max_crowd = 15 ; These are used only in the iteration processes used ; to refine some of the lsq values (as args to ; smei_star_enough_bins). min_psf_bins = 3 min_bkgnd_bins = 3 IF silent LE 0 THEN $ message, /info, strcompress(count,/rem)+' stars' FOR ibright=0,count-1 DO BEGIN ; Loop over all remaining stars IF star_fit[ibright].done NE 0 THEN $ ; Star already fitted; skip to next one continue IF silent LE 0 THEN BEGIN print print ENDIF ;=============== ; Star selection ; ; Construct a group of stars near ibright closer than ; incl_dist. The bright star itself is the first in the group. ; This group is fitted and subtracted as a group, unless ; kill_one is set. If /kill_one is set the group is fitted ; as a whole, but only the brightest star is subtracted, and ; the other stars will be fitted separately later. n_crowd = 1 ; Nr of stars in group crowd = ibright ; Subset of valid stars ; Look for nearby stars that have not been subtracted yet. ; Exclude the bright star itself not_done = star_fit.done EQ 0 not_done[crowd] = 0 member = 0 REPEAT BEGIN ; Loop over members ; (increasing elongation to ibright) boss = crowd[member] near = where(not_done, n_near) IF n_near GT 0 THEN BEGIN ; 'grunts' are all stars that have not been subtracted ; yet, excluding stars already member of the group. grunts = near elo = sphere_distance(rastar[*,grunts],rastar[*,boss],degrees=degrees) ; Grunts within incl_dist near = where(elo LT incl_dist, n_near) IF n_near GT 0 THEN BEGIN grunts = grunts[near] ; Sort on elongation to the brightest star elo = sphere_distance(rastar[*,grunts],rastar[*,ibright],degrees=degrees) grunts = grunts[sort(elo)] ; Make sure not to exceed the maximum group size max_crowd. i = max_crowd-n_crowd ; i stars to fill up group IF n_near GT i THEN BEGIN grunts = grunts[0:i-1] n_near = i ENDIF ; Extend group with nearby 'grunts' ; Update the not_done flag for the new members not_done[grunts] = 0 crowd = [crowd,grunts] n_crowd += n_near ENDIF ENDIF ; Members in crowd are ordered on increasing ; elongation from the brightest star. member += 1 ; Next member ENDREP UNTIL member GE n_crowd OR n_crowd GE max_crowd ; Sort 'crowd' on brightness. ; The brightest star ibright should stay in first position. crowd = crowd[sort(crowd)] ; Now trim down the group by removing stars that are within ; 'skip_dist' of a brighter star in the group. These are too ; close to the brighter stars to be separated from it. ; The fainter star will be marked as 'done', and is never fitted. ; Note that star_fit[crowd].done = 0 for all members, ; so don't have to test for it. member = 0 ; Loop over members (decreasing brightness) WHILE member LT n_crowd-1 DO BEGIN boss = crowd[member] ; Stars in group fainter than 'boss' grunts = crowd[member+1:*] elo = sphere_distance(rastar[*,grunts],rastar[*,boss],degrees=degrees) ; Stars fainter than 'boss' and closer than skip_dist near = where(elo LT skip_dist, n_near, complement=far, ncomplement=n_far) IF n_near GT 0 THEN BEGIN ; Stars to be removed from group ; Mark the 'done' and 'name_close' fields in star_fit star_fit[grunts[near]].done = 1 ; Mark fainter star as 'done' star_fit[grunts[near]].name_close = names[boss] FOR i=0,(silent LE 1)*n_near-1 DO $ message, /info, cstar[grunts[near[i]]] + $ ' Elo='+string(elo[near[i]]*dpm,format='(F4.2)') + $ cskip_dist+' -'+strtrim(names[boss]) ; Removing the nearby faint stars from the group CASE n_far GT 0 OF 0: crowd = crowd[0:member] 1: crowd = [crowd[0:member],grunts[far]] ENDCASE n_crowd -= n_near ENDIF member += 1 ENDWHILE ; 'crowd' now contains the group of remaining stars ; sorted in order of decreasing brightness. ; If the group is bigger than the max allowed group size ; of max_stars, then shrink the group by removing stars ; farthest away from ibright. elo = sphere_distance(rastar[*,crowd],rastar[*,ibright],degrees=degrees) IF n_crowd GT max_stars THEN BEGIN ; Sort in order of increasing elongations member = sort(elo) crowd = crowd[member] elo = elo[member] FOR member=max_stars,(silent LE 1)*n_crowd-1 DO $ message, /info, cstar[crowd[member]] + $ ' Elo='+string(elo[member]*dpm,format='(F4.2)') + $ ' K' n_crowd = max_stars crowd = crowd[0:n_crowd-1] elo = elo[0:n_crowd-1] ; Sort back into order of decreasing brightness member = sort(crowd) crowd = crowd[member] elo = elo[member] ENDIF FOR member=0,(silent LE 1)*n_crowd-1 DO $ message, /info, cstar[crowd[member]] + $ ' Elo='+string(elo[member]*dpm,format='(F4.2)') ;=============== ; 'crowd' now contains a subset of entries from 'valid' ; These are the stars to be fitted simultaneously. ; Note that crowd[0] = ibright one_crowd = replicate(1,n_crowd) ; Pick up coordinates for part of skymap around star as array ; indices into original sky map. ipix = smei_star_box( $ fullmap , $ maptype , $ std_ndim , $ ndim , $ xystar[*,ibright] , $ ibeg , $ iend , $ jbeg , $ jend , $ rr , $ ibeg0 = ibeg0 , $ jbeg0 = jbeg0 ) ; Shift origin to location of brightest star rr0 = rr-xystar[*,ibright]#replicate(1.0,(iend-ibeg+1)*(jend-jbeg+1)) ; The box is now (iend-ibeg+1) x (jend-jbeg+1) pixels ; Extract the box around the star to be fitted from the skymap original_star = final_sky[ipix,jbeg:jend] original_flat = flat_sky [ipix,jbeg:jend] IF IsType(weights,/defined) THEN $ original_weights = weights[ipix,jbeg:jend] ; Extract the same box from the map storing the subtracted stars. ; This is passed to smei_star_stdmaps to exclude skybins that ; are part of an already-subtracted star to become part of the ; background in_bkgnd. ;IF NOT no_stars_mask THEN $ ; stars_mask = star_sky[ipix,jbeg:jend] ; corepsf_mask: ; 0: bins in core PSF of stars brighter than m_smei+2 ; 1: bins in background and core PSF of stars fainter than m_smei+2 IF NOT no_corepsf_mask THEN $ corepsf_mask = corepsf_sky[ipix,jbeg:jend] GT (m_smei[ibright]+2) ; Get the brightness at the locations in the standard star ; equivalent to locations rr. Note that if multiple stars are ; fitted then array standard_stars will have a trailing ; dimension equal to the number of stars fitted simultaneously. ; The PSF and background area are defined in in_fitpsf and ; in_bkgnd, based on the location of bins relative to the ; standard star. psfangle = star_fit[crowd].psfangle psfstretch = cos(star_fit[crowd].fovangle*rpm) ; wing_delta is added to wing_radius in smei_star_standard, ; and decreases the effective wing_radius by 1 bin (0.1 deg) ; per unit of magnitude. IF auto_wing THEN $ wing_delta = scale*(3-m_smei[crowd]) standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin , $ scale , $ rastar [*,crowd] , $ psfangle , $ psfstretch , $ degrees = degrees , $ in_fullpsf = in_fullpsf , $ in_fitpsf = in_fitpsf , $ in_corepsf = in_corepsf , $ in_bkgnd = in_bkgnd , $ std_rr = std_rr , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea , $ stdadus = stdadus ) ; Refine the PSF and background areas, based on properties of the ; original star and its environment. smei_star_lsqbins,original_star,rr,rr0,in_fitpsf,in_bkgnd,bkgnd_count=bkgnd_count in_fullpsf = where(in_fullpsf, n_fullpsf) in_fitpsf = where(in_fitpsf , n_fitpsf ) in_bkgnd = where(in_bkgnd , n_bkgnd ) ;IF n_fitpsf GE 25 AND (1-fit_bkgnd OR n_bkgnd GE 5) THEN BEGIN IF smei_star_enough_bins(fit_bkgnd, n_fitpsf, n_bkgnd) THEN BEGIN only = -1 star_model = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf , $ in_fitpsf , $ in_bkgnd , $ back , $ bright , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model , $ one_star_model = one_star_model, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights,$ fixed_bright = fixed_bright[crowd]) cvfit = correlate((bkgnd_model+star_model)[in_fitpsf],original_star[in_fitpsf]) ; Corrections for parameters fitted iteratively ; (set by iteration loop below) ; Corrections from first round of iterations: ; dxystar, dpsangle and dpsfstretch. ; dxystar is determined by moving the n_crowd stars around ; as a group keeping relative locations fixed. This will ; primarily catch pointing errors which affect all stars ; in the group the same way (if it doesn't than this ; should be caught in the second round of iterations, ; see below). ; dradec is determined from dxystar. dxystar = [0.0,0.0] dpsfangle = 0.0 dpsfstretch = [0.0,0.0] dradec = [0.0,0.0] ; Corrections from second round of iterations: ; here the relative locations of the stars are allowed ; to vary individually. This should catch small errors ; in the catalogue locations. ddxystar = dxystar#one_crowd ddradec = dradec #one_crowd ; Also fit the centroid and/or psf width of the brightest star IF fix_something THEN BEGIN ; Max number of parameters fitted. ; Currently x,y of centroid and fovx and fovy angles nfix = 5 ; Array index in iteration vectors (0,..,nfix-1) pos_centroid = [0,1] ; x,y of centroid pos_psfangle = 2 ; psf angle pos_psfstretch = [3,4] ; fovx and fovy angle pos = [pos_centroid,pos_psfangle,pos_psfstretch] ; Differential correction vector dfix_best = fltarr(nfix) dfix_best[pos] = [dxystar,dpsfangle,dpsfstretch] ; Set size of neighbourhood to be explored in each iteration ; Should be one for each parameter fitted ; Step size for centroid fitting in x and y direction ; (in skybins) dstep_centroid = [0.10,0.10]*fix_centroid ; Step size for psfangle fitting and for psfstretch fitting in long ; and narrow dimension (angular units) dstep_psfangle = 0.01/dpm*fix_psfangle dstep_psfstretch = [0.001,0.001]*fix_fovangle dstep = fltarr(nfix) dstep[pos] = [dstep_centroid,dstep_psfangle,dstep_psfstretch] ; Indices in dfix_best,dstep of parameters that need fixing fit_bfix = dstep NE 0.0 fit_nfix = where(fit_bfix) ; Left and right side indices in cv (cv[0] is center value) left = 1+2*indgen(n_elements(fit_nfix)) right = left+1 ; Set up dsteps (array[nfix,2*nfix]) array defining ; neighbourhood around current best fit. For each parameter ; pick points left and right of the current best fit. ; The following creates an array of structure: ; dsteps = [[-a, 0, 0, 0] , $ ; [ a, 0, 0, 0] , $ ; [ 0,-a, 0, 0] , $ ; [ 0, a, 0, 0] , $ ; [ 0, 0,-b, 0] , $ ; [ 0, 0, b, 0] , $ ; [ 0, 0, 0,-b] , $ ; [ 0, 0, 0, b] ] dsteps = fltarr(nfix,2*nfix) ; All zeroes FOR i=0,nfix-1 DO $ dsteps[pos[i],2*pos[i]+[0,1]] = dstep[pos[i]]*[-1,1] ; Initialize "newbest" values cv_newbest = cvfit ; Best fit so far dfix_newbest = dfix_best ; Zeroes max_pass = 400 pass = 0 fixed = 0 ; This switch is used to abort the iteration process, if ; something goes really wrong (i.e. if smei_star_stdmaps ; returns a near-zero number of sky bins for the lsq fit. just_died = 0 REPEAT BEGIN ; The hard way: repeat until convergence ; Save best values so far cv_best = cv_newbest dfix_best = dfix_newbest cv = cv_best ; Used to accumulate corr coeff FOR i=0,2*nfix-1 DO BEGIN ; Check neigbourhood (left, right) IF NOT fit_bfix[i/2] THEN continue dfix_tmp = dfix_best+dsteps[*,i] origin_tmp = origin +dfix_tmp[pos_centroid ] psfangle_tmp = psfangle +dfix_tmp[pos_psfangle ]*one_crowd psfstretch_tmp = psfstretch+dfix_tmp[pos_psfstretch]#one_crowd standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin_tmp , $ ; Iterate on centroid location scale , $ rastar[*,crowd] , $ psfangle_tmp , $ ; Iterate on PSF angle psfstretch_tmp , $ ; Iterate on PSF width degrees = degrees , $ in_fullpsf = in_fullpsf_tmp, $ in_fitpsf = in_fitpsf_tmp , $ in_corepsf = in_corepsf_tmp, $ in_bkgnd = in_bkgnd_tmp , $ std_rr = std_rr_tmp , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea_tmp , $ stdadus = stdadus_tmp ) smei_star_lsqbins,original_star,rr,rr0,in_fitpsf_tmp,in_bkgnd_tmp,bkgnd_count=bkgnd_count in_fullpsf_tmp = where(in_fullpsf_tmp, m_fullpsf) in_fitpsf_tmp = where(in_fitpsf_tmp , m_fitpsf ) ;in_corepsf_tmp = where(in_corepsf_tmp, m_corepsf) in_bkgnd_tmp = where(in_bkgnd_tmp , m_bkgnd ) just_died = 1-smei_star_enough_bins( $ fit_bkgnd, m_fitpsf, m_bkgnd, $ min_psf_bins = min_psf_bins , $ min_bkgnd_bins = min_bkgnd_bins) IF just_died THEN break ; Breaks out off the neighbourhood FOR loop only = -1 star_model_tmp = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf_tmp , $ in_fitpsf_tmp , $ in_bkgnd_tmp , $ back_tmp , $ bright_tmp , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model_tmp , $ one_star_model = one_star_model_tmp, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights , $ fixed_bright = fixed_bright[crowd]) ; Correlation skymap vs model star cv_tmp = correlate((bkgnd_model_tmp+star_model_tmp)[in_fitpsf_tmp],original_star[in_fitpsf_tmp]) ; Save the best correlation IF cv_tmp GT cv_best THEN BEGIN cv_newbest = cv_tmp dfix_newbest = dfix_best+dsteps[*,i] fixed = 1 ENDIF ; Accumulate correlation coefficients cv = [cv,cv_tmp] ENDFOR ; This happens if 2*nfix neighbourhood points ; does not have enough points for an lsq fit anymore. ; In this case we just use the values used to initialize ; the iteration process. IF just_died THEN break ; Breaks out off the iteration REPEAT loop tmp = max(cv,i_max) cv_left = cv[left ]-cv[0] cv_right = cv[right]-cv[0] ; Fit parabole in each of the parameters to improve ; the location of maximum correlation top = 0.5*(cv_left-cv_right)/(cv_left+cv_right) top = (top < 2.0) > (-2.0) ntop = 1 IF i_max NE 0 THEN BEGIN ; At least one of the sides is higher than the center. ; Step closer to solution using a steepest ascent estimate ; If the center correlation cv[0] is lower (or higher) ; than both sides, then set the difference between both ; sides to zero (this presumably means that the best ; fit in that dimension has been found, so we don't want ; to step in that direction). ; Find direction of steepest ascent uphill = (cv_right-cv_left)*((cv_right GT 0.0) NE (cv_left GT 0.0)) IF total(uphill*uphill) NE 0 THEN BEGIN uphill /= sqrt(total(uphill*uphill)) ; Unit vector top = [[top], [0.5*uphill], [5*uphill]] ntop = 3 ENDIF ENDIF FOR i=0,ntop-1 DO BEGIN dfix_top = dfix_best dfix_top[fit_nfix] += top[*,i]*dstep[fit_nfix];?? psfangle_tmp = psfangle +dfix_top[pos_psfangle ]*one_crowd psfstretch_tmp = psfstretch+dfix_top[pos_psfstretch]#one_crowd standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin+dfix_top[pos_centroid] , $ scale , $ rastar[*,crowd] , $ psfangle_tmp , $ psfstretch_tmp , $ degrees = degrees , $ in_fullpsf = in_fullpsf_tmp, $ in_fitpsf = in_fitpsf_tmp , $ in_corepsf = in_corepsf_tmp, $ in_bkgnd = in_bkgnd_tmp , $ std_rr = std_rr_tmp , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea_tmp , $ stdadus = stdadus_tmp ) smei_star_lsqbins,original_star,rr,rr0,in_fitpsf_tmp,in_bkgnd_tmp,bkgnd_count=bkgnd_count in_fullpsf_tmp = where(in_fullpsf_tmp, m_fullpsf) in_fitpsf_tmp = where(in_fitpsf_tmp , m_fitpsf ) ;in_corepsf_tmp = where(in_corepsf_tmp, m_corepsf) in_bkgnd_tmp = where(in_bkgnd_tmp , m_bkgnd ) just_died = 1-smei_star_enough_bins( $ fit_bkgnd, m_fitpsf, m_bkgnd, $ min_psf_bins = min_psf_bins , $ min_bkgnd_bins = min_bkgnd_bins) IF just_died THEN BEGIN ; Skip the corresponding step just_died = 0 ; Don't want to abort because of this continue ENDIF only = -1 star_model_tmp = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf_tmp , $ in_fitpsf_tmp , $ in_bkgnd_tmp , $ back_tmp , $ bright_tmp , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model_tmp , $ one_star_model = one_star_model_tmp, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights , $ fixed_bright = fixed_bright[crowd]) ; Correlation star model and skymap cv_tmp = correlate((bkgnd_model_tmp+star_model_tmp)[in_fitpsf_tmp],original_star[in_fitpsf_tmp]) IF cv_tmp GT cv_best THEN BEGIN cv_newbest = cv_tmp dfix_newbest = dfix_top fixed = 1 ENDIF ENDFOR ; Stop if there is no improvement, or if the iteration limit is reached pass += 1 ENDREP UNTIL pass EQ max_pass OR (where(dfix_newbest NE dfix_best))[0] EQ -1 IF silent LE 0 THEN $ message, /info, strcompress(pass,/rem)+'/'+strcompress(max_pass,/rem)+' iterations' IF pass EQ max_pass THEN $ message, /info, 'iteration limit reached for centroid fix; no convergence' IF just_died THEN BEGIN IF silent LE 2 THEN $ message, /info, 'iteration aborted: not enough bins in neighbourhood point(s)' ENDIF ELSE BEGIN psfangle_tmp = psfangle +dfix_best[pos_psfangle ]*one_crowd psfstretch_tmp = psfstretch+dfix_best[pos_psfstretch]#one_crowd standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin+dfix_best[pos_centroid] , $ scale , $ rastar[*,crowd] , $ psfangle_tmp , $ psfstretch_tmp , $ degrees = degrees , $ in_fullpsf = in_fullpsf_tmp, $ in_fitpsf = in_fitpsf_tmp , $ in_corepsf = in_corepsf_tmp, $ in_bkgnd = in_bkgnd_tmp , $ std_rr = std_rr_tmp , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea_tmp , $ stdadus = stdadus_tmp ) smei_star_lsqbins,original_star,rr,rr0,in_fitpsf_tmp,in_bkgnd_tmp,bkgnd_count=bkgnd_count in_fullpsf_tmp = where(in_fullpsf_tmp, m_fullpsf) in_fitpsf_tmp = where(in_fitpsf_tmp , m_fitpsf ) ;in_corepsf_tmp = where(in_corepsf_tmp, m_corepsf) in_bkgnd_tmp = where(in_bkgnd_tmp , m_bkgnd ) ; Not likely to happen at this point just_died = 1-smei_star_enough_bins( $ fit_bkgnd, m_fitpsf, m_bkgnd, $ min_psf_bins = min_psf_bins , $ min_bkgnd_bins = min_bkgnd_bins) IF just_died AND silent LE 2 THEN $ message, /info, 'iteration aborted: step determination for next iteration failed' ENDELSE IF NOT just_died THEN BEGIN only = -1 star_model_tmp = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf_tmp , $ in_fitpsf_tmp , $ in_bkgnd_tmp , $ back_tmp , $ bright_tmp , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model_tmp , $ one_star_model = one_star_model_tmp, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights , $ fixed_bright = fixed_bright[crowd]) cv = correlate((bkgnd_model_tmp+star_model_tmp)[in_fitpsf_tmp],original_star[in_fitpsf_tmp]) IF fixed AND cv LT cvfit THEN message, 'oops: worse' IF cv GT cvfit THEN BEGIN ; If the iteration process above was successfull, this is ; where the final values for the best fit are copied from ; the temp array into the final arrays, replacing the ; values from the initial values set at the start of the ; iterations. cvfit = cv star_model = star_model_tmp bkgnd_model = bkgnd_model_tmp one_star_model = one_star_model_tmp in_fullpsf = in_fullpsf_tmp in_fitpsf = in_fitpsf_tmp in_corepsf = in_corepsf_tmp in_bkgnd = in_bkgnd_tmp std_rr = std_rr_tmp back = back_tmp bright = bright_tmp stdarea = stdarea_tmp stdadus = stdadus_tmp dxystar = dfix_best[pos_centroid] dpsfangle = dfix_best[pos_psfangle] dpsfstretch = dfix_best[pos_psfstretch] ; The vector dxystar is the shift in centroid in skymap ; pixels. Convert to a shift in RA,dec using the inverse ; of the map projection. ; xystar+dxystar in the map corresponds to rastar+dradec in the sky dradec = cvsmei( $ from_map= xystar[*,ibright]+dxystar , $ mode = mapcode , $ /to_equatorial , $ degrees = degrees , $ /silent) - rastar[*,ibright] ; RA,dec shift for whole groupd dradec[0] = AngleRange(dradec[0],/pi,degrees=degrees) ddxystar = dxystar#one_crowd ddradec = dradec #one_crowd ; RA,dec shifts for individual stars ;IF total(dxystar*dxystar) GT 0.25 THEN $ ; IF silent LE 0 THEN $ ; message, /info, 'bad quaternion ??' IF in_fullpsf[0] NE -1 THEN n_fullpsf = n_elements(in_fullpsf) ELSE n_fullpsf = 0 IF in_fitpsf [0] NE -1 THEN n_fitpsf = n_elements(in_fitpsf ) ELSE n_fitpsf = 0 IF in_bkgnd [0] NE -1 THEN n_bkgnd = n_elements(in_bkgnd ) ELSE n_bkgnd = 0 ENDIF ENDIF ; If multiple stars are fitted simultaneously, and fix_pattern is set ; then start a second iteration process to refine the relative ; locations of all the stars. IF NOT just_died AND fix_pattern AND n_crowd GT 1 THEN BEGIN ; Number of parameters fitted: x,y for each star nfix = 2*n_crowd ; Array index in iteration vectors (0,..,nfix-1) pos_pattern = indgen(2*n_crowd) ; x,y of centroid pos = pos_pattern ; Differential correction vector dfix_best = fltarr(nfix) dfix_best[pos] = ddxystar ; Set size of neighbourhood to be explored in each iteration ; Should be one for each parameter fitted ; Step size for centroid fitting in x and y direction ; (in skybins) dstep_pattern = replicate(0.01,2*n_crowd) dstep = fltarr(nfix) dstep[pos] = dstep_pattern ; Indices in dfix_best,dstep of parameters that need fixing fit_bfix = dstep NE 0.0 fit_nfix = where(fit_bfix) ; Left and right side indices in cv (cv[0] is center value) left = 1+2*indgen(n_elements(fit_nfix)) right = left+1 ; Set up dsteps (array[nfix,2*nfix]) array defining ; neighbourhood around current best fit. For each parameter ; pick points left and right of the current best fit. ; The following creates an array of structure: ; dsteps = [[-a, 0, 0, 0] , $ ; [ a, 0, 0, 0] , $ ; [ 0,-a, 0, 0] , $ ; [ 0, a, 0, 0] , $ ; [ 0, 0,-b, 0] , $ ; [ 0, 0, b, 0] , $ ; [ 0, 0, 0,-b] , $ ; [ 0, 0, 0, b] ] dsteps = fltarr(nfix,2*nfix) ; All zeroes FOR i=0,nfix-1 DO $ dsteps[pos[i],2*pos[i]+[0,1]] = dstep[pos[i]]*[-1,1] ; Initialize "newbest" values cv_newbest = cvfit ; Best fit so far dfix_newbest = dfix_best ; Copy of ddxystar max_pass = 400 pass = 0 fixed = 0 xystar_best = xystar[*,crowd] psfangle_best = psfangle +dpsfangle *one_crowd psfstretch_best = psfstretch+dpsfstretch#one_crowd REPEAT BEGIN ; The hard way: repeat until convergence ; Save best values so far cv_best = cv_newbest dfix_best = dfix_newbest cv = cv_best ; Used to accumulate corr coeff FOR i=0,2*nfix-1 DO BEGIN ; Check neigbourhood (left, right) IF NOT fit_bfix[i/2] THEN continue dfix_tmp = dfix_best+dsteps[*,i] rastar_tmp = cvsmei( $ from_map= xystar_best+dfix_tmp , $ mode = mapcode , $ /to_equatorial , $ degrees = degrees , $ /silent) standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin , $ scale , $ rastar_tmp , $ ; Iterate on position psfangle_best , $ psfstretch_best , $ degrees = degrees , $ in_fullpsf = in_fullpsf_tmp, $ in_fitpsf = in_fitpsf_tmp , $ in_corepsf = in_corepsf_tmp, $ in_bkgnd = in_bkgnd_tmp , $ std_rr = std_rr_tmp , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea_tmp , $ stdadus = stdadus_tmp ) smei_star_lsqbins,original_star,rr,rr0,in_fitpsf_tmp,in_bkgnd_tmp,bkgnd_count=bkgnd_count in_fullpsf_tmp = where(in_fullpsf_tmp, m_fullpsf) in_fitpsf_tmp = where(in_fitpsf_tmp , m_fitpsf ) ;in_corepsf_tmp = where(in_corepsf_tmp, m_corepsf) in_bkgnd_tmp = where(in_bkgnd_tmp , m_bkgnd ) just_died = 1-smei_star_enough_bins( $ fit_bkgnd, m_fitpsf, m_bkgnd, $ min_psf_bins = min_psf_bins , $ min_bkgnd_bins = min_bkgnd_bins) IF just_died THEN break ; Breaks out off the neighbourhood FOR loop only = -1 star_model_tmp = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf_tmp , $ in_fitpsf_tmp , $ in_bkgnd_tmp , $ back_tmp , $ bright_tmp , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model_tmp , $ one_star_model = one_star_model_tmp, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights , $ fixed_bright = fixed_bright[crowd]) ; Correlation skymap vs model star cv_tmp = correlate((bkgnd_model_tmp+star_model_tmp)[in_fitpsf_tmp],original_star[in_fitpsf_tmp]) ; Save the best correlation IF cv_tmp GT cv_best THEN BEGIN cv_newbest = cv_tmp dfix_newbest = dfix_best+dsteps[*,i] fixed = 1 ENDIF ; Accumulate correlation coefficients cv = [cv,cv_tmp] ENDFOR ; This happens if 2*nfix neighbourhood points ; does not have enough points for an lsq fit anymore. ; In this case we just use the values used to initialize ; the iteration process. IF just_died THEN break ; Breaks out off the iteration REPEAT loop tmp = max(cv,i_max) cv_left = cv[left ]-cv[0] cv_right = cv[right]-cv[0] ; Fit parabole in each of the parameters to improve ; the location of maximum correlation top = 0.5*(cv_left-cv_right)/(cv_left+cv_right) top = (top < 2.0) > (-2.0) ntop = 1 IF i_max NE 0 THEN BEGIN ; At least one of the sides is higher than the center. ; Step closer to solution using a steepest ascent estimate ; If the center correlation cv[0] is lower (or higher) ; than both sides, then set the difference between both ; sides to zero (this presumably means that the best ; fit in that dimension has been found, so we don't want ; to step in that direction). ; Find direction of steepest ascent uphill = (cv_right-cv_left)*((cv_right GT 0.0) NE (cv_left GT 0.0)) IF total(uphill*uphill) NE 0 THEN BEGIN uphill /= sqrt(total(uphill*uphill)) ; Unit vector top = [[top], [0.5*uphill], [5*uphill]] ntop = 3 ENDIF ENDIF FOR i=0,ntop-1 DO BEGIN dfix_top = dfix_best dfix_top[fit_nfix] += top[*,i]*dstep[fit_nfix];??? rastar_tmp = cvsmei( $ from_map= xystar_best+dfix_top , $ mode = mapcode , $ /to_equatorial , $ degrees = degrees , $ /silent) standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin , $ scale , $ rastar_tmp , $ psfangle_best , $ psfstretch_best , $ degrees = degrees , $ in_fullpsf = in_fullpsf_tmp, $ in_fitpsf = in_fitpsf_tmp , $ in_corepsf = in_corepsf_tmp, $ in_bkgnd = in_bkgnd_tmp , $ std_rr = std_rr_tmp , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea_tmp , $ stdadus = stdadus_tmp ) smei_star_lsqbins,original_star,rr,rr0,in_fitpsf_tmp,in_bkgnd_tmp,bkgnd_count=bkgnd_count in_fullpsf_tmp = where(in_fullpsf_tmp, m_fullpsf) in_fitpsf_tmp = where(in_fitpsf_tmp , m_fitpsf ) ;in_corepsf_tmp = where(in_corepsf_tmp, m_corepsf) in_bkgnd_tmp = where(in_bkgnd_tmp , m_bkgnd ) just_died = 1-smei_star_enough_bins( $ fit_bkgnd, m_fitpsf, m_bkgnd, $ min_psf_bins = min_psf_bins , $ min_bkgnd_bins = min_bkgnd_bins) IF just_died THEN BEGIN just_died = 0 ; Don't want to abort because of this continue ; Skip the corresponding step ENDIF only = -1 star_model_tmp = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf_tmp , $ in_fitpsf_tmp , $ in_bkgnd_tmp , $ back_tmp , $ bright_tmp , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model_tmp , $ one_star_model = one_star_model_tmp, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights , $ fixed_bright = fixed_bright[crowd]) ; Correlation star model and skymap cv_tmp = correlate((bkgnd_model_tmp+star_model_tmp)[in_fitpsf_tmp],original_star[in_fitpsf_tmp]) IF cv_tmp GT cv_best THEN BEGIN cv_newbest = cv_tmp dfix_newbest = dfix_top fixed = 1 ENDIF ENDFOR ; Stop if there is no improvement, or if the iteration limit is reached pass += 1 ENDREP UNTIL pass EQ max_pass OR (where(dfix_newbest NE dfix_best))[0] EQ -1 IF silent LE 0 THEN $ message, /info, strcompress(pass,/rem)+'/'+strcompress(max_pass,/rem)+' iterations' IF pass EQ max_pass THEN $ message, /info, 'iteration limit reached for pattern fix; no convergence' IF just_died THEN BEGIN IF silent LE 2 THEN $ message, /info, 'pattern fix iteration aborted: not enough bins in neighbourhood point(s)' ENDIF ELSE BEGIN rastar_tmp = cvsmei( $ from_map= xystar_best+dfix_best , $ mode = mapcode , $ /to_equatorial , $ degrees = degrees , $ /silent) standard_stars = smei_star_stdmaps( $ camera , $ maptype , $ rr , $ origin , $ scale , $ rastar_tmp , $ psfangle_best , $ psfstretch_best , $ degrees = degrees , $ in_fullpsf = in_fullpsf_tmp, $ in_fitpsf = in_fitpsf_tmp , $ in_corepsf = in_corepsf_tmp, $ in_bkgnd = in_bkgnd_tmp , $ std_rr = std_rr_tmp , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ use_filled_psf = use_filled_psf, $ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea_tmp , $ stdadus = stdadus_tmp ) smei_star_lsqbins,original_star,rr,rr0,in_fitpsf_tmp,in_bkgnd_tmp,bkgnd_count=bkgnd_count in_fullpsf_tmp = where(in_fullpsf_tmp, m_fullpsf) in_fitpsf_tmp = where(in_fitpsf_tmp , m_fitpsf ) ;in_corepsf_tmp = where(in_corepsf_tmp, m_corepsf) in_bkgnd_tmp = where(in_bkgnd_tmp , m_bkgnd ) ; Not likely to happen at this point just_died = 1-smei_star_enough_bins( $ fit_bkgnd, m_fitpsf, m_bkgnd, $ min_psf_bins = min_psf_bins , $ min_bkgnd_bins = min_bkgnd_bins) IF just_died AND silent LE 2 THEN $ message, /info, 'pattern fix iteration aborted: step determination for next iteration failed' ENDELSE IF NOT just_died THEN BEGIN only = -1 star_model_tmp = smei_star_lsqfit( $ rr0 , $ standard_stars , $ original_star , $ in_fullpsf_tmp , $ in_fitpsf_tmp , $ in_bkgnd_tmp , $ back_tmp , $ bright_tmp , $ sigma = sigma , $ only = only , $ noplane = noplane , $ bkgnd_model = bkgnd_model_tmp , $ one_star_model = one_star_model_tmp, $ noyzero = noyzero , $ use_weights = use_weights , $ original_weights= original_weights , $ fixed_bright = fixed_bright[crowd]) cv = correlate((bkgnd_model_tmp+star_model_tmp)[in_fitpsf_tmp],original_star[in_fitpsf_tmp]) IF fixed AND cv LT cvfit THEN message, 'oops: worse in pattern fix' IF cv GT cvfit THEN BEGIN ; If the iteration process triggered by setting ; fix_patter above was successfull, this is ; where the final values for the best fit are copied from ; the temp array into the final arrays, replacing the ; values from first iteration process. cvfit = cv star_model = star_model_tmp bkgnd_model = bkgnd_model_tmp one_star_model = one_star_model_tmp in_fullpsf = in_fullpsf_tmp in_fitpsf = in_fitpsf_tmp in_corepsf = in_corepsf_tmp in_bkgnd = in_bkgnd_tmp std_rr = std_rr_tmp back = back_tmp bright = bright_tmp stdarea = stdarea_tmp stdadus = stdadus_tmp ddxystar = reform(dfix_best,2,n_crowd) ddradec = cvsmei( $ from_map= xystar[*,crowd]+ddxystar , $ mode = mapcode , $ /to_equatorial , $ degrees = degrees , $ /silent)-rastar[*,crowd] IF in_fullpsf[0] NE -1 THEN n_fullpsf = n_elements(in_fullpsf) ELSE n_fullpsf = 0 IF in_fitpsf [0] NE -1 THEN n_fitpsf = n_elements(in_fitpsf ) ELSE n_fitpsf = 0 IF in_bkgnd [0] NE -1 THEN n_bkgnd = n_elements(in_bkgnd ) ELSE n_bkgnd = 0 ENDIF ENDIF ENDIF ; END fix_pattern ENDIF ; END fix_something ; The vector back[1:2] is a gradient vector in the xy-plane of the ; skymap indicating the 'uphill' direction of the background plane. ; For the polar maps transform the gradient vector into a coordinate ; frame with the x-axis in the direction of increasing RA and the ; y-axis in the the direction of increasing declination (same as for ; the equatorial map). For the north polar map this amounts to a ; rotation of the coordinate frame over angle 90+atan(y/x), where (x,y) ; is the distance to the north pole in bins (units of array index). ; For the south polar map the y-axis needs to be flipped in addition. back_slope = back[1:2] CASE maptype EQ 0 OF 0: BEGIN tmp = xystar[*,ibright] ; Distance to pole ;tmp = xystar[*,ibright]+ddxystar[*,0] tmp -= origin ; Distance to pole tmp = atan(tmp[1],tmp[0]) ; atan(y/x) rpole = sqrt(total(tmp*tmp)) ; Polar distance ; Set up rotation matrix for rotation over 90+atan(y/x) tmp = -[[sin(tmp),maptype*cos(tmp)],[-cos(tmp),maptype*sin(tmp)]] ; Rotate back_slope = tmp#back_slope ; back_slope still is in adus/bin. Multiply by the map scaling factor ; in bins/angle. In the x-direction (RA) this is rpole bins/radian or ; rpole*rpm bins/angle. In the y-direction (declination) the scaling is ; 1.0/scale bins/angle back_slope *= [rpole*rpm,1.0/scale] END 1: back_slope /= scale ; From adus/bin to adus/degree ENDCASE ; The width of the PSF in bins is calculated as a vector with ; two components: the weighted mean of the x,y distance to the ; centroid xystar. The weight factor is the star intensity. psfwidth = fltarr(2,n_crowd) ; PSF width bright_raw = fltarr( n_crowd) ; Average adus/sqdeg in PSF in raw star bright_mdl = fltarr( n_crowd) ; Average adus/sqdeg in PSF in model star FOR i=0,n_crowd-1 DO BEGIN ; Loop over stars fitted simultaneously ; Bins inside PSF in_fullpsf_tmp = where(one_star_model[*,i] GT 0.0) ;in_fullpsf_tmp = where(in_fitpsf) ;?? IF in_fullpsf_tmp[0] NE -1 THEN BEGIN ; Subtract the background from the skymap intensities tmp = (original_star-bkgnd_model)[in_fullpsf_tmp] ; Absolute value of RA/dec in standard star ; Remember that the standard star has its centroid at RA=dec=0, ; with the symmetry axis in the vertical (declination; ; 2nd dimension) direction rr0 = abs(std_rr[*,in_fullpsf_tmp,i]) ; PSF width measured in angular units perpendicular and along the ; the symmetry axis of the PSF: mean of RA,dec weighted with the ; brightness of the original star (after background subtraction) CASE n_elements(tmp) EQ 1 OF 0: psfwidth[*,i] = total(rr0*(replicate(1,2)#(tmp/total(tmp,/nan))),2,/nan) 1: psfwidth[*,i] = rr0 ENDCASE ; Mean intensity over all bins inside PSF for original star, and ; the same for the modeled representation of the star. bright_raw[i] = (scale*scale)*total(tmp,/nan)/stdadus[i] tmp = one_star_model[in_fullpsf_tmp,i] bright_mdl[i] = (scale*scale)*total(tmp,/nan)/stdadus[i] ENDIF ENDFOR ; Difference between skymap and lsq fit (bkgnd+star) ; This is what should be close to zero. residual_sky = original_star-(bkgnd_model+star_model) ; Standard deviations of residual background and residual star stdev_psf = stddev(residual_sky[in_fitpsf],/nan) IF fit_bkgnd THEN stdev_bkgnd = stddev(residual_sky[in_bkgnd],/nan) ELSE stdev_bkgnd = 0.0 ; Results of fit IF silent LE 0 THEN $ print, format='(I4,2I5,2F9.3,2F6.2,4F10.4)' , $ ibright, n_fitpsf, n_bkgnd , $ back[0], bright[0] , $ ddxystar[*,0] , $ mean(residual_sky[in_fitpsf],/nan) , $ stdev_psf , $ cvfit , $ correlate(residual_sky[in_fitpsf],original_star[in_fitpsf]) CASE kill_one OF 0: BEGIN ; Store results of fit for all stars fit simultaneously IF has_sigma THEN BEGIN j = where(star_model GT 0.0) ; Inside PSF (same as in_fitpsf?) i = where(residual_sky[j] GT sigma*stdev_psf) IF i[0] NE -1 THEN BEGIN ; Outliers inside bright star PSF i = j[i] star_model[i] += residual_sky[i] ENDIF ENDIF final_out = star_model flat_out = star_model ; Combine all the core psfs from all stars in the crowd. tmp = in_corepsf[*,0] FOR i=1,n_crowd-1 DO tmp OR= in_corepsf[*,i] tmp = where(tmp) IF tmp[0] NE -1 THEN flat_out[tmp] = (original_flat-bkgnd_model)[tmp] ; The > 0 avoids subtraction of stars fitted with a negative ; brightness. This should not happen, but it does. final_out >= 0 ; Safety belt flat_out >= 0 IF star_fit[ibright].remove THEN BEGIN final_sky[ipix,jbeg:jend] -= final_out flat_sky [ipix,jbeg:jend] -= flat_out star_sky [ipix,jbeg:jend] += flat_out ;final_out ENDIF star_fit[crowd].total_area = stdarea star_fit[crowd].total_adus = stdadus star_fit[crowd].bright_raw = bright_raw star_fit[crowd].bright_mdl = bright_mdl star_fit[crowd].bright = bright star_fit[crowd].back_const = back[0] star_fit[crowd].back_slope = back_slope #one_crowd star_fit[crowd].n_bkgnd = n_bkgnd star_fit[crowd].n_fitpsf = n_fitpsf star_fit[crowd].stdev_psf = stdev_psf star_fit[crowd].stdev_bkgnd = stdev_bkgnd star_fit[crowd].cvfit = cvfit star_fit[crowd].dradec = ddradec star_fit[crowd].dpsfangle = dpsfangle *one_crowd star_fit[crowd].dpsfstretch = dpsfstretch#one_crowd star_fit[crowd].psfwidth = psfwidth IF show NE '' AND show NE '*' THEN $ IF (where(strpos(names[crowd],show) NE -1))[0] NE -1 THEN $ show = strtrim(names[ibright]) END 1: BEGIN IF has_sigma THEN BEGIN j = where(one_star_model[*,0] GT 0.0); Inside bright star PSF i = where(residual_sky[j] GT sigma*stdev_psf) IF i[0] NE -1 THEN BEGIN ; Outliers inside bright star PSF i = j[i] one_star_model[i,0] += residual_sky[i] ENDIF ENDIF ; Only subtract the bright star itself. ; Fainter stars will be fit and subtracted later. final_out = one_star_model[*,0] flat_out = final_out ; In flat_out replace the core PSF area by the background model. tmp = where(in_corepsf[*,0]) IF tmp[0] NE -1 THEN flat_out[tmp] = (original_flat-bkgnd_model)[tmp] ; The > 0 avoids subtraction of stars fitted with a negative ; brightness. This should not happen, but it does. final_out >= 0 ; Safety belt flat_out >= 0 IF star_fit[ibright].remove THEN BEGIN final_sky[ipix,jbeg:jend] -= final_out flat_sky [ipix,jbeg:jend] -= flat_out star_sky [ipix,jbeg:jend] += flat_out ;final_out ENDIF ; Store results of fit for brightest star in star_fit structure star_fit[ibright].total_area = stdarea [0] star_fit[ibright].total_adus = stdadus [0] star_fit[ibright].bright_raw = bright_raw[0] star_fit[ibright].bright_mdl = bright_mdl[0] star_fit[ibright].bright = bright[0] star_fit[ibright].back_const = back [0] star_fit[ibright].back_slope = back_slope star_fit[ibright].n_bkgnd = n_bkgnd star_fit[ibright].n_fitpsf = n_fitpsf star_fit[ibright].stdev_psf = stdev_psf star_fit[ibright].stdev_bkgnd= stdev_bkgnd star_fit[ibright].cvfit = cvfit star_fit[ibright].dradec = ddradec[*,0] star_fit[ibright].dpsfangle = dpsfangle star_fit[ibright].dpsfstretch= dpsfstretch star_fit[ibright].psfwidth = psfwidth[*,0] END ENDCASE IF show EQ '*' OR show EQ strtrim(names[ibright]) THEN BEGIN ; Original star view_star *= 0.0 view_star[ibeg-ibeg0:iend-ibeg0,jbeg-jbeg0:jend-jbeg0] = original_star IF IsType(min_show,/defined) THEN $ view_star = view_star > min_show[0] IF IsType(max_show,/defined) THEN $ view_star = view_star < max_show[0] view, in=MagnifyArray(view_star,5), /stretch, /bl, title="Original star" area_show = view_star ; Same area after subtraction, i.e. ; the residual after stars have been taken out view_star *= 0.0 view_star[ibeg-ibeg0:iend-ibeg0,jbeg-jbeg0:jend-jbeg0] = final_sky[ipix,jbeg:jend] IF IsType(min_show,/defined) THEN $ view_star = view_star > min_show[n_elements(min_show)-1] IF IsType(max_show,/defined) THEN $ view_star = view_star < max_show[n_elements(max_show)-1] boost, area_show, view_star, push=2 view, in=MagnifyArray(view_star,5), /stretch, /br, title="Residual star" ; Same area after subtraction, i.e. ; the residual after stars have been taken out with star replaced by flat background view_star *= 0.0 view_star[ibeg-ibeg0:iend-ibeg0,jbeg-jbeg0:jend-jbeg0] = flat_sky[ipix,jbeg:jend] IF IsType(min_show,/defined) THEN $ view_star = view_star > min_show[n_elements(min_show)-1] IF IsType(max_show,/defined) THEN $ view_star = view_star < max_show[n_elements(max_show)-1] boost, area_show, view_star, push=2 view, in=MagnifyArray(view_star,5), /stretch, /tc, title="Flat background" ; Same area indicating which bins were used for PSF and background fitting tmp = 0.0*original_star tmp[in_fitpsf] = 1.0 IF in_bkgnd[0] NE -1 THEN tmp[in_bkgnd] = -1.0 CASE n_crowd NE 1 OF 0: i = in_corepsf 1: i = total(in_corepsf,2) GT 0 ENDCASE i = where(i) IF i[0] NE -1 THEN tmp[i] += 1.0 view_star = 0.0*view_star view_star[ibeg-ibeg0:iend-ibeg0,jbeg-jbeg0:jend-jbeg0] = tmp boost, area_show, view_star, push=2 view, in=MagnifyArray(view_star,5), /stretch, /tl, title="PSF Mask" ; Scatter plot of star_model vs original star (with background subtracted) tmp = (original_star-bkgnd_model)[in_fitpsf] frange = 1.05* $ [ min( [star_model[in_fitpsf], tmp] ) < 0, $ max( [star_model[in_fitpsf], tmp] ) ] twin, /show, xsize=300, ysize=300, /tr, title="Model vs Data" plot, star_model[in_fitpsf]/frange[1] , $ tmp/frange[1] , $ psym = 4 , $ xstyle = 1 , $ ystyle = 1 , $ xrange = frange/frange[1] , $ yrange = frange/frange[1] , $ xtitle = 'standard' , $ ytitle = 'map' oplot, frange, frange IF NOT fullmap THEN $ view, in=(final_sky > (-100)) < 100,/stretch, /bc ; Histogram of residuals ;IF show NE '*' THEN qbar, /histo, residual_sky[in_fitpsf] ENDIF CASE kill_one OF 0: star_fit[crowd ].done = 1 1: star_fit[ibright].done = 1 ENDCASE ENDIF ELSE BEGIN ; Not enough points to fit IF show EQ '*' OR show EQ strtrim(names[ibright]) THEN BEGIN view_star = 0.0*view_star view_star[ibeg-ibeg0:iend-ibeg0,jbeg-jbeg0:jend-jbeg0] = original_star view, in=MagnifyArray(view_star,5), /stretch, /bl ENDIF IF silent LE 1 THEN message, /info, $ names[ibright]+', only' + $ strcompress(n_fitpsf)+ ' PSF and' + $ strcompress(n_bkgnd )+' background points used for fit' CASE kill_one OF 0: BEGIN star_fit[crowd].n_bkgnd = n_bkgnd star_fit[crowd].n_fitpsf = n_fitpsf star_fit[crowd ].done = 2 END 1: BEGIN star_fit[ibright].n_bkgnd = n_bkgnd star_fit[ibright].n_fitpsf = n_fitpsf star_fit[ibright].done = 2 END ENDCASE ENDELSE ; Quit if the requested star has just been displayed IF show EQ strtrim(names[ibright]) THEN break ENDFOR CASE count NE 0 OF 0: BEGIN IF silent LE 1 THEN message, /info, error star_index = -1 END 1: star_index = valid ENDCASE IF count NE 0 THEN IF (where(star_fit.done EQ 0))[0] NE -1 THEN message, 'oops' RETURN & END