LSST Solar System Simulations, June 2021
========================================

Juric, Eggl, Jones, Fedorets, Cornwall, Berres, Chernyavskaya, Moeyens,
Schwamb, et many many al.

This notebook illustrates what is available in the June 2021 version of
the simulated LSST Solar System dataset. Use it as a starting point for
your own experiments.

Getting this notebook
---------------------

This notebook is available from
https://github.com/lsst-sssc/lsst-simulation. Open a terminal and clone
it into your home directory to run it:

::

       git clone https://github.com/lsst-sssc/lsst-simulation

Connecting and inspecting available tables
------------------------------------------

.. code:: ipython3

    import psycopg2 as pg
    import pandas as pd
    import numpy as np
    import matplotlib.pyplot as plt

.. code:: ipython3

    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")

See which tables are available:

.. code:: ipython3

    tables = pd.read_sql("SELECT table_name FROM information_schema.tables WHERE table_schema = 'public'", con)
    tables = tables.set_index("table_name")
    tables




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
        </tr>
        <tr>
          <th>table_name</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>ssobjects</th>
        </tr>
        <tr>
          <th>sssources</th>
        </tr>
        <tr>
          <th>mpcorb</th>
        </tr>
        <tr>
          <th>diasources</th>
        </tr>
      </tbody>
    </table>
    </div>



See how many rows does each table have (this may take a minute or two).

.. code:: ipython3

    %%time
    tables["nrows"] = np.zeros(len(tables), dtype=int)
    for table in tables.index:
        df = pd.read_sql(f"SELECT COUNT(*) FROM {table}", con, params=dict(table=table))
        tables["nrows"].loc[table] = df["count"].iloc[0]


.. parsed-literal::

    CPU times: user 11.1 ms, sys: 1.17 ms, total: 12.3 ms
    Wall time: 1min 11s


.. code:: ipython3

    tables




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>nrows</th>
        </tr>
        <tr>
          <th>table_name</th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>ssobjects</th>
          <td>10556741</td>
        </tr>
        <tr>
          <th>sssources</th>
          <td>1043415800</td>
        </tr>
        <tr>
          <th>mpcorb</th>
          <td>14600302</td>
        </tr>
        <tr>
          <th>diasources</th>
          <td>1043415800</td>
        </tr>
      </tbody>
    </table>
    </div>



Let’s get a feel for the available data, by grabbing the top 5 rows of
each table:

.. code:: ipython3

    pd.read_sql("SELECT * FROM diasources LIMIT 5", con)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>diasourceid</th>
          <th>ccdvisitid</th>
          <th>diaobjectid</th>
          <th>ssobjectid</th>
          <th>_name</th>
          <th>ssobjectreassoctime</th>
          <th>midpointtai</th>
          <th>ra</th>
          <th>rasigma</th>
          <th>decl</th>
          <th>declsigma</th>
          <th>ra_decl_cov</th>
          <th>snr</th>
          <th>filter</th>
          <th>mag</th>
          <th>magsigma</th>
          <th>_v</th>
          <th>_magtrue</th>
          <th>_ratrue</th>
          <th>_dectrue</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>1952839109036039781</td>
          <td>122231</td>
          <td>1537215218520783399</td>
          <td>819889643154482779</td>
          <td>S1005CYCa</td>
          <td>60039.037419</td>
          <td>60039.037419</td>
          <td>129.590684</td>
          <td>0.000009</td>
          <td>12.003733</td>
          <td>0.000009</td>
          <td>0.0</td>
          <td>16.698526</td>
          <td>i</td>
          <td>21.043230</td>
          <td>0.063147</td>
          <td>21.521744</td>
          <td>21.066744</td>
          <td>129.590684</td>
          <td>12.003725</td>
        </tr>
        <tr>
          <th>1</th>
          <td>1486543917817498729</td>
          <td>123098</td>
          <td>-7037468574229910592</td>
          <td>5786782821283451137</td>
          <td>S1005D4qa</td>
          <td>60040.005246</td>
          <td>60040.005246</td>
          <td>124.049298</td>
          <td>0.000010</td>
          <td>-5.136787</td>
          <td>0.000010</td>
          <td>0.0</td>
          <td>13.551331</td>
          <td>i</td>
          <td>21.699242</td>
          <td>0.077302</td>
          <td>22.062727</td>
          <td>21.607727</td>
          <td>124.049304</td>
          <td>-5.136777</td>
        </tr>
        <tr>
          <th>2</th>
          <td>8744810120905494086</td>
          <td>124777</td>
          <td>8686809464008266824</td>
          <td>-9093662608188820544</td>
          <td>S1005D66a</td>
          <td>60041.348960</td>
          <td>60041.348960</td>
          <td>283.045411</td>
          <td>0.000008</td>
          <td>-33.137396</td>
          <td>0.000008</td>
          <td>0.0</td>
          <td>11.022760</td>
          <td>y</td>
          <td>21.476673</td>
          <td>0.094285</td>
          <td>21.766966</td>
          <td>21.463966</td>
          <td>283.045401</td>
          <td>-33.137406</td>
        </tr>
        <tr>
          <th>3</th>
          <td>3452285572810631019</td>
          <td>124834</td>
          <td>-1859933907867084349</td>
          <td>-9093662608188820544</td>
          <td>S1005D66a</td>
          <td>60041.375001</td>
          <td>60041.375001</td>
          <td>283.050372</td>
          <td>0.000008</td>
          <td>-33.138294</td>
          <td>0.000008</td>
          <td>0.0</td>
          <td>11.308365</td>
          <td>y</td>
          <td>21.415295</td>
          <td>0.092001</td>
          <td>21.766676</td>
          <td>21.463676</td>
          <td>283.050357</td>
          <td>-33.138322</td>
        </tr>
        <tr>
          <th>4</th>
          <td>-1529909586803787453</td>
          <td>124861</td>
          <td>3661636550541166030</td>
          <td>-9093662608188820544</td>
          <td>S1005D66a</td>
          <td>60041.387444</td>
          <td>60041.387444</td>
          <td>283.052716</td>
          <td>0.000007</td>
          <td>-33.138754</td>
          <td>0.000007</td>
          <td>0.0</td>
          <td>12.064641</td>
          <td>y</td>
          <td>21.430834</td>
          <td>0.086458</td>
          <td>21.766537</td>
          <td>21.463537</td>
          <td>283.052719</td>
          <td>-33.138756</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    pd.read_sql("SELECT * FROM sssources LIMIT 5", con)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>ssobjectid</th>
          <th>diasourceid</th>
          <th>mpcuniqueid</th>
          <th>eclipticlambda</th>
          <th>eclipticbeta</th>
          <th>galacticl</th>
          <th>galacticb</th>
          <th>phaseangle</th>
          <th>heliocentricdist</th>
          <th>topocentricdist</th>
          <th>...</th>
          <th>heliocentricz</th>
          <th>heliocentricvx</th>
          <th>heliocentricvy</th>
          <th>heliocentricvz</th>
          <th>topocentricx</th>
          <th>topocentricy</th>
          <th>topocentricz</th>
          <th>topocentricvx</th>
          <th>topocentricvy</th>
          <th>topocentricvz</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>819889643154482779</td>
          <td>1952839109036039781</td>
          <td>0</td>
          <td>128.836575</td>
          <td>-6.258677</td>
          <td>213.949948</td>
          <td>29.104339</td>
          <td>18.748798</td>
          <td>2.842328</td>
          <td>2.284481</td>
          <td>...</td>
          <td>0.374356</td>
          <td>-0.004101</td>
          <td>-0.008415</td>
          <td>-0.000889</td>
          <td>-1.424062</td>
          <td>1.721966</td>
          <td>0.475116</td>
          <td>-0.008015</td>
          <td>0.007075</td>
          <td>0.005753</td>
        </tr>
        <tr>
          <th>1</th>
          <td>5786782821283451137</td>
          <td>1486543917817498729</td>
          <td>0</td>
          <td>127.700247</td>
          <td>-24.230133</td>
          <td>227.644419</td>
          <td>16.234911</td>
          <td>17.612328</td>
          <td>3.106055</td>
          <td>2.617551</td>
          <td>...</td>
          <td>-0.341533</td>
          <td>-0.004793</td>
          <td>-0.007634</td>
          <td>-0.003896</td>
          <td>-1.459697</td>
          <td>2.160078</td>
          <td>-0.234359</td>
          <td>-0.008954</td>
          <td>0.007755</td>
          <td>0.002716</td>
        </tr>
        <tr>
          <th>2</th>
          <td>-9093662608188820544</td>
          <td>8744810120905494086</td>
          <td>0</td>
          <td>281.072220</td>
          <td>-10.198162</td>
          <td>2.838461</td>
          <td>-14.577027</td>
          <td>19.033407</td>
          <td>3.053244</td>
          <td>2.785494</td>
          <td>...</td>
          <td>-1.638718</td>
          <td>0.009758</td>
          <td>-0.001617</td>
          <td>-0.000669</td>
          <td>0.526492</td>
          <td>-2.272270</td>
          <td>-1.522687</td>
          <td>0.004808</td>
          <td>0.013617</td>
          <td>0.005900</td>
        </tr>
        <tr>
          <th>3</th>
          <td>-9093662608188820544</td>
          <td>3452285572810631019</td>
          <td>0</td>
          <td>281.076340</td>
          <td>-10.199436</td>
          <td>2.839294</td>
          <td>-14.581200</td>
          <td>19.032654</td>
          <td>3.053253</td>
          <td>2.785145</td>
          <td>...</td>
          <td>-1.638735</td>
          <td>0.009758</td>
          <td>-0.001616</td>
          <td>-0.000669</td>
          <td>0.526617</td>
          <td>-2.271916</td>
          <td>-1.522533</td>
          <td>0.004791</td>
          <td>0.013579</td>
          <td>0.005900</td>
        </tr>
        <tr>
          <th>4</th>
          <td>-9093662608188820544</td>
          <td>-1529909586803787453</td>
          <td>0</td>
          <td>281.078284</td>
          <td>-10.200073</td>
          <td>2.839653</td>
          <td>-14.583186</td>
          <td>19.032290</td>
          <td>3.053257</td>
          <td>2.784979</td>
          <td>...</td>
          <td>-1.638743</td>
          <td>0.009758</td>
          <td>-0.001616</td>
          <td>-0.000669</td>
          <td>0.526676</td>
          <td>-2.271747</td>
          <td>-1.522460</td>
          <td>0.004785</td>
          <td>0.013561</td>
          <td>0.005899</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 29 columns</p>
    </div>



.. code:: ipython3

    pd.read_sql("SELECT * FROM ssobjects LIMIT 5", con)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>ssobjectid</th>
          <th>discoverysubmissiondate</th>
          <th>firstobservationdate</th>
          <th>arc</th>
          <th>numobs</th>
          <th>moid</th>
          <th>moidtrueanomaly</th>
          <th>moideclipticlongitude</th>
          <th>moiddeltav</th>
          <th>uh</th>
          <th>...</th>
          <th>yg12</th>
          <th>yherr</th>
          <th>yg12err</th>
          <th>yh_yg12_cov</th>
          <th>ychi2</th>
          <th>yndata</th>
          <th>maxextendedness</th>
          <th>minextendedness</th>
          <th>medianextendedness</th>
          <th>flags</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>668135024989</td>
          <td>60797.136495</td>
          <td>60790.136495</td>
          <td>0.000000</td>
          <td>1</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>...</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>3033569589766</td>
          <td>63480.233685</td>
          <td>63473.233685</td>
          <td>1.975244</td>
          <td>3</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>...</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>3148445770109</td>
          <td>59937.159828</td>
          <td>59930.159828</td>
          <td>2693.060000</td>
          <td>3</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>...</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>3369984299447</td>
          <td>60203.374981</td>
          <td>60196.374981</td>
          <td>2582.777000</td>
          <td>107</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>...</td>
          <td>0.345189</td>
          <td>0.087608</td>
          <td>0.133565</td>
          <td>0.010005</td>
          <td>0.571543</td>
          <td>7</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>3643818542061</td>
          <td>59917.313961</td>
          <td>59910.313961</td>
          <td>46.793457</td>
          <td>10</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>...</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 55 columns</p>
    </div>



.. code:: ipython3

    pd.read_sql("SELECT * FROM mpcorb LIMIT 5", con)




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>mpcdesignation</th>
          <th>mpcnumber</th>
          <th>ssobjectid</th>
          <th>mpch</th>
          <th>mpcg</th>
          <th>epoch</th>
          <th>tperi</th>
          <th>peri</th>
          <th>node</th>
          <th>incl</th>
          <th>...</th>
          <th>arc</th>
          <th>arcstart</th>
          <th>arcend</th>
          <th>rms</th>
          <th>pertsshort</th>
          <th>pertslong</th>
          <th>computer</th>
          <th>flags</th>
          <th>fulldesignation</th>
          <th>lastincludedobservation</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>SR000001a</td>
          <td>0</td>
          <td>-2658575675934308610</td>
          <td>7.95</td>
          <td>0.15</td>
          <td>54800.0</td>
          <td>52354.06796</td>
          <td>47.02487</td>
          <td>66.70177</td>
          <td>0.0</td>
          <td>...</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>None</td>
          <td>None</td>
          <td>0</td>
          <td>None</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>SR000002a</td>
          <td>0</td>
          <td>1638696702905544284</td>
          <td>11.90</td>
          <td>0.15</td>
          <td>54800.0</td>
          <td>54255.95952</td>
          <td>111.14332</td>
          <td>149.76906</td>
          <td>0.0</td>
          <td>...</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>None</td>
          <td>None</td>
          <td>0</td>
          <td>None</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>SR000003a</td>
          <td>0</td>
          <td>-6148727610239885338</td>
          <td>7.58</td>
          <td>0.15</td>
          <td>54800.0</td>
          <td>32718.27775</td>
          <td>86.93460</td>
          <td>157.72377</td>
          <td>0.0</td>
          <td>...</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>None</td>
          <td>None</td>
          <td>0</td>
          <td>None</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>SR000004a</td>
          <td>0</td>
          <td>2476339031007217136</td>
          <td>12.99</td>
          <td>0.15</td>
          <td>54800.0</td>
          <td>49162.23730</td>
          <td>357.05184</td>
          <td>318.01916</td>
          <td>0.0</td>
          <td>...</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>None</td>
          <td>None</td>
          <td>0</td>
          <td>None</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>SR000006a</td>
          <td>0</td>
          <td>-6439077407999146266</td>
          <td>9.28</td>
          <td>0.15</td>
          <td>54800.0</td>
          <td>-13278.77515</td>
          <td>302.03901</td>
          <td>292.49508</td>
          <td>0.0</td>
          <td>...</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>None</td>
          <td>None</td>
          <td>None</td>
          <td>0</td>
          <td>None</td>
          <td>0.0</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 27 columns</p>
    </div>



Plot a visit
~~~~~~~~~~~~

.. code:: ipython3

    %%time
    sql = """
        SELECT
            ssObjectId, ra, decl, ccdVisitId, midPointTai
        FROM
            diaSources
        WHERE
            ccdVisitId = 126272
    """
    df = pd.read_sql(sql, con)


.. parsed-literal::

    CPU times: user 5.96 ms, sys: 2.46 ms, total: 8.42 ms
    Wall time: 45.4 ms


.. code:: ipython3

    plt.figure(figsize=(10, 10))
    plt.scatter(df["ra"], df["decl"], s=1)
    plt.gca().invert_xaxis()
    plt.xlabel("R.A. (deg)")
    plt.ylabel("Dec (deg)")
    plt.title(f"Visit {df['ccdvisitid'].iloc[0]} (MJD {df['midpointtai'].iloc[0]:.5f})")




.. parsed-literal::

    Text(0.5, 1.0, 'Visit 126272 (MJD 60043.17665)')




.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_20_1.png


Now grab the same data in ecliptic coordinates:

.. code:: ipython3

    %%time
    sql = """
        SELECT
            eclipticLambda as lon, eclipticBeta as lat, ccdVisitId, midPointTAI
        FROM
            diaSources JOIN ssSources USING(diaSourceId)
        WHERE
            ccdVisitId = 126272
    """
    df = pd.read_sql(sql, con)


.. parsed-literal::

    CPU times: user 59 ms, sys: 6.12 ms, total: 65.1 ms
    Wall time: 25.3 s


.. code:: ipython3

    plt.figure(figsize=(10, 10))
    plt.scatter(df["lon"], df["lat"], s=1)
    plt.gca().invert_xaxis()
    plt.xlabel(r"$\lambda$ (deg)")
    plt.ylabel(r"$\beta$ (deg)")
    plt.title(f"Visit {df['ccdvisitid'].iloc[0]} (MJD {df['midpointtai'].iloc[0]:.5f})");



.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_23_0.png


Plot a night
~~~~~~~~~~~~

.. code:: ipython3

    %%time
    sql = """
        SELECT
            ssObjectId, ra, decl, ccdVisitId, midPointTai
        FROM
            diaSources
        WHERE
            midPointTai BETWEEN 60000-0.5 AND 60000+0.5
    """
    df = pd.read_sql(sql, con)


.. parsed-literal::

    CPU times: user 836 ms, sys: 148 ms, total: 984 ms
    Wall time: 4.36 s


.. code:: ipython3

    plt.figure(figsize=(14, 7))
    plt.scatter(df["ra"], df["decl"], s=1)
    plt.gca().invert_xaxis()
    plt.xlabel("R.A. (deg)")
    plt.ylabel("Dec (deg)")
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.title(f"Night MJD 60000")




.. parsed-literal::

    Text(0.5, 1.0, 'Night MJD 60000')




.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_26_1.png


First month of the survey
~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    t0 = pd.read_sql("SELECT MIN(midPointTai) from diaSources", con)["min"].iloc[0]
    t0




.. parsed-literal::

    59853.98564424414



.. code:: ipython3

    %%time
    sql = """
        SELECT
            ra, decl
        FROM
            diaSources
        WHERE
            midPointTai BETWEEN %(start)s AND %(end)s
    """
    df = pd.read_sql(sql, con, params=dict(start=t0, end=t0+30))


.. parsed-literal::

    CPU times: user 12.6 s, sys: 1.97 s, total: 14.6 s
    Wall time: 23.1 s


.. code:: ipython3

    plt.figure(figsize=(14, 7))
    plt.scatter(df["ra"], df["decl"], s=.1)
    plt.gca().invert_xaxis()
    plt.xlabel("R.A. (deg)")
    plt.ylabel("Dec (deg)")
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.title(f"First month of the survey");



.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_30_0.png


Plot a phase curve
~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    sql = """
        SELECT
            mpcdesignation, ssObjects.ssObjectId, mag, magSigma, filter, midPointTai as mjd, ra, decl, phaseAngle,
            topocentricDist, heliocentricDist
        FROM
            mpcorb
            JOIN ssObjects USING (ssobjectid)
            JOIN diaSources USING (ssobjectid)
            JOIN ssSources USING (diaSourceid)
        WHERE
            mpcdesignation = 'S00001vAa' and filter='r'
    """
    df = pd.read_sql(sql, con)
    # Distance correction
    df["cmag"] = df["mag"] - 5*np.log10(df["topocentricdist"]*df["heliocentricdist"])
    df.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>mpcdesignation</th>
          <th>ssobjectid</th>
          <th>mag</th>
          <th>magsigma</th>
          <th>filter</th>
          <th>mjd</th>
          <th>ra</th>
          <th>decl</th>
          <th>phaseangle</th>
          <th>topocentricdist</th>
          <th>heliocentricdist</th>
          <th>cmag</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>S00001vAa</td>
          <td>8050632269120433289</td>
          <td>20.652550</td>
          <td>0.013050</td>
          <td>r</td>
          <td>62352.235856</td>
          <td>338.662138</td>
          <td>-15.892054</td>
          <td>14.838730</td>
          <td>0.595123</td>
          <td>1.578381</td>
          <td>20.788458</td>
        </tr>
        <tr>
          <th>1</th>
          <td>S00001vAa</td>
          <td>8050632269120433289</td>
          <td>20.670551</td>
          <td>0.016261</td>
          <td>r</td>
          <td>62352.244968</td>
          <td>338.657501</td>
          <td>-15.893365</td>
          <td>14.830055</td>
          <td>0.595161</td>
          <td>1.578451</td>
          <td>20.806224</td>
        </tr>
        <tr>
          <th>2</th>
          <td>S00001vAa</td>
          <td>8050632269120433289</td>
          <td>20.605118</td>
          <td>0.010641</td>
          <td>r</td>
          <td>62360.313185</td>
          <td>334.722045</td>
          <td>-16.968608</td>
          <td>7.679678</td>
          <td>0.635753</td>
          <td>1.639855</td>
          <td>20.514649</td>
        </tr>
        <tr>
          <th>3</th>
          <td>S00001vAa</td>
          <td>8050632269120433289</td>
          <td>20.660650</td>
          <td>0.018832</td>
          <td>r</td>
          <td>62251.421038</td>
          <td>332.847465</td>
          <td>-15.415406</td>
          <td>96.779530</td>
          <td>0.302469</td>
          <td>0.924456</td>
          <td>23.427812</td>
        </tr>
        <tr>
          <th>4</th>
          <td>S00001vAa</td>
          <td>8050632269120433289</td>
          <td>22.143944</td>
          <td>0.040766</td>
          <td>r</td>
          <td>61286.975565</td>
          <td>242.153201</td>
          <td>-24.872374</td>
          <td>51.970360</td>
          <td>0.893512</td>
          <td>1.272925</td>
          <td>21.864428</td>
        </tr>
      </tbody>
    </table>
    </div>



Now make a plot:

.. code:: ipython3

    plt.figure(figsize=(14, 6))
    plt.errorbar(df["phaseangle"], df["cmag"], df["magsigma"], ls='none')
    plt.gca().invert_yaxis()
    plt.xlabel("Phase Angle (deg)")
    plt.ylabel("Distance corrected mag_r (mag)")
    plt.title(f'Phase curve for {df["mpcdesignation"].iloc[0]}, r band');



.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_34_0.png


Now grab our (H, G) fit, and overplot it

.. code:: ipython3

    ssoId = int(df["ssobjectid"].iloc[0])
    hg = pd.read_sql("SELECT rH, rG12, rHErr, rG12Err, rChi2 FROM ssObjects WHERE ssObjectId=%(ssoId)s", con, params=dict(ssoId=ssoId))
    hg




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>rh</th>
          <th>rg12</th>
          <th>rherr</th>
          <th>rg12err</th>
          <th>rchi2</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>19.9568</td>
          <td>0.150825</td>
          <td>0.009669</td>
          <td>0.004122</td>
          <td>1.167215</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    plt.figure(figsize=(14, 6))
    plt.errorbar(df["phaseangle"], df["cmag"], df["magsigma"], ls='none')
    plt.gca().invert_yaxis()
    plt.xlabel("Phase Angle (deg)")
    plt.ylabel("Distance corrected mag_r (mag)")
    plt.title(f'Phase curve for {df["mpcdesignation"].iloc[0]}, r band')
    
    from sbpy.photometry import HG
    H, G, sigmaH, sigmaG, chi2dof = hg.iloc[0]
    _ph = sorted(df["phaseangle"])
    _mag = HG.evaluate(np.deg2rad(_ph), H, G)
    plt.plot(_ph, _mag)
    print(f"H={H:.2f}±{sigmaH:.3}, G={G:.2f}±{sigmaG:.3}, χ2/dof={chi2dof:.3f}")


.. parsed-literal::

    H=19.96±0.00967, G=0.15±0.00412, χ2/dof=1.167



.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_37_1.png


Look at the input population
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This is just S3M, so it should correspond to the plots from the Grav et
al. 2011 paper.

.. code:: ipython3

    sql = """
        SELECT
            FLOOR(mpcH*10)/10 AS binH, count(*)
        FROM
            mpcorb
        GROUP BY binH
    """
    df = pd.read_sql(sql, con)
    plt.figure(figsize=(10, 5))
    plt.plot(df["binh"], df["count"])
    plt.xlabel("H (mag)")
    plt.ylabel("Objects in bin")




.. parsed-literal::

    Text(0, 0.5, 'Objects in bin')




.. image:: solsys-sims-2021-demo_files/solsys-sims-2021-demo_40_1.png