Source code for lsst.sims.maf.stackers.NEODistStacker

import numpy as np
from .baseStacker import BaseStacker
from .generalStackers import FiveSigmaStacker

__all__ = ['NEODistStacker']


[docs]class NEODistStacker(BaseStacker): """ For each observation, find the max distance to a ~144 km NEO, also stack on the x,y position of the object. """ colsAdded = ['MaxGeoDist', 'NEOHelioX', 'NEOHelioY'] def __init__(self, stepsize=.001, maxDist=3., minDist=.3, H=22, elongCol='solarElong', filterCol='filter', sunAzCol='sunAz', azCol='azimuth', m5Col='fiveSigmaDepth'): """ stepsize: The stepsize to use when solving (in AU) maxDist: How far out to try and measure (in AU) H: Asteroid magnitude Adds columns: MaxGeoDist: Geocentric distance to the NEO NEOHelioX: Heliocentric X (with Earth at x,y,z (0,1,0)) NEOHelioY: Heliocentric Y (with Earth at (0,1,0)) Note that both opsim v3 and v4 report solarElongation in degrees. """ self.units = ['AU', 'AU', 'AU'] # Also grab things needed for the HA stacker self.colsReq = [elongCol, filterCol, sunAzCol, azCol, m5Col] self.sunAzCol = sunAzCol self.elongCol = elongCol self.filterCol = filterCol self.azCol = azCol self.m5Col = m5Col self.H = H # Magic numbers (Ivezic '15, private comm.)that convert an asteroid # V-band magnitude to LSST filters: # V_5 = m_5 + (adjust value) self.limitingAdjust = {'u': -2.1, 'g': -0.5, 'r': 0.2, 'i': 0.4, 'z': 0.6, 'y': 0.6} self.deltas = np.arange(minDist, maxDist+stepsize, stepsize) self.G = 0.15 # Magic numbers from http://adsabs.harvard.edu/abs/2002AJ....124.1776J self.a1 = 3.33 self.b1 = 0.63 self.a2 = 1.87 self.b2 = 1.22 def _run(self, simData, cols_present=False): if cols_present: # This is a pretty rare stacker. Assume we need to rerun pass elongRad = np.radians(simData[self.elongCol]) v5 = np.zeros(simData.size, dtype=float) + simData[self.m5Col] for filterName in self.limitingAdjust: fmatch = np.where(simData[self.filterCol] == filterName) v5[fmatch] += self.limitingAdjust[filterName] for i, elong in enumerate(elongRad): # Law of cosines: # Heliocentric Radius of the object R = np.sqrt(1.+self.deltas**2-2.*self.deltas*np.cos(elong)) # Angle between sun and earth as seen by NEO alphas = np.arccos((1.-R**2-self.deltas**2)/(-2.*self.deltas*R)) ta2 = np.tan(alphas/2.) phi1 = np.exp(-self.a1*ta2**self.b1) phi2 = np.exp(-self.a2*ta2**self.b2) alpha_term = 2.5*np.log10((1. - self.G)*phi1+self.G*phi2) appmag = self.H+5.*np.log10(R*self.deltas)-alpha_term # There can be some local minima/maxima when solving, so # need to find the *1st* spot where it is too faint, not the # last spot it is bright enough. tooFaint = np.where(appmag > v5[i]) # Check that there is a minimum if np.size(tooFaint[0]) == 0: simData['MaxGeoDist'][i] = 0 else: simData['MaxGeoDist'][i] = np.min(self.deltas[tooFaint]) # Make coords in heliocentric interior = np.where(elongRad <= np.pi/2.) outer = np.where(elongRad > np.pi/2.) simData['NEOHelioX'][interior] = simData['MaxGeoDist'][interior]*np.sin(elongRad[interior]) simData['NEOHelioY'][interior] = -simData['MaxGeoDist'][interior]*np.cos(elongRad[interior]) + 1. simData['NEOHelioX'][outer] = simData['MaxGeoDist'][outer]*np.sin(np.pi-elongRad[outer]) simData['NEOHelioY'][outer] = simData['MaxGeoDist'][outer]*np.cos(np.pi-elongRad[outer]) + 1. # Flip the X coord if sun az is negative? if simData[self.azCol].min() < - np.pi/2.0: halfval = 180. else: halfval = np.pi flip = np.where(((simData[self.sunAzCol] > halfval) & (simData[self.azCol] > halfval)) | ((simData[self.sunAzCol] < halfval) & (simData[self.azCol] > halfval))) simData['NEOHelioX'][flip] *= -1. return simData