diff --git a/notebooks/sim_ephem/simulated_asteroid_ephem_rsp.ipynb b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp.ipynb new file mode 100644 index 0000000..b890d7f --- /dev/null +++ b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp.ipynb @@ -0,0 +1,635 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing DP0.3 ephemerides and coordinates\n", + "By Jamie Robinson" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the output from propagating the orbit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:30:46.535448Z", + "iopub.status.busy": "2024-07-05T15:30:46.535291Z", + "iopub.status.idle": "2024-07-05T15:30:49.002871Z", + "shell.execute_reply": "2024-07-05T15:30:49.002322Z", + "shell.execute_reply.started": "2024-07-05T15:30:46.535434Z" + } + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord, GCRS\n", + "from sbpy.data import Orbit\n", + "from astropy.time import Time\n", + "from astropy.coordinates import solar_system_ephemeris\n", + "from sbpy.photometry import HG\n", + "from astropy.table import QTable\n", + "\n", + "from lsst.rsp import get_tap_service" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:30:49.004023Z", + "iopub.status.busy": "2024-07-05T15:30:49.003641Z", + "iopub.status.idle": "2024-07-05T15:30:49.082021Z", + "shell.execute_reply": "2024-07-05T15:30:49.081508Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.004005Z" + } + }, + "outputs": [], + "source": [ + "service = get_tap_service(\"ssotap\")\n", + "assert service is not None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:30:49.082908Z", + "iopub.status.busy": "2024-07-05T15:30:49.082732Z", + "iopub.status.idle": "2024-07-05T15:30:49.085508Z", + "shell.execute_reply": "2024-07-05T15:30:49.085037Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.082894Z" + } + }, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"6098332225018\" # good test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:30:49.086275Z", + "iopub.status.busy": "2024-07-05T15:30:49.086098Z", + "iopub.status.idle": "2024-07-05T15:30:49.809760Z", + "shell.execute_reply": "2024-07-05T15:30:49.808044Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.086256Z" + } + }, + "outputs": [], + "source": [ + "query = \"\"\"\n", + "SELECT\n", + " *\n", + "FROM\n", + " dp03_catalogs_10yr.DiaSource as dia\n", + "INNER JOIN\n", + " dp03_catalogs_10yr.SSSource as sss\n", + "ON\n", + " dia.diaSourceId = sss.diaSourceId\n", + "WHERE\n", + " dia.ssObjectId={}\n", + "ORDER by dia.ssObjectId\n", + "\"\"\".format(\n", + " ssoid\n", + ")\n", + "\n", + "df = service.search(query).to_table().to_pandas()\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.810323Z", + "iopub.status.idle": "2024-07-05T15:30:49.810571Z", + "shell.execute_reply": "2024-07-05T15:30:49.810457Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.810448Z" + } + }, + "outputs": [], + "source": [ + "results = service.search(\"SELECT * FROM dp03_catalogs_10yr.MPCORB \" \"WHERE ssObjectId={}\".format(ssoid))\n", + "df_orb = results.to_table().to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.811051Z", + "iopub.status.idle": "2024-07-05T15:30:49.811253Z", + "shell.execute_reply": "2024-07-05T15:30:49.811163Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.811156Z" + } + }, + "outputs": [], + "source": [ + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.812059Z", + "iopub.status.idle": "2024-07-05T15:30:49.812280Z", + "shell.execute_reply": "2024-07-05T15:30:49.812190Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.812182Z" + } + }, + "outputs": [], + "source": [ + "# put obs in time order\n", + "df = df.sort_values(\"midPointMjdTai\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample the sparse observations by propagating orbit with sbpy & oorb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.812815Z", + "iopub.status.idle": "2024-07-05T15:30:49.813013Z", + "shell.execute_reply": "2024-07-05T15:30:49.812926Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.812918Z" + } + }, + "outputs": [], + "source": [ + "# check orbit\n", + "e = df_orb[\"e\"]\n", + "incl = df_orb[\"incl\"]\n", + "q = df_orb[\"q\"]\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.813642Z", + "iopub.status.idle": "2024-07-05T15:30:49.813847Z", + "shell.execute_reply": "2024-07-05T15:30:49.813758Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.813750Z" + } + }, + "outputs": [], + "source": [ + "df_orb[\"a\"] = a\n", + "df_orb[\"Q\"] = Q\n", + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.814629Z", + "iopub.status.idle": "2024-07-05T15:30:49.814844Z", + "shell.execute_reply": "2024-07-05T15:30:49.814751Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.814743Z" + } + }, + "outputs": [], + "source": [ + "# Calculate some extra orbital elements\n", + "df_orb[\"P\"] = df_orb[\"a\"] ** (3.0 / 2.0) # orbital period in years\n", + "df_orb[\"n\"] = 360.0 / (df_orb[\"P\"] * 365.25) # mean motion in deg/day\n", + "df_orb[\"M\"] = (\n", + " df_orb[\"n\"] * (df_orb[\"epoch\"] - df_orb[\"tperi\"])\n", + ") % 360 # angles must be in correct range otherwise sbpy/pyoorb freak out\n", + "df_orb[[\"a\", \"P\", \"n\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.815488Z", + "iopub.status.idle": "2024-07-05T15:30:49.815705Z", + "shell.execute_reply": "2024-07-05T15:30:49.815616Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.815608Z" + } + }, + "outputs": [], + "source": [ + "# rename columns for consistency\n", + "df_orb = df_orb.rename(columns={\"node\": \"Omega\", \"peri\": \"w\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.816210Z", + "iopub.status.idle": "2024-07-05T15:30:49.816409Z", + "shell.execute_reply": "2024-07-05T15:30:49.816322Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.816315Z" + } + }, + "outputs": [], + "source": [ + "df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.817161Z", + "iopub.status.idle": "2024-07-05T15:30:49.817360Z", + "shell.execute_reply": "2024-07-05T15:30:49.817271Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.817264Z" + } + }, + "outputs": [], + "source": [ + "# create an sbpy oorb object from dataframe via QTable\n", + "tab = QTable.from_pandas(\n", + " df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]],\n", + " units={\"a\": u.au, \"incl\": u.deg, \"Omega\": u.deg, \"w\": u.deg, \"M\": u.deg},\n", + ")\n", + "orbit = Orbit.from_table(tab)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.817935Z", + "iopub.status.idle": "2024-07-05T15:30:49.818130Z", + "shell.execute_reply": "2024-07-05T15:30:49.818043Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.818036Z" + } + }, + "outputs": [], + "source": [ + "# oorb requires certain extra fields\n", + "orbit[\"epoch\"] = Time(Time(df_orb[\"epoch\"], format=\"mjd\").jd, format=\"jd\")\n", + "orbit[\"targetname\"] = np.array(df_orb[\"ssObjectId\"]).astype(str)\n", + "orbit[\"H\"] = df_orb[\"mpcH\"] * u.mag\n", + "orbit[\"G\"] = df_orb[\"mpcG\"] * u.dimensionless_unscaled" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.818889Z", + "iopub.status.idle": "2024-07-05T15:30:49.819093Z", + "shell.execute_reply": "2024-07-05T15:30:49.819005Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.818998Z" + } + }, + "outputs": [], + "source": [ + "orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.819774Z", + "iopub.status.idle": "2024-07-05T15:30:49.819972Z", + "shell.execute_reply": "2024-07-05T15:30:49.819883Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.819876Z" + } + }, + "outputs": [], + "source": [ + "# define a set of JD times to propagate the orbital elements to\n", + "N = 1000\n", + "times = Time(\n", + " Time(np.linspace(np.amin(df[\"midPointMjdTai\"]), np.amax(df[\"midPointMjdTai\"]), N), format=\"mjd\").jd,\n", + " format=\"jd\",\n", + ")\n", + "times[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.820507Z", + "iopub.status.idle": "2024-07-05T15:30:49.820733Z", + "shell.execute_reply": "2024-07-05T15:30:49.820637Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.820629Z" + } + }, + "outputs": [], + "source": [ + "# create an empty dataframe to hold resampled observations\n", + "df_dense = pd.DataFrame()\n", + "df_dense[\"midPointMjdTai\"] = times.mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.821352Z", + "iopub.status.idle": "2024-07-05T15:30:49.821568Z", + "shell.execute_reply": "2024-07-05T15:30:49.821465Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.821457Z" + } + }, + "outputs": [], + "source": [ + "# propagate the orbit forward in time.\n", + "# probably a better way to do this but I can't get oo_propagate to work with a Time list right now\n", + "# see: https://github.com/NASA-Planetary-Science/sbpy/issues/341\n", + "\n", + "df_pos = pd.DataFrame() # empty dataframe to hold cartesian coordinates\n", + "\n", + "for i in range(len(times)):\n", + " print(i)\n", + " prop_elem = orbit.oo_propagate(times[i]) # propagate the orbit to the selected time step\n", + " del prop_elem.table[\n", + " \"orbtype\"\n", + " ] # orbtype is added as int, sbpy freaks out so delete the orbtype and then _to_oo works it out\n", + " print(\"propagate\")\n", + " statevec = prop_elem.oo_transform(\"CART\") # transform from orbital elements to cartesian\n", + " print(\"transform\")\n", + "\n", + " # append new cartesian coordinates to the dataframe\n", + " _df_statevec = statevec.table.to_pandas()\n", + " df_pos = pd.concat((df_pos, _df_statevec))\n", + "\n", + "df_pos.reset_index(drop=True, inplace=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.822161Z", + "iopub.status.idle": "2024-07-05T15:30:49.822363Z", + "shell.execute_reply": "2024-07-05T15:30:49.822271Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.822264Z" + } + }, + "outputs": [], + "source": [ + "df_pos" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Use astropy coordinates to transform between all the coordinate systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.822933Z", + "iopub.status.idle": "2024-07-05T15:30:49.823129Z", + "shell.execute_reply": "2024-07-05T15:30:49.823041Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.823034Z" + } + }, + "outputs": [], + "source": [ + "# define heliocentric cartesian coordinates\n", + "c_xyz_hel = SkyCoord(\n", + " x=np.array(df_pos[\"x\"]),\n", + " y=np.array(df_pos[\"y\"]),\n", + " z=np.array(df_pos[\"z\"]),\n", + " unit=\"AU\",\n", + " representation_type=\"cartesian\",\n", + " frame=\"heliocentrictrueecliptic\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.823966Z", + "iopub.status.idle": "2024-07-05T15:30:49.824165Z", + "shell.execute_reply": "2024-07-05T15:30:49.824078Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.824071Z" + } + }, + "outputs": [], + "source": [ + "# transform to heliocentric ecliptic coords\n", + "c_ecl_hel = c_xyz_hel.copy()\n", + "c_ecl_hel.representation_type = \"spherical\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.824676Z", + "iopub.status.idle": "2024-07-05T15:30:49.824877Z", + "shell.execute_reply": "2024-07-05T15:30:49.824783Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.824776Z" + } + }, + "outputs": [], + "source": [ + "# transform to geocentric equatorial coords (times required to calculate Earth position)\n", + "with solar_system_ephemeris.set(\"jpl\"):\n", + " c_eq_geo = c_xyz_hel.transform_to(GCRS(obstime=times))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.825499Z", + "iopub.status.idle": "2024-07-05T15:30:49.825709Z", + "shell.execute_reply": "2024-07-05T15:30:49.825623Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.825615Z" + } + }, + "outputs": [], + "source": [ + "# transform to geocentric cartesian coords\n", + "c_xyz_geo = c_eq_geo.copy()\n", + "c_xyz_geo.representation_type = \"cartesian\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.826270Z", + "iopub.status.idle": "2024-07-05T15:30:49.826473Z", + "shell.execute_reply": "2024-07-05T15:30:49.826380Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.826373Z" + } + }, + "outputs": [], + "source": [ + "# transform from geo equatorial (ra, dec) to geo ecliptic (lon, lat)\n", + "c_ecl_geo = c_eq_geo.transform_to(\"geocentrictrueecliptic\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.826974Z", + "iopub.status.idle": "2024-07-05T15:30:49.827178Z", + "shell.execute_reply": "2024-07-05T15:30:49.827091Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.827084Z" + } + }, + "outputs": [], + "source": [ + "# plot the propagated cartesian positions against the database values\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot1 = \"heliocentricX\"\n", + "y_plot2 = \"heliocentricY\"\n", + "y_plot3 = \"heliocentricZ\"\n", + "df_plot = df\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot3], label=y_plot3)\n", + "\n", + "ax1.plot(times.mjd, c_xyz_hel.x)\n", + "ax1.plot(times.mjd, c_xyz_hel.y)\n", + "ax1.plot(times.mjd, c_xyz_hel.z)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"distance\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are some deviations of x and y positions, probably due to slightly different reference frames and methods of propagating orbits.\n", + "\n", + "There is something wrong with database z positions, will be fixed soon!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:30:49.827679Z", + "iopub.status.idle": "2024-07-05T15:30:49.827879Z", + "shell.execute_reply": "2024-07-05T15:30:49.827788Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.827780Z" + } + }, + "outputs": [], + "source": [ + "# the ecliptic coordinates look good!\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot1 = \"eclipticLambda\"\n", + "y_plot2 = \"eclipticBeta\"\n", + "df_plot = df\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "\n", + "ax1.plot(times.mjd, c_ecl_geo.lon.degree)\n", + "ax1.plot(times.mjd, c_ecl_geo.lat.degree)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"angle\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "lsst-sbpy", + "language": "python", + "name": "lsst-sbpy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb new file mode 100644 index 0000000..fb45ddb --- /dev/null +++ b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb @@ -0,0 +1,910 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing DP0.3 ephemerides and coordinates\n", + "By Jamie Robinson" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the JPL Horizons values." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:34.945521Z", + "iopub.status.busy": "2024-07-05T15:33:34.945104Z", + "iopub.status.idle": "2024-07-05T15:33:39.271163Z", + "shell.execute_reply": "2024-07-05T15:33:39.270462Z", + "shell.execute_reply.started": "2024-07-05T15:33:34.945496Z" + } + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord, GCRS\n", + "from sbpy.data import Orbit\n", + "from astropy.time import Time\n", + "from astropy.coordinates import solar_system_ephemeris\n", + "from sbpy.photometry import HG\n", + "from astropy.table import QTable\n", + "\n", + "from lsst.rsp import get_tap_service\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:39.272663Z", + "iopub.status.busy": "2024-07-05T15:33:39.272094Z", + "iopub.status.idle": "2024-07-05T15:33:39.333026Z", + "shell.execute_reply": "2024-07-05T15:33:39.332159Z", + "shell.execute_reply.started": "2024-07-05T15:33:39.272639Z" + } + }, + "outputs": [], + "source": [ + "service = get_tap_service(\"ssotap\")\n", + "assert service is not None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:39.334212Z", + "iopub.status.busy": "2024-07-05T15:33:39.333990Z", + "iopub.status.idle": "2024-07-05T15:33:39.337089Z", + "shell.execute_reply": "2024-07-05T15:33:39.336551Z", + "shell.execute_reply.started": "2024-07-05T15:33:39.334195Z" + } + }, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "# ssoid = \"6098332225018\" # good test object\n", + "# ssoid = \"8268570668335894776\" # NEO\n", + "ssoid = \"-3649348068486548794\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:39.338054Z", + "iopub.status.busy": "2024-07-05T15:33:39.337859Z", + "iopub.status.idle": "2024-07-05T15:33:40.267531Z", + "shell.execute_reply": "2024-07-05T15:33:40.264933Z", + "shell.execute_reply.started": "2024-07-05T15:33:39.338038Z" + } + }, + "outputs": [], + "source": [ + "query = \"\"\"\n", + "SELECT\n", + " *\n", + "FROM\n", + " dp03_catalogs_10yr.DiaSource as dia\n", + "INNER JOIN\n", + " dp03_catalogs_10yr.SSSource as sss\n", + "ON\n", + " dia.diaSourceId = sss.diaSourceId\n", + "WHERE\n", + " dia.ssObjectId={}\n", + "ORDER by dia.ssObjectId\n", + "\"\"\".format(\n", + " ssoid\n", + ")\n", + "\n", + "df_obs = service.search(query).to_table().to_pandas()\n", + "df_obs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.268400Z", + "iopub.status.idle": "2024-07-05T15:33:40.268689Z", + "shell.execute_reply": "2024-07-05T15:33:40.268574Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.268562Z" + } + }, + "outputs": [], + "source": [ + "# calculate elongation angle\n", + "R = df_obs[\"heliocentricDist\"]\n", + "Delta = df_obs[\"topocentricDist\"]\n", + "alpha = np.radians(df_obs[\"phaseAngle\"])\n", + "\n", + "R_E = np.sqrt((R * R) + (Delta * Delta) - (2.0 * R * Delta * np.cos(alpha)))\n", + "df_obs[\"elong\"] = np.degrees(np.arccos(((R_E * R_E) + (Delta * Delta) - (R * R)) / (2.0 * R_E * Delta)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.269564Z", + "iopub.status.idle": "2024-07-05T15:33:40.269819Z", + "shell.execute_reply": "2024-07-05T15:33:40.269697Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.269688Z" + } + }, + "outputs": [], + "source": [ + "# put obs in time order\n", + "df_obs = df_obs.sort_values(\"midPointMjdTai\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.270771Z", + "iopub.status.idle": "2024-07-05T15:33:40.271049Z", + "shell.execute_reply": "2024-07-05T15:33:40.270945Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.270933Z" + } + }, + "outputs": [], + "source": [ + "# get orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.271857Z", + "iopub.status.idle": "2024-07-05T15:33:40.272093Z", + "shell.execute_reply": "2024-07-05T15:33:40.271983Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.271974Z" + } + }, + "outputs": [], + "source": [ + "results = service.search(\"SELECT * FROM dp03_catalogs_10yr.MPCORB \" \"WHERE ssObjectId={}\".format(ssoid))\n", + "df_orb = results.to_table().to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.273422Z", + "iopub.status.idle": "2024-07-05T15:33:40.273664Z", + "shell.execute_reply": "2024-07-05T15:33:40.273559Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.273548Z" + } + }, + "outputs": [], + "source": [ + "df_orb" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample the sparse observations by propagating orbit with sbpy & oorb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.274352Z", + "iopub.status.idle": "2024-07-05T15:33:40.274586Z", + "shell.execute_reply": "2024-07-05T15:33:40.274472Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.274463Z" + } + }, + "outputs": [], + "source": [ + "# check orbit\n", + "e = df_orb[\"e\"]\n", + "incl = df_orb[\"incl\"]\n", + "q = df_orb[\"q\"]\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.275180Z", + "iopub.status.idle": "2024-07-05T15:33:40.275404Z", + "shell.execute_reply": "2024-07-05T15:33:40.275301Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.275292Z" + } + }, + "outputs": [], + "source": [ + "df_orb[\"a\"] = a\n", + "df_orb[\"Q\"] = Q\n", + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.276437Z", + "iopub.status.idle": "2024-07-05T15:33:40.276674Z", + "shell.execute_reply": "2024-07-05T15:33:40.276573Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.276563Z" + } + }, + "outputs": [], + "source": [ + "# Calculate some extra orbital elements\n", + "df_orb[\"P\"] = df_orb[\"a\"] ** (3.0 / 2.0) # orbital period in years\n", + "df_orb[\"n\"] = 360.0 / (df_orb[\"P\"] * 365.25) # mean motion in deg/day\n", + "df_orb[\"M\"] = (\n", + " df_orb[\"n\"] * (df_orb[\"epoch\"] - df_orb[\"tperi\"])\n", + ") % 360 # angles must be in correct range otherwise sbpy/pyoorb freak out\n", + "df_orb[[\"a\", \"P\", \"n\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.277508Z", + "iopub.status.idle": "2024-07-05T15:33:40.277727Z", + "shell.execute_reply": "2024-07-05T15:33:40.277634Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.277626Z" + } + }, + "outputs": [], + "source": [ + "# rename columns for consistency\n", + "df_orb = df_orb.rename(columns={\"node\": \"Omega\", \"peri\": \"w\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.278432Z", + "iopub.status.idle": "2024-07-05T15:33:40.278645Z", + "shell.execute_reply": "2024-07-05T15:33:40.278549Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.278541Z" + } + }, + "outputs": [], + "source": [ + "df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.289637Z", + "iopub.status.idle": "2024-07-05T15:33:40.290438Z", + "shell.execute_reply": "2024-07-05T15:33:40.290280Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.290263Z" + } + }, + "outputs": [], + "source": [ + "# # create an sbpy oorb object from dataframe via QTable\n", + "# tab = QTable.from_pandas(\n", + "# df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]],\n", + "# units={\"a\": u.au, \"incl\": u.deg, \"Omega\": u.deg, \"w\": u.deg, \"M\": u.deg},\n", + "# )\n", + "# orbit = Orbit.from_table(tab)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.291224Z", + "iopub.status.idle": "2024-07-05T15:33:40.291730Z", + "shell.execute_reply": "2024-07-05T15:33:40.291609Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.291597Z" + } + }, + "outputs": [], + "source": [ + "# # oorb requires certain extra fields\n", + "# orbit[\"epoch\"] = Time(Time(df_orb[\"epoch\"], format=\"mjd\").jd, format=\"jd\")\n", + "# orbit[\"targetname\"] = np.array(df_orb[\"ssObjectId\"]).astype(str)\n", + "# orbit[\"H\"] = df_orb[\"mpcH\"] * u.mag\n", + "# orbit[\"G\"] = df_orb[\"mpcG\"] * u.dimensionless_unscaled" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.292524Z", + "iopub.status.idle": "2024-07-05T15:33:40.292763Z", + "shell.execute_reply": "2024-07-05T15:33:40.292661Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.292652Z" + } + }, + "outputs": [], + "source": [ + "# orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.293660Z", + "iopub.status.idle": "2024-07-05T15:33:40.293932Z", + "shell.execute_reply": "2024-07-05T15:33:40.293825Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.293812Z" + } + }, + "outputs": [], + "source": [ + "# # define a set of JD times to propagate the orbital elements to\n", + "# N = 1000\n", + "# times = Time(\n", + "# Time(np.linspace(np.amin(df[\"midPointMjdTai\"]), np.amax(df[\"midPointMjdTai\"]), N), format=\"mjd\").jd,\n", + "# format=\"jd\",\n", + "# )\n", + "# times[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.294706Z", + "iopub.status.idle": "2024-07-05T15:33:40.294965Z", + "shell.execute_reply": "2024-07-05T15:33:40.294855Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.294845Z" + } + }, + "outputs": [], + "source": [ + "# # create an empty dataframe to hold resampled observations\n", + "# df_dense = pd.DataFrame()\n", + "# df_dense[\"midPointMjdTai\"] = times.mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.296032Z", + "iopub.status.idle": "2024-07-05T15:33:40.296273Z", + "shell.execute_reply": "2024-07-05T15:33:40.296163Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.296154Z" + } + }, + "outputs": [], + "source": [ + "# # propagate the orbit forward in time.\n", + "# # probably a better way to do this but I can't get oo_propagate to work with a Time list right now\n", + "# # see: https://github.com/NASA-Planetary-Science/sbpy/issues/341\n", + "\n", + "# df_pos = pd.DataFrame() # empty dataframe to hold cartesian coordinates\n", + "\n", + "# for i in range(len(times)):\n", + "# print(i)\n", + "# prop_elem = orbit.oo_propagate(times[i]) # propagate the orbit to the selected time step\n", + "# del prop_elem.table[\n", + "# \"orbtype\"\n", + "# ] # orbtype is added as int, sbpy freaks out so delete the orbtype and then _to_oo works it out\n", + "# print(\"propagate\")\n", + "# statevec = prop_elem.oo_transform(\"CART\") # transform from orbital elements to cartesian\n", + "# print(\"transform\")\n", + "\n", + "# # append new cartesian coordinates to the dataframe\n", + "# _df_statevec = statevec.table.to_pandas()\n", + "# df_pos = pd.concat((df_pos, _df_statevec))\n", + "\n", + "# df_pos.reset_index(drop=True, inplace=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.296898Z", + "iopub.status.idle": "2024-07-05T15:33:40.297138Z", + "shell.execute_reply": "2024-07-05T15:33:40.297029Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.297018Z" + } + }, + "outputs": [], + "source": [ + "# df_pos" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Use astropy coordinates to transform between all the coordinate systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.298026Z", + "iopub.status.idle": "2024-07-05T15:33:40.298250Z", + "shell.execute_reply": "2024-07-05T15:33:40.298152Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.298143Z" + } + }, + "outputs": [], + "source": [ + "# # define heliocentric cartesian coordinates\n", + "# c_xyz_hel = SkyCoord(\n", + "# x=np.array(df_pos[\"x\"]),\n", + "# y=np.array(df_pos[\"y\"]),\n", + "# z=np.array(df_pos[\"z\"]),\n", + "# unit=\"AU\",\n", + "# representation_type=\"cartesian\",\n", + "# frame=\"heliocentrictrueecliptic\",\n", + "# )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.298942Z", + "iopub.status.idle": "2024-07-05T15:33:40.299172Z", + "shell.execute_reply": "2024-07-05T15:33:40.299073Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.299064Z" + } + }, + "outputs": [], + "source": [ + "# # transform to heliocentric ecliptic coords\n", + "# c_ecl_hel = c_xyz_hel.copy()\n", + "# c_ecl_hel.representation_type = \"spherical\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.299978Z", + "iopub.status.idle": "2024-07-05T15:33:40.300217Z", + "shell.execute_reply": "2024-07-05T15:33:40.300117Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.300108Z" + } + }, + "outputs": [], + "source": [ + "# # transform to geocentric equatorial coords (times required to calculate Earth position)\n", + "# with solar_system_ephemeris.set(\"jpl\"):\n", + "# c_eq_geo = c_xyz_hel.transform_to(GCRS(obstime=times))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.301111Z", + "iopub.status.idle": "2024-07-05T15:33:40.301359Z", + "shell.execute_reply": "2024-07-05T15:33:40.301243Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.301233Z" + } + }, + "outputs": [], + "source": [ + "# # transform to geocentric cartesian coords\n", + "# c_xyz_geo = c_eq_geo.copy()\n", + "# c_xyz_geo.representation_type = \"cartesian\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.301927Z", + "iopub.status.idle": "2024-07-05T15:33:40.302147Z", + "shell.execute_reply": "2024-07-05T15:33:40.302042Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.302034Z" + } + }, + "outputs": [], + "source": [ + "# # transform from geo equatorial (ra, dec) to geo ecliptic (lon, lat)\n", + "# c_ecl_geo = c_eq_geo.transform_to(\"geocentrictrueecliptic\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.302981Z", + "iopub.status.idle": "2024-07-05T15:33:40.303210Z", + "shell.execute_reply": "2024-07-05T15:33:40.303111Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.303102Z" + } + }, + "outputs": [], + "source": [ + "# # plot the propagated cartesian positions against the database values\n", + "\n", + "# x_plot = \"midPointMjdTai\"\n", + "# y_plot1 = \"heliocentricX\"\n", + "# y_plot2 = \"heliocentricY\"\n", + "# y_plot3 = \"heliocentricZ\"\n", + "# df_plot = df\n", + "\n", + "# fig = plt.figure()\n", + "# gs = gridspec.GridSpec(1, 1)\n", + "# ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot3], label=y_plot3)\n", + "\n", + "# ax1.plot(times.mjd, c_xyz_hel.x)\n", + "# ax1.plot(times.mjd, c_xyz_hel.y)\n", + "# ax1.plot(times.mjd, c_xyz_hel.z)\n", + "\n", + "# ax1.set_xlabel(x_plot)\n", + "# ax1.set_ylabel(\"distance\")\n", + "# ax1.legend()\n", + "\n", + "# plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are some deviations of x and y positions, probably due to slightly different reference frames and methods of propagating orbits.\n", + "\n", + "There is something wrong with database z positions, will be fixed soon!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.303807Z", + "iopub.status.idle": "2024-07-05T15:33:40.304051Z", + "shell.execute_reply": "2024-07-05T15:33:40.303949Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.303940Z" + } + }, + "outputs": [], + "source": [ + "# # the ecliptic coordinates look good!\n", + "\n", + "# x_plot = \"midPointMjdTai\"\n", + "# y_plot1 = \"eclipticLambda\"\n", + "# y_plot2 = \"eclipticBeta\"\n", + "# df_plot = df\n", + "\n", + "# fig = plt.figure()\n", + "# gs = gridspec.GridSpec(1, 1)\n", + "# ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "\n", + "# ax1.plot(times.mjd, c_ecl_geo.lon.degree)\n", + "# ax1.plot(times.mjd, c_ecl_geo.lat.degree)\n", + "\n", + "# ax1.set_xlabel(x_plot)\n", + "# ax1.set_ylabel(\"angle\")\n", + "# ax1.legend()\n", + "\n", + "# plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.304969Z", + "iopub.status.idle": "2024-07-05T15:33:40.305201Z", + "shell.execute_reply": "2024-07-05T15:33:40.305099Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.305090Z" + } + }, + "outputs": [], + "source": [ + "# query JPL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.305853Z", + "iopub.status.idle": "2024-07-05T15:33:40.306082Z", + "shell.execute_reply": "2024-07-05T15:33:40.305973Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.305966Z" + } + }, + "outputs": [], + "source": [ + "target = df_orb.iloc[0][\"fullDesignation\"].split(\"2011 \")[-1]\n", + "print(target)\n", + "\n", + "site = \"X05\" # Roques de los Muchachos\n", + "times = {\"start\": \"2023-10-01 04:00\", \"stop\": \"2033-10-01 04:00\", \"step\": \"1day\"} # dates to query\n", + "obj = Horizons(id=target, location=site, epochs=times)\n", + "eph = obj.ephemerides()\n", + "df_eph = eph.to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.307034Z", + "iopub.status.idle": "2024-07-05T15:33:40.307249Z", + "shell.execute_reply": "2024-07-05T15:33:40.307156Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.307148Z" + } + }, + "outputs": [], + "source": [ + "obj = Horizons(id=target, epochs=times, location=\"399\")\n", + "vec = obj.vectors()\n", + "df_vec_earth = vec.to_pandas()\n", + "\n", + "obj = Horizons(id=target, epochs=times, location=\"@10\")\n", + "vec = obj.vectors()\n", + "df_vec_sun = vec.to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.308110Z", + "iopub.status.idle": "2024-07-05T15:33:40.308328Z", + "shell.execute_reply": "2024-07-05T15:33:40.308227Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.308220Z" + } + }, + "outputs": [], + "source": [ + "df_vec_earth = df_vec_earth.rename({\"x\": \"topocentricX\", \"y\": \"topocentricY\", \"z\": \"topocentricZ\"}, axis=1)\n", + "df_vec_sun = df_vec_sun.rename({\"x\": \"heliocentricX\", \"y\": \"heliocentricY\", \"z\": \"heliocentricZ\"}, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.309065Z", + "iopub.status.idle": "2024-07-05T15:33:40.309290Z", + "shell.execute_reply": "2024-07-05T15:33:40.309190Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.309182Z" + } + }, + "outputs": [], + "source": [ + "df_vec = df_vec_earth[\n", + " [\"targetname\", \"datetime_jd\", \"datetime_str\", \"topocentricX\", \"topocentricY\", \"topocentricZ\"]\n", + "].merge(\n", + " df_vec_sun[\n", + " [\"targetname\", \"datetime_jd\", \"datetime_str\", \"heliocentricX\", \"heliocentricY\", \"heliocentricZ\"]\n", + " ],\n", + " on=[\"targetname\", \"datetime_jd\", \"datetime_str\"],\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.310032Z", + "iopub.status.idle": "2024-07-05T15:33:40.310264Z", + "shell.execute_reply": "2024-07-05T15:33:40.310165Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.310156Z" + } + }, + "outputs": [], + "source": [ + "df_eph[\"datetime_mjd\"] = Time(df_eph[\"datetime_jd\"], format=\"jd\").mjd\n", + "df_vec[\"datetime_mjd\"] = Time(df_vec[\"datetime_jd\"], format=\"jd\").mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.311155Z", + "iopub.status.idle": "2024-07-05T15:33:40.311405Z", + "shell.execute_reply": "2024-07-05T15:33:40.311286Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.311277Z" + } + }, + "outputs": [], + "source": [ + "night_mask = (~np.isin(df_eph[\"solar_presence\"], [\"*\", \"N\", \"C\"])) & (df_eph[\"EL\"] > 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.312011Z", + "iopub.status.idle": "2024-07-05T15:33:40.312264Z", + "shell.execute_reply": "2024-07-05T15:33:40.312149Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.312141Z" + } + }, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"elong\"\n", + "df_plot = df_obs\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_plot[x_plot], df_plot[y_plot], label=\"DP0.3\")\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"elong\"], c=\"r\", label=\"JPL\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.313021Z", + "iopub.status.idle": "2024-07-05T15:33:40.313245Z", + "shell.execute_reply": "2024-07-05T15:33:40.313148Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.313139Z" + } + }, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"phaseAngle\"\n", + "df_plot = df_obs\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_plot[x_plot], df_plot[y_plot], label=\"DP0.3\")\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"alpha\"], c=\"r\", label=\"JPL\")\n", + "ax1.scatter(df_eph[night_mask][\"datetime_mjd\"], df_eph[night_mask][\"alpha\"], edgecolor=\"C1\", facecolor=\"none\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.314008Z", + "iopub.status.idle": "2024-07-05T15:33:40.314240Z", + "shell.execute_reply": "2024-07-05T15:33:40.314140Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.314130Z" + } + }, + "outputs": [], + "source": [ + "x_plot = \"heliocentricX\"\n", + "y_plot = \"heliocentricY\"\n", + "c_plot = \"midPointMjdTai\"\n", + "df_plot = df_obs\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "s1 = ax1.scatter(df_plot[x_plot], df_plot[y_plot], c=df_plot[c_plot])\n", + "cbar1 = plt.colorbar(s1)\n", + "\n", + "ax1.scatter(0, 0, marker=\"+\", c=\"k\")\n", + "circle1 = plt.Circle((0, 0), 1.0, edgecolor=\"r\", facecolor=\"none\", label=\"1au\")\n", + "ax1.add_patch(circle1)\n", + "\n", + "mask = df_plot[\"phaseAngle\"] > 120\n", + "_df_plot = df_plot[mask]\n", + "ax1.scatter(_df_plot[x_plot], _df_plot[y_plot], facecolor=\"none\", edgecolor=\"r\")\n", + "\n", + "ax1.set_aspect(\"equal\")\n", + "ax1.legend()\n", + "\n", + "cbar1.set_label(c_plot)\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "lsst-sbpy", + "language": "python", + "name": "lsst-sbpy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}