;+ ; NAME: ; smei_star_stdmaps ; PURPOSE: ; Get brightness in standard star at indicated locations ; CATEGORY: ; camera/idl/star ; CALLING SEQUENCE: FUNCTION smei_star_stdmaps, camera, $ maptype , $ rr , $ origin , $ scale , $ radec , $ psfangle , $ psfstretch , $ degrees = degrees , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ in_fullpsf = in_fullpsf , $ in_fitpsf = in_fitpsf , $ in_corepsf = in_corepsf , $ in_bkgnd = in_bkgnd , $ std_rr = std_rr , $ use_filled_psf=use_filled_psf,$ stars_mask = stars_mask , $ corepsf_mask = corepsf_mask , $ stdarea = stdarea , $ stdadus = stdadus ; INPUTS: ; camera scalar; type: integer ; camera id (1,2,3) ; maptype scalar; type: integer ; -1: map of south pole ; 0: equatorial map ; +1: map of north pole ; rr array[2,n]; type: integer ; position vectors for all bins in the box around ; the star as column and row indices into the original skymap ; for all bins in box around star ; origin array[2]; type: float ; array indices of RA=0, dec=0 (equatorial map), ; position of north pole, dec=90 (north polar map) or ; position of south pole, dec=-90 (south polar map) ; scale scalar; type: float ; map resolution (radians/degrees per bin) ; radec array[2,m]; type: float ; RA and dec of the star centroids for m stars ; psfangle array[m]; type: float ; orientations of PSF relative to equatorial north ; for m stars ; psfstretch array[2,m]; type: float ; 'stretch' factors along long and narrow dimension of the ; SMEI field of view for m stars. Larger stretch factors ; will shrink the PSF size (similar to the geometric effect ; of the shrinking aperture size when viewed from a direction ; away from the optical axis (the fov direction cosine angle)). ; OPTIONAL INPUTS: ; /degrees if set, all angles are in degrees (default: radians) ; ; stars_mask=stars_mask ; array[n]; type: float ; box around star with results for previously- ; subtracted stars. ; This is a box extracted from the "stars_sky" ; array in smei_star_fit: it is zero in skybins ; where nothing has been subtracted yet, and ; positive-definite when subtraction has taken ; place (the content is the intensity that was ; subtracted from the original skymap). ; If present this affects 'in_fit_psf' and 'in_bkgnd': ; in_bkgnd: skybins that are part of an already-subtracted ; star are excluded. ; in_fitpsf: skybins that that are part of an already-subtracted ; star and are in the core PSF of one of the m stars ; are excluded. ; ; corepsf_mask=corepsf_mask ; array[n]; type: byte ; 0 where bins are inside the core PSF of bright stars. ; This mask affects 'in_fit_psf': ; skybins inside the core PSF of bright stars, but outside ; the core PSF of the m stars are excluded ; ; Passed to href=smei_star_standard=: ; ; bkgnd_origin=bkgnd_origin ; array[2]; type: float ; bkgnd_radius=bkgnd_radius ; scalar or array[2]; type: float ; wing_origin=wing_origin ; array[2]; type: float ; wing_radius=wing_radius ; scalar or array[2]; type: float ; wing_delta=wing_delta ; scalar or array[2]; type float ; ; /use_filled_psf if set, use the alternative PSF map ; OUTPUTS: ; Result array[n,m]; type: float ; 'standard star' brightnesses at the locations rr for m stars. ; The 'standard star' gives the brightness of a star located ; at RA=0, dec=0 with PSF angle of zero. ; The sky locations rr are transformed to locations in the ; 'standard star' reference frame; 'standard star' brightnesses ; are obtained by interpolation. ; in_fullpsf=in_fullpsf ; array[n]; type: byte ; 1B identifies locations inside combined PSF of ; m stars (i.e. the union of the non-zero location ; for the individual stars). ; This includes locations inside the "core" PSF, and ; (if wing_radius is used) additional locations where ; the core PSF is extended out to the specified ; wing_radius. ; in_fitpsf=in_fitpsf ; array[n]; type: byte ; subset of in_fullpsf ; excludes bins set to 1B in in_fullpsf based ; on info provided in keyword stars_mask ; and corepsf_mask ; in_corepsf=in_corepsf[n,m] ; array[n]; type: byte ; 1B identifies locations inside the "core" PSF ; (stored in the standard star map) for each of the ; m fitted stars individually. ; If wing_radius is set this is NOT the same as in_fitpsf ; (in_corepsf is a subset of in_fitpsf, identifying the ; PSF before it was extended out to wing_radius) ; in_bkgnd=in_bkgnd ; array[n]; type: byte ; 1B identifies locations in background area ; Is affected by stars_mask, if specified ; ; std_rr=std_rr array[2,n,m]; type: float ; equivalent RA,dec in the 'standard star' reference frame ; (centered on RA=0, dec=0) with abs(RA) < 180 degrees ; stdarea=stdarea array[m]; type: float ; total angular size of the PSF for all m stars. ; Calculated in the 'standard star' reference frame, and ; includes all bins with non-zero brightness (i.e. if ; wing_radius is used this includes the bins outside the ; core PSF filled with an exponential dropoff out to ; wing_radius. (!!! Not sure this is the right way to ; do it: maybe should use core PSF only) ; stdadus=stdadus array[m]; type: float ; Total adus in the 'standard star' in the same set of ; bins used to calculate 'stdarea'. ; The sum is defined as the sum of the ADUs in the ; standard star skybins, times the binarea (0.025^2 deg). ; (this makes totaladus independent of bin size, making ; it easier to compare with the same quantity calculated ; in the fullsky maps (with 0.1 deg resolution). ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; ToRadians, ToDegrees, EulerRotate, AngleRange, bilinear ; smei_star_standard, IsType ; PROCEDURE: ; See Section 6 in href=http://supercat.ucsd.edu/reports/smei.pdf= ; MODIFICATION HISTORY: ; JUN-2006, Paul Hick (UCSD/CASS) ; JUL-2006, Jordan Vaughan (jtvaugha@ucsd.edu) ; Added costhetax argument ; JUL-2006, Paul Hick (UCSD/CASS) ; Modified the linear interpolation on the standard map. ; The 'missing data' value now is -1.0 instead of 0.0. ; The value 0.0 is now used to indicate bins used to sample ; the background. ; JUL-2006, Paul Hick (UCSD/CASS) ; Added /use_filled_psf keyword ; FEB-2007, Jordan Vaughan, Paul Hick (UCSD/CASS) ; Added stars_mask keyword to be able to exclude skybins ; in already-subtracted stars from the background. ; JUL-2007, Paul Hick (UCSD/CASS) ; Added keyword in_corepsf ; NOV-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added more documentation. ;- rpm = ToRadians(degrees=degrees) dpm = ToDegrees(degrees=degrees) ; Use inverse of map projection to convert from array indices to RA, dec. ; For the equatorial map the map projection is a linear transformation ; for both angles (equi-rectangular projection). For the south- and ; northpole maps the map projection is a polar transformation. np = n_elements(rr)/2 ; Number of x,y pairs one = replicate(1.0,np) pp = rr-origin#one ; Array indices relative to origin CASE maptype EQ 0 OF 0: pp = transpose([[atan(reform(pp[1,*]),reform(pp[0,*]))/rpm],[maptype*(90.0/dpm-sqrt(total(pp*pp,1))*scale)]]) 1: pp *= scale ENDCASE ; pp now contains RA and dec ns = n_elements(psfangle) ; Number of stars fit simultaneously stdfit = fltarr(np,ns) stdval = fltarr(np,ns) std_rr = fltarr(2,np,ns) stdarea = fltarr(ns) stdadus = fltarr(ns) in_corepsf = bytarr(np,ns) nomansland_id = smei_star_standard(camera , $ use_filled_psf = use_filled_psf , $ wing_radius = wing_radius , $ /get_nomansland_id) background_id = smei_star_standard(camera , $ use_filled_psf = use_filled_psf , $ wing_radius = wing_radius , $ /get_background_id) wayoutside_id = smei_star_standard(camera , $ use_filled_psf = use_filled_psf , $ wing_radius = wing_radius , $ /get_wayoutside_id) nogrid_id = min([nomansland_id,background_id,wayoutside_id])-1.0 FOR i=0,ns-1 DO BEGIN ; Loop over psf angles (will be more than one ; ... when fitting multiple stars at once). qq = pp ; RA and dec in skymap ; Transform coordinates into equivalent RA,dec in standard star. ; EulerRotate rotates the PSF with the tail pointing north, ; and translates from the real star location to the standard ; star located at RA=0,dec=0 tmp = radec[1,i]*rpm sdec = sin(tmp) cdec = cos(tmp) tmp = psfangle[i]*rpm spsf = sin(tmp) cpsf = cos(tmp) tmp = [ -(atan(sdec*cpsf,spsf)+!pi/2.0-radec[0,i]*rpm), $ acos(cdec*cpsf), atan(sdec,cdec*spsf)+!pi/2.0 ]/rpm qq = EulerRotate(tmp,qq,degrees=degrees) qq[0,*] = AngleRange(qq[0,*],/pi,degrees=degrees) ; Stretch the RA in the standard star by psfstretch. ; This increases the x-coordinates so the standard star will ; be sampled further away from the centroid. This takes care ; of the shrinking of the PSF with psfstretch. ; Note that the equivalent brightness interpolated on the ; standard map below needs the same correction to conserve ; the brightness inside the PSF. qq[0,*] /= psfstretch[0] qq[1,*] /= psfstretch[1] ; qq[0,*] is RA in standard star ; qq[1,*] is dec in standard star std_rr[*,*,i] = qq ; Find interpolated intensities in the standard star. ; The map returned by smei_star_standard is a 200x200 ; array with several zones: ; 1. The PSF proper (or 'core PSF') with values GT zero. ; This is Aarons original standard star. ; If the bkgnd_* keywords are used: ; 2. An area between PSF and the inner radius specified in ; bkgnd_radius with values EQ -1.0 (may not be present) ; 3. A surrounding area with values EQ -2.0 between inner ; and outer radius. These are used in the fit and are ; supposed to sample the background. ; 4. A third zone with values EQ -3.0 outside the ; background or wing radius area. These just fill out the ; square standard map and are not used in the fit. ; If the wing_* keywords are used: ; 2. A surrounding area with values EQ -2.0 between 'core PSF' ; and wing radius. The combined bins in 1. and 2. are used ; to fit star(s) and background simultaneously. ; 3. A third zone with values EQ -3.0 outside the ; background or wing radius area. These just fill out the ; square standard map and are not used in the fit. stdmap = smei_star_standard(camera , $ bkgnd_origin = bkgnd_origin , $ bkgnd_radius = bkgnd_radius , $ wing_origin = wing_origin , $ wing_radius = wing_radius , $ wing_delta = wing_delta , $ psfstretch = psfstretch[*,i],$ use_filled_psf = use_filled_psf, $ in_psf = in_stdpsf , $ degrees = degrees , $ scale = std_scale , $ origin = std_origin ) stddim = size(stdmap,/dim) ; wing_radius SET: ; stdmap > 0 includes all bins inside the elllipse defined by wing_radius ; wing_radius NOT set: ; stdmap > 0 includes all bins inside the core PSF inside = where(stdmap GT 0.0,tmp) stdarea[i] = (std_scale*std_scale)*tmp ; Integrated area stdadus[i] = (std_scale*std_scale)*total(stdmap[inside]); Integrated adus ; Convert to pixel coordinates in the standard star qq = std_origin#one+qq/std_scale xx = reform(qq[0,*]) yy = reform(qq[1,*]) ; Some locations will fall outside the standard star map inside = where( 0 LE xx AND xx LE stddim[0]-1 AND $ 0 LE yy AND yy LE stddim[1]-1, complement=outside ) IF inside[0] NE -1 THEN BEGIN ; Locations inside standard map ; The locations that still fall inside the 200x200 ; standard map are interpolated. xx = xx[inside] yy = yy[inside] ; Pick up the value in the standard star for the bins nearest ; the equivalent locations qq. These will be used below to ; select the actual bins to be used in the lsq fitting ; procedure (see smei_star_lsqfit) by testing against ; nomansland_id, background_id and wayoutside_id. ; Why not use linear interpolation here? stdval[inside,i] = stdmap[round(xx),round(yy)] ; in_corepsf identifies bins inside the original core PSF ; (excluding the wing added by keyword wing_radius) in_corepsf[inside,i] = in_stdpsf[round(xx),round(yy)] ; Find the value at the equivalent locations qq by linear ; interpolation on the standard star. Note that we ; interpolate on stdmap > 0, so all values in stdfit will ; be >= 0 with zero values occurring outside the PSF. ; The correction psfstretch makes sure the integrated ; brightness is conserved. ; The transpose of xx,yy is really necessary! Read up on the ; IDL bilinear procedure if you're not buying it! stdfit[inside,i] = bilinear(stdmap > 0.0,transpose(xx), $ transpose(yy),missing=nogrid_id)/(psfstretch[0]*psfstretch[1]) ; The missing keyword to bilinear should be redundant since ; the standard map 'stdmap' does not have any NaNs. IF (where(stdfit[inside,i] LT 0.0))[0] NE -1 THEN $ message, 'oops' ENDIF IF outside[0] NE -1 THEN BEGIN ; Locations outside map: way outside stdval [outside,i] = wayoutside_id stdfit [outside,i] = 0.0 in_corepsf[outside,i] = 0 ENDIF ENDFOR ; Use array stdval to select the bins to be used to fit the star ; and to fit the background CASE ns EQ 1 OF 0: BEGIN ; Multiple stars fitted simultaneously in_fullpsf = total(stdval GT 0.0,2) GT 0 ; Selects bins inside (wing-extended) PSF for at least one star in_bkgnd = total(stdval EQ background_id,2) GT 0 AND $; Selects bins in background for at least one star total(stdval EQ nomansland_id,2) EQ 0 ; Selects bins not between PSF and background END 1: BEGIN ; Just one star fitted in_fullpsf = stdval GT 0.0 in_bkgnd = stdval EQ background_id END ENDCASE tmp = in_corepsf[*,0] FOR i=1,ns-1 DO tmp OR= in_corepsf[*,i] ; Union of all core PSFs in_fitpsf = in_fullpsf IF IsType(stars_mask,/defined) THEN BEGIN in_fitpsf AND= stars_mask EQ 0 OR 1-tmp ; Avoid bins in core PSF that have been subtracted already in_bkgnd AND= stars_mask EQ 0 ; Avoid bins that have been subtracted already ENDIF ; PPH (2011-11-04) ; I think, the inclusion of 'tmp' below is unnecessary (always true): ; the union of all core PSFs should be a subset of in_fitpsf at this point. IF IsType(corepsf_mask,/defined) THEN $ in_fitpsf AND= corepsf_mask OR tmp ; Avoid bins outside core PSF of m stars, that are ; inside the core PSF of brighter stars RETURN, stdfit & END