End-to-End Validation Demo
This notebook implements an independent verification of the results from Sorcha on a one month LSST simulation for two objects.
The reference ephemerides for these objects (as CSV files) come from JPL Horizons queries. The reference query uses dates from 2023-10-08 UT to 2023-11-08, with one minute timesteps, for the objects 2011 OB60 (a TNO) and 2010 TU149 (MBA), as seen from the X05 observatory.
[1]:
import numpy as np
import astropy.table as tb
import pandas as pd
import sqlite3
[2]:
import matplotlib.pyplot as pl
[3]:
#read pointing database and restrict to our time baseline
con = sqlite3.connect('../../src/sorcha/data/demo/baseline_v2.0_1yr.db')
pointing = pd.read_sql("SELECT observationId, observationStartMJD as observationStartMJD_TAI, visitTime, visitExposureTime, filter, seeingFwhmGeom as seeingFwhmGeom_arcsec, seeingFwhmEff as seeingFwhmEff_arcsec, fiveSigmaDepth as fieldFiveSigmaDepth_mag , fieldRA as fieldRA_deg, fieldDec as fieldDec_deg, rotSkyPos as fieldRotSkyPos_deg FROM observations order by observationId", con=con, )
pointing
[3]:
| observationId | observationStartMJD_TAI | visitTime | visitExposureTime | filter | seeingFwhmGeom_arcsec | seeingFwhmEff_arcsec | fieldFiveSigmaDepth_mag | fieldRA_deg | fieldDec_deg | fieldRotSkyPos_deg | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 60218.001806 | 34.0 | 30.0 | y | 0.667370 | 0.748626 | 22.370557 | 310.024480 | -60.812928 | 62.750775 |
| 1 | 1 | 60218.002254 | 34.0 | 30.0 | y | 0.678175 | 0.761770 | 22.338327 | 310.601871 | -63.561425 | 61.084250 |
| 2 | 2 | 60218.002703 | 34.0 | 30.0 | y | 0.690895 | 0.777245 | 22.295487 | 311.292611 | -66.317774 | 60.726189 |
| 3 | 3 | 60218.003152 | 34.0 | 30.0 | y | 1.657640 | 1.953333 | 21.270421 | 312.140731 | -69.082666 | 60.656781 |
| 4 | 4 | 60218.003624 | 34.0 | 30.0 | y | 1.713416 | 2.021188 | 21.205276 | 304.170163 | -73.375442 | 49.095537 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 216228 | 216228 | 60583.103691 | 34.0 | 30.0 | g | 0.988866 | 1.139739 | 24.422508 | 8.496364 | -29.615770 | 86.518379 |
| 216229 | 216229 | 60583.104140 | 34.0 | 30.0 | g | 1.003682 | 1.157764 | 24.397505 | 11.808603 | -29.898183 | 86.777460 |
| 216230 | 216230 | 60583.104587 | 34.0 | 30.0 | g | 1.014421 | 1.170829 | 24.374331 | 13.363397 | -27.303043 | 91.147723 |
| 216231 | 216231 | 60583.105032 | 34.0 | 30.0 | g | 1.026912 | 1.186025 | 24.343872 | 14.847339 | -24.711444 | 94.995634 |
| 216232 | 216232 | 60583.105481 | 34.0 | 30.0 | g | 1.046555 | 1.209921 | 24.309309 | 18.100120 | -24.787034 | 94.983283 |
216233 rows × 11 columns
[4]:
import astropy.time as ti
pointing = pointing[pointing['observationStartMJD_TAI'] > 60225]
pointing = pointing[pointing['observationStartMJD_TAI'] < 60255]
pointing['MJD_TDB'] = ti.Time(pointing['observationStartMJD_TAI'] + pointing['visitTime']/(2*86400.0), format='mjd', scale='tai').tdb.mjd
pointing
[4]:
| observationId | observationStartMJD_TAI | visitTime | visitExposureTime | filter | seeingFwhmGeom_arcsec | seeingFwhmEff_arcsec | fieldFiveSigmaDepth_mag | fieldRA_deg | fieldDec_deg | fieldRotSkyPos_deg | MJD_TDB | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 5398 | 5398 | 60225.000341 | 34.0 | 30.0 | r | 0.693873 | 0.780867 | 24.336680 | 325.777117 | -40.076878 | 110.946937 | 60225.000910 |
| 5399 | 5399 | 60225.000800 | 34.0 | 30.0 | r | 0.696867 | 0.784510 | 24.335598 | 325.609267 | -43.105406 | 102.956598 | 60225.001369 |
| 5400 | 5400 | 60225.001260 | 34.0 | 30.0 | r | 0.684560 | 0.769538 | 24.347565 | 332.965064 | -46.169394 | 107.885215 | 60225.001829 |
| 5401 | 5401 | 60225.001709 | 34.0 | 30.0 | r | 0.689428 | 0.775460 | 24.340440 | 333.038252 | -49.111447 | 102.767992 | 60225.002278 |
| 5402 | 5402 | 60225.002157 | 34.0 | 30.0 | r | 0.695373 | 0.782692 | 24.327080 | 333.059351 | -52.057748 | 98.135653 | 60225.002726 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 24474 | 24474 | 60251.366258 | 34.0 | 30.0 | i | 0.699315 | 0.787488 | 23.270283 | 81.912982 | -27.019050 | 234.493172 | 60251.366827 |
| 24475 | 24475 | 60251.366706 | 34.0 | 30.0 | i | 0.706234 | 0.795905 | 23.224385 | 79.796144 | -24.525717 | 228.806467 | 60251.367276 |
| 24476 | 24476 | 60251.367202 | 34.0 | 30.0 | i | 0.702517 | 0.791383 | 23.199502 | 83.540274 | -20.909984 | 218.028671 | 60251.367771 |
| 24477 | 24477 | 60251.367658 | 34.0 | 30.0 | i | 0.745439 | 0.843600 | 23.088602 | 78.410579 | -19.065750 | 218.537314 | 60251.368227 |
| 24478 | 24478 | 60251.368161 | 34.0 | 30.0 | i | 0.742031 | 0.839454 | 23.098674 | 76.814813 | -25.030965 | 230.902390 | 60251.368731 |
19081 rows × 12 columns
Now we have our reference pointings, on to the objects
[5]:
# read our objects
obj = {}
obj['2011 OB60'] = tb.Table.read('2011ob60.txt', format='csv')
obj['2010 TU149'] = tb.Table.read('2010tu149.txt', format='csv')
for i in obj:
obj[i]['MJD'] = ti.Time(obj[i]['JDUT'],format='jd', scale='utc').tdb.mjd
[6]:
obj['2010 TU149']
[6]:
Table length=19908
| JDUT | col1 | _1 | RA | DEC | APmag | S-brt | r | rdot | delta | deldot | S-O-T | /r | S-T-O | O-P-T | PsAng | PsAMV | Sky_motion | Sky_mot_PA | RelVel-ANG | Lun_Sky_Brt | sky_SNR | _2 | MJD |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| float64 | str1 | str1 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | str2 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | str5 | str7 | int64 | float64 |
| 2460225.5 | A | -- | 19.7210206 | 4.856703608 | 21.416 | 6.527 | 1.696279513621 | -20.6503447 | 0.69968470063905 | -22.6844214 | 173.2633 | /L | 3.9648 | 2.49 | 276.388 | 248.881 | 1.6304974 | 246.66601 | -58.70543 | n.a. | n.a. | -- | 60225.00080072161 |
| 2460225.500694444 | A | -- | 19.720603299 | 4.856523883 | 21.416 | 6.527 | 1.696271230638 | -20.6504168 | 0.69967560204437 | -22.6836065 | 173.2643 | /L | 3.9642 | 2.4896 | 276.393 | 248.881 | 1.6307293 | 246.66868 | -58.70123 | n.a. | n.a. | -- | 60225.00149516555 |
| 2460225.501388889 | A | -- | 19.72018593 | 4.856344153 | 21.416 | 6.527 | 1.696262947626 | -20.6504889 | 0.69966650377801 | -22.6827842 | 173.2653 | /L | 3.9637 | 2.4892 | 276.397 | 248.881 | 1.6309611 | 246.67133 | -58.69703 | n.a. | n.a. | -- | 60225.002189610415 |
| 2460225.502083333 | A | -- | 19.719768494 | 4.856164415 | 21.416 | 6.527 | 1.696254664586 | -20.6505609 | 0.69965740584296 | -22.6819545 | 173.2663 | /L | 3.9631 | 2.4888 | 276.402 | 248.881 | 1.6311926 | 246.67398 | -58.69282 | n.a. | n.a. | -- | 60225.00288405482 |
| 2460225.502777778 | A | -- | 19.71935099 | 4.855984672 | 21.415 | 6.527 | 1.696246381516 | -20.650633 | 0.69964830824217 | -22.6811174 | 173.2673 | /L | 3.9625 | 2.4884 | 276.407 | 248.881 | 1.6314239 | 246.67662 | -58.6886 | n.a. | n.a. | -- | 60225.00357849969 |
| 2460225.503472222 | A | -- | 19.718933419 | 4.855804922 | 21.415 | 6.527 | 1.696238098417 | -20.6507051 | 0.69963921097862 | -22.6802729 | 173.2683 | /L | 3.962 | 2.4879 | 276.411 | 248.881 | 1.6316551 | 246.67926 | -58.68438 | n.a. | n.a. | -- | 60225.004272943625 |
| 2460225.504166667 | A | -- | 19.71851578 | 4.855625166 | 21.415 | 6.526 | 1.69622981529 | -20.6507771 | 0.69963011405528 | -22.679421 | 173.2693 | /L | 3.9614 | 2.4875 | 276.416 | 248.881 | 1.631886 | 246.68188 | -58.68016 | n.a. | n.a. | -- | 60225.00496738849 |
| 2460225.504861111 | A | -- | 19.718098075 | 4.855445404 | 21.415 | 6.526 | 1.696221532134 | -20.6508492 | 0.69962101747509 | -22.6785618 | 173.2703 | /L | 3.9608 | 2.4871 | 276.421 | 248.881 | 1.6321167 | 246.68451 | -58.67593 | n.a. | n.a. | -- | 60225.00566183243 |
| 2460225.505555556 | A | -- | 19.717680302 | 4.855265635 | 21.415 | 6.526 | 1.696213248949 | -20.6509213 | 0.69961192124103 | -22.6776952 | 173.2713 | /L | 3.9603 | 2.4867 | 276.425 | 248.881 | 1.6323472 | 246.68712 | -58.67169 | n.a. | n.a. | -- | 60225.006356277765 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2460256.49375 | N | -- | 350.607204357 | -7.190676285 | 21.078 | 7.017 | 1.294842964898 | -24.3811701 | 0.44976224719555 | -5.9900853 | 123.5297 | /T | 39.6396 | 17.1575 | 64.949 | 25.286 | 3.1818073 | 249.62846 | -19.10278 | n.a. | n.a. | -- | 60255.99455072442 |
| 2460256.494444444 | N | -- | 350.606369072 | -7.190983606 | 21.078 | 7.017 | 1.294833186 | -24.3812663 | 0.44975984486785 | -5.9881083 | 123.5281 | /T | 39.6408 | 17.1579 | 64.949 | 25.314 | 3.1818856 | 249.62984 | -19.09659 | n.a. | n.a. | -- | 60255.99524516836 |
| 2460256.495138889 | N | -- | 350.605533758 | -7.191290915 | 21.078 | 7.017 | 1.294823407064 | -24.3813625 | 0.44975744333359 | -5.9861291 | 123.5265 | /T | 39.642 | 17.1582 | 64.95 | 25.342 | 3.1819626 | 249.63121 | -19.0904 | n.a. | n.a. | -- | 60255.99593961369 |
| 2460256.495833333 | N | -- | 350.604698417 | -7.191598211 | 21.078 | 7.017 | 1.294813628089 | -24.3814587 | 0.44975504259367 | -5.9841477 | 123.5249 | /T | 39.6433 | 17.1586 | 64.95 | 25.37 | 3.1820383 | 249.63257 | -19.08422 | n.a. | n.a. | -- | 60255.99663405763 |
| 2460256.496527778 | N | -- | 350.603863048 | -7.191905495 | 21.078 | 7.017 | 1.294803849075 | -24.381555 | 0.44975264264894 | -5.9821642 | 123.5234 | /T | 39.6445 | 17.1589 | 64.95 | 25.398 | 3.1821127 | 249.63393 | -19.07803 | n.a. | n.a. | -- | 60255.9973285025 |
| 2460256.497222222 | N | -- | 350.603027651 | -7.192212766 | 21.078 | 7.017 | 1.294794070023 | -24.3816512 | 0.44975024350027 | -5.9801786 | 123.5218 | /T | 39.6457 | 17.1593 | 64.95 | 25.427 | 3.1821857 | 249.63528 | -19.07184 | n.a. | n.a. | -- | 60255.998022946435 |
| 2460256.497916667 | N | -- | 350.602192228 | -7.192520025 | 21.078 | 7.017 | 1.294784290932 | -24.3817474 | 0.44974784514848 | -5.9781908 | 123.5202 | /T | 39.6469 | 17.1596 | 64.951 | 25.455 | 3.1822575 | 249.63662 | -19.06566 | n.a. | n.a. | -- | 60255.9987173913 |
| 2460256.498611111 | N | -- | 350.601356778 | -7.192827271 | 21.078 | 7.017 | 1.294774511803 | -24.3818436 | 0.44974544759442 | -5.9762011 | 123.5186 | /T | 39.6482 | 17.16 | 64.951 | 25.483 | 3.182328 | 249.63795 | -19.05947 | n.a. | n.a. | -- | 60255.999411835706 |
| 2460256.499305556 | N | -- | 350.600521302 | -7.193134505 | 21.078 | 7.017 | 1.294764732635 | -24.3819399 | 0.4497430508389 | -5.9742093 | 123.517 | /T | 39.6494 | 17.1603 | 64.951 | 25.511 | 3.1823971 | 249.63928 | -19.05329 | n.a. | n.a. | -- | 60256.000106280575 |
| 2460256.5 | N | -- | 350.599685801 | -7.193441726 | 21.078 | 7.017 | 1.294754953429 | -24.3820361 | 0.44974065488271 | -5.9722155 | 123.5155 | /T | 39.6506 | 17.1607 | 64.951 | 25.539 | 3.182465 | 249.64059 | -19.04711 | n.a. | n.a. | -- | 60256.00080072451 |
Because of the very dense reference JPL ephemerides, and because we’re using an asteroid and a TNO, I’m assuming that the functions are ~linear in the <1 minute timescale needed for the interpolation (Taylor’s theorem is our friend here). I think this also means that I don’t have to bother with light time correction - linearity once again.
[7]:
interpolated = {}
for i in obj:
interpolated[i] = tb.Table()
for j in pointing.keys():
interpolated[i][j] = pointing[j]
interpolated[i]['RA_INTERP'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['RA'])
interpolated[i]['DEC_INTERP'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['DEC'])
interpolated[i]['delta'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['delta'])
interpolated[i]['r'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['r'])
interpolated[i]['S-T-O'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['S-T-O'])
interpolated[i]['rdot'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['rdot'])
interpolated[i]['deldot'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['deldot'])
interpolated[i]['Sky_motion'] = np.interp(pointing['MJD_TDB'], obj[i]['MJD'], obj[i]['Sky_motion'])
interpolated[i]['ExpTime'] = pointing['visitExposureTime']
interpolated[i]['seeing'] = pointing['seeingFwhmEff_arcsec']
[8]:
interpolated['2011 OB60']
[8]:
Table length=19081
| observationId | observationStartMJD_TAI | visitTime | visitExposureTime | filter | seeingFwhmGeom_arcsec | seeingFwhmEff_arcsec | fieldFiveSigmaDepth_mag | fieldRA_deg | fieldDec_deg | fieldRotSkyPos_deg | MJD_TDB | RA_INTERP | DEC_INTERP | delta | r | S-T-O | rdot | deldot | Sky_motion | ExpTime | seeing |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| int64 | float64 | float64 | float64 | object | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
| 5398 | 60225.00034115714 | 34.0 | 30.0 | r | 0.6938726423855437 | 0.7808669615395909 | 24.33668024295265 | 325.7771174872412 | -40.076877551247264 | 110.94693749742562 | 60225.00091039727 | 1.9873554445353354 | -11.893477087167534 | 35.97110419504151 | 36.90820805260721 | 0.5412 | -0.4212427842066933 | 8.251323402069083 | 0.05184701062597272 | 30.0 | 0.7808669615395909 |
| 5399 | 60225.00080010896 | 34.0 | 30.0 | r | 0.6968672883684072 | 0.7845100831732448 | 24.335598043680374 | 325.6092673713696 | -43.10540619671949 | 102.95659764864814 | 60225.00136934909 | 1.9873465251489812 | -11.89348088332604 | 35.97110638231554 | 36.90820794094966 | 0.5412 | -0.42124271811758285 | 8.25217581941642 | 0.05184872894284594 | 30.0 | 0.7845100831732448 |
| 5400 | 60225.001259967845 | 34.0 | 30.0 | r | 0.6845598801000351 | 0.7695375670316729 | 24.347564772301894 | 332.9650642821256 | -46.16939439980464 | 107.8852146179335 | 60225.00182920798 | 1.9873375876621897 | -11.893484687464527 | 35.971108574159835 | 36.908207829076346 | 0.5412 | -0.42124260379583894 | 8.253032999167933 | 0.05185049875617452 | 30.0 | 0.7695375670316729 |
| 5401 | 60225.001708804055 | 34.0 | 30.0 | r | 0.6894283327924795 | 0.7754602588716295 | 24.340440054873458 | 333.0382524231313 | -49.11144728559235 | 102.76799182310332 | 60225.00227804419 | 1.9873288641020699 | -11.893488400467604 | 35.971110713623865 | 36.90820771988512 | 0.5412 | -0.42124248726553615 | 8.253871540922693 | 0.051852231096060054 | 30.0 | 0.7754602588716295 |
| 5402 | 60225.002156650444 | 34.0 | 30.0 | r | 0.6953728388032058 | 0.7826920180087662 | 24.32708021818426 | 333.05935059163113 | -52.05774779228485 | 98.13565300639495 | 60225.002725890576 | 1.9873201592575491 | -11.893492104766532 | 35.971112848640665 | 36.908207610929466 | 0.5412 | -0.4212424227756523 | 8.254711392679924 | 0.05185390783303985 | 30.0 | 0.7826920180087662 |
| 5403 | 60225.00261591108 | 34.0 | 30.0 | r | 0.6998916012261921 | 0.7881892958956108 | 24.3168344571622 | 337.26060311437726 | -50.34499990813282 | 106.20159055210118 | 60225.00318515121 | 1.9873112321232955 | -11.89349590347513 | 35.97111503829866 | 36.90820749920124 | 0.5412 | -0.4212423132842921 | 8.25557542423461 | 0.051855670662056576 | 30.0 | 0.7881892958956108 |
| 5404 | 60225.00306292132 | 34.0 | 30.0 | r | 0.7063866806269503 | 0.7960908523442217 | 24.305780110279173 | 337.5200318330988 | -53.24619200825491 | 102.2521992381796 | 60225.00363216145 | 1.9873025428145017 | -11.893499600778838 | 35.97111716970597 | 36.90820739045473 | 0.5412077272999876 | -0.42124219227270004 | 8.256418298540236 | 0.05185740090979968 | 30.0 | 0.7960908523442217 |
| 5405 | 60225.00350911682 | 34.0 | 30.0 | r | 0.7141020500337235 | 0.8054769465130456 | 24.286376821856962 | 337.77910898067853 | -56.15688383383206 | 98.71121658386097 | 60225.00407835695 | 1.987293868767728 | -11.893503290782588 | 35.971119297528354 | 36.90820728190064 | 0.5412719794983109 | -0.4212421280205017 | 8.25726308644379 | 0.05185907146695608 | 30.0 | 0.8054769465130456 |
| 5406 | 60225.00395464612 | 34.0 | 30.0 | r | 0.723203397473208 | 0.8165491453445354 | 24.267378428106024 | 338.06419298018744 | -59.08511327945817 | 95.57397825850781 | 60225.00452388625 | 1.9872852073169596 | -11.893506975635537 | 35.971121422363225 | 36.908207173512324 | 0.5413 | -0.42124202772856767 | 8.25810888894636 | 0.05186073952862049 | 30.0 | 0.8165491453445354 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 24469 | 60251.36396536924 | 34.0 | 30.0 | i | 0.7658574955370132 | 0.8684397755924735 | 23.265012163462085 | 84.92763726650861 | -26.406384933962286 | 231.3041995337763 | 60251.36453461151 | 1.5372784971908364 | -12.052164179387676 | 36.193739703064416 | 36.90183866888542 | 1.0846 | -0.41541992232017727 | 20.5123439581853 | 0.03534005976762277 | 30.0 | 0.8684397755924735 |
| 24470 | 60251.36441506614 | 34.0 | 30.0 | i | 0.688792880588721 | 0.7746872026626777 | 23.371623125512052 | 87.92121132653646 | -25.733660553639123 | 226.89862145028235 | 60251.36498430841 | 1.5372722013531812 | -12.052165759018783 | 36.19374503056942 | 36.9018385610013 | 1.0846 | -0.4154198151276112 | 20.512302352622147 | 0.03533723280744675 | 30.0 | 0.7746872026626777 |
| 24471 | 60251.36486418307 | 34.0 | 30.0 | i | 0.6942765348829685 | 0.781358314942784 | 23.328659898469283 | 85.70114045455082 | -23.33851671640359 | 221.3639781858851 | 60251.365433425344 | 1.5372659139983802 | -12.052167336390221 | 36.19375035119758 | 36.90183845325558 | 1.0846 | -0.41541969289092984 | 20.512258632179012 | 0.03533439430998308 | 30.0 | 0.781358314942784 |
| 24472 | 60251.36531309651 | 34.0 | 30.0 | i | 0.6993078263098574 | 0.7874791074329166 | 23.30321589967497 | 82.75731320279174 | -23.961345452242394 | 225.65656971715728 | 60251.36588233878 | 1.537259630650622 | -12.052168913045074 | 36.19375566939229 | 36.90183834555305 | 1.0846 | -0.4154196282474344 | 20.51221073134888 | 0.035331614639678544 | 30.0 | 0.7874791074329166 |
| 24473 | 60251.36580192323 | 34.0 | 30.0 | i | 0.6892814588848022 | 0.7752815801518275 | 23.31835430181076 | 87.06708091095462 | -28.81907452563399 | 239.52450105117353 | 60251.3663711655 | 1.5372527890608927 | -12.052170629883062 | 36.193761460423815 | 36.90183822827872 | 1.0846421436269913 | -0.415419515712746 | 20.512155410800375 | 0.03532858782403937 | 30.0 | 0.7752815801518275 |
| 24474 | 60251.36625772165 | 34.0 | 30.0 | i | 0.6993148289856683 | 0.7874876265032461 | 23.270283068893956 | 81.91298173932086 | -27.01905023032415 | 234.49317169674072 | 60251.366826963924 | 1.5372464100725505 | -12.052172230720833 | 36.19376686016648 | 36.901838118930804 | 1.0847 | -0.4154193844427362 | 20.512101277017596 | 0.0353257577401957 | 30.0 | 0.7874876265032461 |
| 24475 | 60251.36670644113 | 34.0 | 30.0 | i | 0.7062338782491898 | 0.7959049613737103 | 23.22438540541309 | 79.79614407028939 | -24.52571697093717 | 228.8064668938785 | 60251.3672756834 | 1.5372401307361665 | -12.052173806693403 | 36.19377217601759 | 36.90183801128134 | 1.0847 | -0.41541925521169304 | 20.51204376920342 | 0.035322914657247656 | 30.0 | 0.7959049613737103 |
| 24476 | 60251.36720158276 | 34.0 | 30.0 | i | 0.7025167006889115 | 0.7913828475534203 | 23.199501962605517 | 83.5402742535662 | -20.909983893749537 | 218.0286713384419 | 60251.36777082503 | 1.5372332026379798 | -12.052175545709959 | 36.193778041809665 | 36.90183789249052 | 1.0847 | -0.4154191563054547 | 20.51197712215289 | 0.035319821134552606 | 30.0 | 0.7913828475534203 |
| 24477 | 60251.367658208954 | 34.0 | 30.0 | i | 0.7454393324967294 | 0.8436001611882353 | 23.08860205706436 | 78.41057860401273 | -19.065749943842853 | 218.53731409614824 | 60251.368227451225 | 1.537226814053712 | -12.05217714945507 | 36.19378345131093 | 36.90183778293838 | 1.0847 | -0.41541908110249537 | 20.51191310164384 | 0.0353169937036502 | 30.0 | 0.8436001611882353 |
| 24478 | 60251.36816127136 | 34.0 | 30.0 | i | 0.742030868233932 | 0.839453610990185 | 23.098674365772652 | 76.81481324316428 | -25.03096543470499 | 230.90239004531708 | 60251.36873051363 | 1.5372197764161373 | -12.052178916289659 | 36.193789410895214 | 36.90183766225177 | 1.0847 | -0.4154189362206101 | 20.511837980386353 | 0.03531387874311783 | 30.0 | 0.839453610990185 |
Now we need to see whether the detections are inside the circular footprint we’re using as a proxy for the LSST camera:
[9]:
import numba
@numba.njit
def haversine(ra1, dec1, ra2, dec2):
# adapted from Mario's implementation. in his wise words:
# "SkyCoord is slow AF"
dlon = np.deg2rad(ra2 - ra1)
dlat = np.deg2rad(dec2 - dec1)
a = np.sin(dlat / 2.0) ** 2 + np.cos(np.deg2rad(dec1)) * np.cos(np.deg2rad(dec2)) * np.sin(dlon / 2.0) ** 2
c = 2 * np.arcsin(np.sqrt(a))
return np.degrees(c)
[10]:
for i in interpolated:
interpolated[i]['DISTANCE'] = haversine(interpolated[i]['RA_INTERP'], interpolated[i]['DEC_INTERP'], interpolated[i]['fieldRA_deg'], interpolated[i]['fieldDec_deg'])
[11]:
insideccd = {}
for i in interpolated:
insideccd[i] = interpolated[i][interpolated[i]['DISTANCE'] < 1.75]
insideccd[i].sort('observationStartMJD_TAI')
Now onto magnitudes and colors
[12]:
phot = tb.Table.read('../../demo/sspp_testset_colours.txt', format='ascii')
mags = {}
for i in phot:
b = ' '.join(i['ObjID'].split('_'))
mags[b] = {}
mags[b]['H_r'] = i['H_r']
mags[b]['GS'] = i['GS']
mags[b]['r-r'] = 0
for j in 'ugizy':
mags[b][f'{j}-r'] = i[f'{j}-r']
[13]:
for i in insideccd:
colors = np.array([mags[i][f'{j}-r'] for j in insideccd[i]['filter']])
insideccd[i]['mag'] = mags[i]['H_r'] + colors + 5 * np.log10(insideccd[i]['r'] * insideccd[i]['delta'])
Loading the (pre-computed) Sorcha results
[14]:
verification = tb.Table.read('verification_output.csv')
v = {}
v['2011 OB60'] = verification[verification['ObjID'] == '2011_OB60']
v['2010 TU149'] = verification[verification['ObjID'] == '2010_TU149']
[15]:
for i in v:
v[i].sort('fieldMJD_TAI')
If everything went according to plan, we should have the same number of detections in both datasets
[16]:
for i in v:
print(i,len(v[i]), len(insideccd[i]))
2011 OB60 32 32
2010 TU149 18 18
And they should be at the same times:
[17]:
for i in v:
print(i, np.mean(v[i]['fieldMJD_TAI'] - insideccd[i]['observationStartMJD_TAI']))
2011 OB60 0.00019567419099075778
2010 TU149 0.0001957947562105902
(note that this 0.00019567419099075778 offset is the ~17 seconds we added earlier to correct for exposure midpoints - so these are the same pointings)
Finally, let’s visualize the results:
[18]:
fig, ax = pl.subplots(1,3, sharey=True)
fig.subplots_adjust(wspace=0.2)
ax[0].hist(3600*1000*(v['2011 OB60']['RA_deg'] - insideccd['2011 OB60']['RA_INTERP']), histtype='step', color='xkcd:cerulean blue')
ax[0].hist(3600*1000*(v['2010 TU149']['RA_deg'] - insideccd['2010 TU149']['RA_INTERP']), histtype='step', color='xkcd:scarlet')
ax[0].set_ylabel('Number of simulated observations')
ax[0].set_xlabel(r'RA offset (mas)')
ax[1].hist(3600*1000*(v['2011 OB60']['Dec_deg'] - insideccd['2011 OB60']['DEC_INTERP']), histtype='step', color='xkcd:cerulean blue')
ax[1].hist(3600*1000*(v['2010 TU149']['Dec_deg'] - insideccd['2010 TU149']['DEC_INTERP']), histtype='step',color='xkcd:scarlet')
ax[1].set_xlabel(r'Dec offset (mas)')
ax[2].hist(1000_0000*(v['2011 OB60']['trailedSourceMag'] - insideccd['2011 OB60']['mag']), histtype='step', color='xkcd:cerulean blue', linestyle='-' )
ax[2].hist(1000_0000*(v['2010 TU149']['trailedSourceMag'] - insideccd['2010 TU149']['mag']), histtype='step', color='xkcd:scarlet', linestyle='-')
ax[2].set_xlabel(r'Mag offset ($10^{-7}$ mag)')
ax[0].text(-1, 4.2, r'2010 TU$_{149}$', fontsize=16, color='xkcd:scarlet')
ax[0].text(0.1, 8, r'2011 OB$_{60}$', fontsize=16, color='xkcd:cerulean blue')
fig.set_size_inches(14,6)
pl.tight_layout()
pl.show()
And see the mean and standard deviation:
[19]:
pref = {'RA_deg' : 3600*1000, 'Dec_deg' : 3600*1000, 'trailedSourceMag' : 1e7} # note this means mas and 10^-7 mag!
for i in v:
for j, k in zip(['RA_deg', 'Dec_deg', 'trailedSourceMag'], ['RA_INTERP', 'DEC_INTERP', 'mag']):
print(i, j, np.mean(pref[j]*(v[i][j] - insideccd[i][k])), np.std(pref[j]*(v[i][j] - insideccd[i][k])))
2011 OB60 RA_deg 0.0008525950445648078 0.021086587833108447
2011 OB60 Dec_deg 8.785052685311712e-05 0.0075905242490396
2011 OB60 trailedSourceMag -1.1012382783448516 0.019616769792408627
2010 TU149 RA_deg 0.10147533791027286 1.0202124198587195
2010 TU149 Dec_deg 0.024679328869403605 0.4224011685568308
2010 TU149 trailedSourceMag -0.18841977854498915 0.16053158232152698