Adding lightcurve effects to LSST Solar System Simulations ========================================================== By Jamie Robinson LSST SSSC Sprint June 2022 This notebook demonstrates the work I did at the Sprint. I wanted to add changes in asteroid brightness due to rotation and also the longer term changes of aspect angle that leads to shifts in brightness (on top of the already modelled distance and phase effects). The asteroid is modelled as a triaxial ellipsoid rotating around its shortest axis - axes a > b > c - spin pole orientation (ecliptic longitude and latitude) - rotation period and simple sinusoidal lightcurve This shape model is applied to some asteroid from the simulated database. The database contains only sparse measurements as expected from the LSST cadence. Therefore we take the orbital elements and propagate using `sbpy open_orb `__ to densely sample the asteroid positional coordinates. These positions and and shape parameters are fed into the equations from `Fernández-Valenzuela 2022 `__ and so a shift from the expected brightness due to lightcurve effects can be calculated for any observation made with LSST. Based on the demo notebook by Juric et al., available from https://github.com/lsst-sssc/lsst-simulation .. code:: ipython3 %matplotlib notebook .. code:: ipython3 %matplotlib inline .. code:: ipython3 import psycopg2 as pg import pandas as pd import numpy as np import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from astropy import units as u from astropy.coordinates import SkyCoord, GCRS from sbpy.data import Orbit from astropy.time import Time from astropy.coordinates import solar_system_ephemeris from sbpy.photometry import HG from astropy.table import QTable .. code:: ipython3 # connect to database pwd = open("/home/shared/sssc-db-pass.txt").read() con = pg.connect(database="lsst_solsys", user="sssc", password=pwd, host="epyc.astro.washington.edu", port="5432") Select an asteroid from the database ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 # select some asteroids with a reasonable number of observations df_ssobj = pd.read_sql("SELECT * FROM ssobjects where numobs>100 LIMIT 5", con) .. code:: ipython3 df_ssobj .. raw:: html
ssobjectid discoverysubmissiondate firstobservationdate arc numobs moid moidtrueanomaly moideclipticlongitude moiddeltav uh ... yg12 yherr yg12err yh_yg12_cov ychi2 yndata maxextendedness minextendedness medianextendedness flags
0 3369984299447 60203.374981 60196.374981 2582.7770 107 0.0 0.0 0.0 0.0 NaN ... 0.345189 0.087608 0.133565 0.010005 0.571543 7 0.0 0.0 0.0 0
1 5992863104062 60040.337008 60033.337008 3241.8570 145 0.0 0.0 0.0 0.0 NaN ... NaN NaN NaN NaN NaN 0 0.0 0.0 0.0 0
2 10660550296725 59873.208132 59866.208132 3442.9424 116 0.0 0.0 0.0 0.0 NaN ... NaN NaN NaN NaN NaN 0 0.0 0.0 0.0 0
3 14169511631100 59963.234375 59956.234375 3333.8674 289 0.0 0.0 0.0 0.0 NaN ... 0.105138 0.038684 0.090257 0.002581 0.736188 13 0.0 0.0 0.0 0
4 14369040295367 59865.265384 59858.265384 2975.7960 308 0.0 0.0 0.0 0.0 20.448154 ... 0.185439 0.020225 0.025261 0.000455 0.792384 47 0.0 0.0 0.0 0

5 rows × 55 columns

.. code:: ipython3 # Choose some object by its identifier "ssobjectid" ssobjectid = 14169511631100 .. code:: ipython3 # Query to find the mpc designation identifier that corresponds to ssobjectid qry = """ SELECT ssobjectid, mpcorb.mpcdesignation FROM ssobjects JOIN mpcorb USING (ssobjectid) WHERE ssobjectid='{}' """.format(ssobjectid) df = pd.read_sql(qry, con) .. code:: ipython3 mpc_des = df.iloc[0]["mpcdesignation"] mpc_des .. parsed-literal:: 'S100doiQa ' Retrieve observations ~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 sql = """ SELECT mpcdesignation, ssObjects.ssObjectId, mag, magSigma, filter, midPointTai as mjd, ra, decl, phaseAngle, topocentricDist, heliocentricDist, eclipticlambda, eclipticbeta, heliocentricx, heliocentricy, heliocentricz, topocentricx, topocentricy, topocentricz, heliocentricvx, heliocentricvy, heliocentricvz FROM mpcorb JOIN ssObjects USING (ssobjectid) JOIN diaSources USING (ssobjectid) JOIN ssSources USING (diaSourceid) WHERE mpcdesignation = '{}' """.format(mpc_des) df = pd.read_sql(sql, con) # Distance correction to get reduced magnitude df["cmag"] = df["mag"] - 5*np.log10(df["topocentricdist"]*df["heliocentricdist"]) .. code:: ipython3 # put obs in time order df=df.sort_values("mjd") .. code:: ipython3 df .. raw:: html
mpcdesignation ssobjectid mag magsigma filter mjd ra decl phaseangle topocentricdist ... heliocentricx heliocentricy heliocentricz topocentricx topocentricy topocentricz heliocentricvx heliocentricvy heliocentricvz cmag
221 S100doiQa 14169511631100 21.747236 0.117441 i 59956.234375 135.646361 11.233583 8.342505 1.823498 ... -1.638028 2.090419 0.719361 -1.278892 1.250356 0.355236 -0.008227 -0.005801 -0.002566 18.244902
227 S100doiQa 14169511631100 21.791485 0.057519 i 59957.229467 135.464059 11.268270 7.959714 1.818000 ... -1.646204 2.084631 0.716803 -1.270909 1.250487 0.355243 -0.008204 -0.005831 -0.002576 18.295850
228 S100doiQa 14169511631100 21.711332 0.073672 z 59957.245586 135.460958 11.268851 7.953384 1.817913 ... -1.646336 2.084537 0.716761 -1.270777 1.250493 0.355243 -0.008204 -0.005831 -0.002576 18.215804
222 S100doiQa 14169511631100 21.821917 0.049366 r 59958.189412 135.284311 11.303410 7.587669 1.812944 ... -1.654069 2.079020 0.714325 -1.263301 1.250825 0.355345 -0.008182 -0.005859 -0.002586 18.332467
231 S100doiQa 14169511631100 21.917894 0.046935 r 59958.189858 135.284238 11.303432 7.587493 1.812941 ... -1.654073 2.079018 0.714324 -1.263298 1.250825 0.355345 -0.008182 -0.005859 -0.002586 18.428446
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
282 S100doiQa 14169511631100 22.553791 0.043438 g 63268.088620 108.756170 17.076115 9.944378 1.848784 ... -1.282065 2.297494 0.813358 -0.568254 1.673430 0.542880 -0.009096 -0.004531 -0.002122 19.019655
277 S100doiQa 14169511631100 22.459470 0.034845 g 63268.177852 108.740690 17.079102 9.977975 1.849339 ... -1.282876 2.297089 0.813168 -0.567963 1.674059 0.543135 -0.009095 -0.004534 -0.002123 18.924694
278 S100doiQa 14169511631100 22.048769 0.039528 r 63268.193938 108.737908 17.079614 9.984020 1.849440 ... -1.283023 2.297016 0.813134 -0.567912 1.674174 0.543181 -0.009094 -0.004534 -0.002123 18.513877
276 S100doiQa 14169511631100 22.643286 0.037484 g 63276.115874 107.612004 17.338086 12.743919 1.906033 ... -1.354484 2.260079 0.795958 -0.550504 1.734147 0.568016 -0.008946 -0.004790 -0.002213 19.044015
220 S100doiQa 14169511631100 22.440725 0.132068 i 63290.101771 106.641982 17.732335 16.685580 2.036539 ... -1.477661 2.189978 0.763915 -0.555536 1.858531 0.620270 -0.008665 -0.005233 -0.002368 18.699518

289 rows × 23 columns

.. code:: ipython3 list(df) .. parsed-literal:: ['mpcdesignation', 'ssobjectid', 'mag', 'magsigma', 'filter', 'mjd', 'ra', 'decl', 'phaseangle', 'topocentricdist', 'heliocentricdist', 'eclipticlambda', 'eclipticbeta', 'heliocentricx', 'heliocentricy', 'heliocentricz', 'topocentricx', 'topocentricy', 'topocentricz', 'heliocentricvx', 'heliocentricvy', 'heliocentricvz', 'cmag'] .. code:: ipython3 # query to retrieve the HG12 phase curve parameters in all filters ssoId = int(df["ssobjectid"].iloc[0]) qry = """SELECT uH, uG12, uHErr, uG12Err, uChi2, gH, gG12, gHErr, gG12Err, gChi2, rH, rG12, rHErr, rG12Err, rChi2, iH, iG12, iHErr, iG12Err, iChi2, zH, zG12, zHErr, zG12Err, zChi2, yH, yG12, yHErr, yG12Err, yChi2 FROM ssObjects WHERE ssObjectId={}""".format(ssoId) df_hg = pd.read_sql(qry, con) df_hg .. raw:: html
uh ug12 uherr ug12err uchi2 gh gg12 gherr gg12err gchi2 ... zh zg12 zherr zg12err zchi2 yh yg12 yherr yg12err ychi2
0 None None None None None 18.362034 0.219264 0.011466 0.020526 1.318255 ... 17.785637 0.24849 0.028324 0.042346 1.183112 17.742504 0.105138 0.038684 0.090257 0.736188

1 rows × 30 columns

.. code:: ipython3 # choose some colours to represent the different filters col_dict = {"u":"C0","g":"C2","r":"C3","i":"C5","z":"C6","y":"C7"} .. code:: ipython3 # plot the phase curve in multiple filters x_plot = "phaseangle" y_plot = "cmag" c_plot = "filter" df_plot = df fig = plt.figure() gs = gridspec.GridSpec(1,1) ax1 = plt.subplot(gs[0,0]) alpha = np.linspace(np.amin(df_plot[x_plot]),np.amax(df_plot[x_plot])) * u.deg for f in col_dict.keys(): _df_plot = df_plot[df_plot["filter"]==f] x=ax1.scatter(_df_plot[x_plot],_df_plot[y_plot], c= col_dict[f], label = f) h = "{}h".format(f) g12 = "{}g12".format(f) try: model = HG(H = df_hg[h] * u.mag, G = df_hg[g12]) except: continue ax1.plot(alpha,model(alpha), c= col_dict[f]) ax1.set_xlabel(x_plot) ax1.set_ylabel(y_plot) ax1.invert_yaxis() ax1.legend() plt.show() .. image:: simulated_asteroid_lightcurves_files/simulated_asteroid_lightcurves_22_0.png .. code:: ipython3 # what is the typical time difference between observations? x_plot = "mjd" df_plot = df.sort_values(x_plot) # make sure df is in time order... bins = 100 fig = plt.figure() gs = gridspec.GridSpec(1,1) ax1 = plt.subplot(gs[0,0]) x=ax1.hist(np.diff(df_plot[x_plot]),bins = bins, log = True) ax1.set_xlabel(x_plot) ax1.set_ylabel("number") plt.show() .. image:: simulated_asteroid_lightcurves_files/simulated_asteroid_lightcurves_23_0.png Define some physical properties of our asteroid shape and rotation model ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 # calculate the spherical radius of the asteroid, given its absolute magnitude from the simulated database C = 1330 # km albedo = 0.15 # pick a typical asteroid albedo abs_mag = df_hg.iloc[0]["rh"] radius = 0.5*((C/np.sqrt(albedo)) * (10**(-0.2*abs_mag))) print("radius = {} km".format(radius)) .. parsed-literal:: radius = 0.4556819391654253 km .. code:: ipython3 print("circular area = {} km^2".format(np.pi*radius*radius)) .. parsed-literal:: circular area = 0.6523392413946846 km^2 .. code:: ipython3 volume = (4.0/3.0)*np.pi*(radius**3.0) print("volume = {} km^3".format(volume)) .. parsed-literal:: volume = 0.3963456140165764 km^3 .. code:: ipython3 def ellipsoid_axes(V,b_a,c_a): """ Find the axes of an ellipsoid given a volume V and axes ratios. b_a is the b/a ratio c_a is the c/a ratio """ a = (((3.0*V)/(4.0*np.pi)) * (1.0/b_a) *(1.0/c_a))**(1.0/3.0) b = b_a * a c = c_a * a return a,b,c .. code:: ipython3 # define the triaxial ellipsoid axes a,b,c = ellipsoid_axes(volume,0.75,0.5) print(a,b,c) print("volume = {} km^3".format((4.0/3.0)*np.pi*a*b*c)) .. parsed-literal:: 0.6319044200766154 0.4739283150574616 0.3159522100383077 volume = 0.3963456140165764 km^3 .. code:: ipython3 # set rotational properties period = 5.0/24 # rotation period (days) t0 = 0 # reference time for lightcurve phase lambda_pole = 0 # spin pole ecliptic longitude beta_pole = np.pi/2.0 # spin pole ecliptic latitude Resample the sparse observations by propagating orbit with sbpy & oorb ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 # get the orbital elements from the database sql = """ SELECT * FROM mpcorb WHERE mpcdesignation = '{}' """.format(mpc_des) df_orb = pd.read_sql(sql, con) df_orb.head() .. raw:: html
mpcdesignation mpcnumber ssobjectid mpch mpcg epoch tperi peri node incl ... arc arcstart arcend rms pertsshort pertslong computer flags fulldesignation lastincludedobservation
0 S100doiQa 0 14169511631100 18.03 0.15 54800.0 51911.34966 339.31535 222.04643 3.59192 ... None None None None None None None 0 2011 S100doiQ None

1 rows × 27 columns

.. code:: ipython3 list(df_orb) .. parsed-literal:: ['mpcdesignation', 'mpcnumber', 'ssobjectid', 'mpch', 'mpcg', 'epoch', 'tperi', 'peri', 'node', 'incl', 'e', 'n', 'q', 'uncertaintyparameter', 'reference', 'nobs', 'nopp', 'arc', 'arcstart', 'arcend', 'rms', 'pertsshort', 'pertslong', 'computer', 'flags', 'fulldesignation', 'lastincludedobservation'] .. code:: ipython3 # Calculate some extra orbital elements df_orb["a"] = df_orb["q"]/(1.0 - df_orb["e"]) df_orb["P"] = df_orb["a"]**(3.0/2.0) # orbital period in years df_orb["n"] = 360.0 / (df_orb["P"] * 365.25) # mean motion in deg/day df_orb["M"] = (df_orb["n"] * (df_orb["epoch"] - df_orb["tperi"])) % 360 # angles must be in correct range otherwise sbpy/pyoorb freak out df_orb[["a","P","n","M"]] .. raw:: html
a P n M
0 2.761968 4.590162 0.214726 260.267776
.. code:: ipython3 # rename columns for consistency df_orb = df_orb.rename(columns={"node":"Omega","peri":"w"}) .. code:: ipython3 df_orb[["a","e","incl","Omega","w","M"]] .. raw:: html
a e incl Omega w M
0 2.761968 0.01803 3.59192 222.04643 339.31535 260.267776
.. code:: ipython3 # create an sbpy oorb object from dataframe via QTable tab = QTable.from_pandas(df_orb[["a","e","incl","Omega","w","M"]], units={'a': u.au, 'incl': u.deg, 'Omega': u.deg, 'w': u.deg, 'M': u.deg}) orbit = Orbit.from_table(tab) .. code:: ipython3 # oorb requires certain extra fields orbit["epoch"] = Time(Time(df_orb["epoch"], format="mjd").jd, format = "jd") orbit["targetname"] = np.array(df_orb["ssobjectid"]).astype(str) orbit["H"] = df_orb["mpch"] * u.mag orbit["G"] = df_orb["mpcg"] * u.dimensionless_unscaled .. code:: ipython3 orbit .. parsed-literal:: a e incl Omega ... targetname H G AU deg deg ... mag float64 float64 float64 float64 ... str21 float64 float64 ------------------ ------- ------- --------- ... -------------- ------- ------- 2.7619682882369117 0.01803 3.59192 222.04643 ... 14169511631100 18.03 0.15 .. code:: ipython3 # define a set of JD times to propagate the orbital elements to N=1000 times = Time(Time(np.linspace(np.amin(df["mjd"]),np.amax(df["mjd"]),N), format = "mjd").jd, format = "jd") times[0] .. parsed-literal::