diff --git a/notebooks/colour_functions_testing.ipynb b/notebooks/colour_functions_testing.ipynb new file mode 100644 index 0000000..b1a2947 --- /dev/null +++ b/notebooks/colour_functions_testing.ipynb @@ -0,0 +1,2057 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "d1ef57c2", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61f50e0d", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82537a61", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.dataclasses.AdlerData import AdlerData\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "from adler.science.Colour import col_obs_ref\n", + "from adler.utilities.plotting_utilities import plot_errorbar\n", + "from adler.utilities.science_utilities import apparition_gap_finder, get_df_obs_filt\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import astropy.units as u\n", + "from scipy.interpolate import splrep, BSpline, CubicSpline\n", + "from astropy.modeling.fitting import SLSQPLSQFitter\n", + "from astropy.time import Time\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9f6dcda", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "# ssoid = \"8268570668335894776\" # NEO\n", + "ssoid = \"6098332225018\" # good MBA test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6f705c1", + "metadata": {}, + "outputs": [], + "source": [ + "# fname = \"../tests/data/testing_database.db\"\n", + "fname = \"/Users/jrobinson/lsst-adler/notebooks/gen_test_data/adler_demo_testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b756aac", + "metadata": {}, + "outputs": [], + "source": [ + "# known issue with mpcDesignation\n", + "# https://dp0-3.lsst.io/data-products-dp0-3/data-simulation-dp0-3.html#known-issues\n", + "planetoid.MPCORB.mpcDesignation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73ea4e74", + "metadata": {}, + "outputs": [], + "source": [ + "planetoid.MPCORB.fullDesignation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "644381d2", + "metadata": {}, + "outputs": [], + "source": [ + "mpc_des = planetoid.MPCORB.fullDesignation\n", + "mpc_des = mpc_des.split(\"2011 \")[-1]\n", + "planetoid.MPCORB.fullDesignation = mpc_des\n", + "mpc_des" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21e0b93e", + "metadata": {}, + "outputs": [], + "source": [ + "# check orbit\n", + "e = planetoid.MPCORB.e\n", + "incl = planetoid.MPCORB.incl\n", + "q = planetoid.MPCORB.q\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8939e778", + "metadata": {}, + "outputs": [], + "source": [ + "dir(planetoid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c4defd8e", + "metadata": {}, + "outputs": [], + "source": [ + "planetoid.SSObject" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee4064e7", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "adler_cols = AdlerData(ssoid, planetoid.filter_list)\n", + "\n", + "for filt in [\"g\", \"r\"]:\n", + " # get observations and define a phase curve model to fit\n", + "\n", + " sso = planetoid.SSObject_in_filter(filt)\n", + " obs = planetoid.observations_in_filter(filt)\n", + "\n", + " H = sso.H\n", + " G12 = sso.G12\n", + " pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name=\"HG12_Pen16\")\n", + "\n", + " # H = sso.H\n", + " # S = 0.2\n", + " # pc = PhaseCurve(H=H * u.mag, phase_parameter_1=S* (u.mag/u.deg), model_name=\"LinearPhaseFunc\")\n", + "\n", + " # H = sso.H\n", + " # G1 = 0.2\n", + " # G2 = 0.2\n", + " # pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G1, phase_parameter_2=G2, model_name=\"HG1G2\")\n", + "\n", + " alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg\n", + " # red_mag = pc.ReducedMag(alpha)\n", + " # pc.model_function.G12.fixed = True\n", + " # pc.FixParam(\"S\")\n", + " # pc.FixParam(\"phase_parameter_1\")\n", + " # pc.model_function.H.fixed = True\n", + " # pc.model_function.G1.fixed = True\n", + " # pc.model_function.G2.fixed = True\n", + "\n", + " pc_fit = pc.FitModel(\n", + " np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n", + " np.array(getattr(obs, \"reduced_mag\")) * u.mag,\n", + " # np.array(getattr(obs, \"magErr\")) * u.mag,\n", + " # fitter = SLSQPLSQFitter()\n", + " )\n", + " pc = pc.InitModelSbpy(pc_fit)\n", + " red_mag = pc.ReducedMag(alpha)\n", + "\n", + " adler_cols.populate_phase_parameters(filt, **pc.ReturnModelDict())\n", + "\n", + " # add this phase curve to the figure\n", + " fig = plot_errorbar(planetoid, filt_list=[filt], fig=fig)\n", + " ax1 = fig.axes[0]\n", + " ax1.plot(\n", + " alpha.value,\n", + " pc.ReducedMag(alpha).value,\n", + " label=\"{}: H={:.2f}, G12={:.2f}\".format(filt, pc.H, pc.phase_parameter_1),\n", + " )\n", + " ax1.legend()\n", + "\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0df3c51", + "metadata": {}, + "outputs": [], + "source": [ + "adler_cols.get_phase_parameters_in_filter(\"r\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c46fb9eb", + "metadata": {}, + "outputs": [], + "source": [ + "adler_cols.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0977c37e", + "metadata": {}, + "outputs": [], + "source": [ + "planetoid.SSObject_in_filter(\"g\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "637fce52", + "metadata": {}, + "outputs": [], + "source": [ + "# # add this phase curve to the figure\n", + "# fig = plot_errorbar(planetoid, filt_list=[filt])\n", + "# ax1 = fig.axes[0]\n", + "# ax1.plot(alpha.value, pc.ReducedMag(alpha).value, label=\"{} {}\".format(filt, pc.model_name))\n", + "# ax1.legend()\n", + "# plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2306ca40", + "metadata": {}, + "outputs": [], + "source": [ + "# calculate data - model residuals\n", + "res = pc.ModelResiduals(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value\n", + "res" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67450145", + "metadata": {}, + "outputs": [], + "source": [ + "# calculate the absolute mag\n", + "abs_mag = pc.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value\n", + "abs_mag" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff80d360", + "metadata": {}, + "outputs": [], + "source": [ + "# x_plot = \"phaseAngle\"\n", + "x_plot = \"midPointMjdTai\"\n", + "x = getattr(obs, x_plot)\n", + "# y = res\n", + "y = abs_mag\n", + "yerr = obs.magErr\n", + "\n", + "sort_mask = np.argsort(x)\n", + "x = x[sort_mask]\n", + "y = y[sort_mask]\n", + "\n", + "tck1 = splrep(x, y, s=0)\n", + "tck2 = splrep(x, y, s=len(x))\n", + "tck3 = splrep(x, y, w=1.0 / yerr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "305adc60", + "metadata": {}, + "outputs": [], + "source": [ + "tck3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40564ff0", + "metadata": {}, + "outputs": [], + "source": [ + "len(tck1[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9919ee9", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the data - model residuals\n", + "# x_plot = \"phaseAngle\"\n", + "x_plot = \"midPointMjdTai\"\n", + "\n", + "x = getattr(obs, x_plot)\n", + "yerr = obs.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# ax1.errorbar(x, res, yerr, fmt=\"o\")\n", + "ax1.errorbar(x, abs_mag, yerr, fmt=\"o\")\n", + "\n", + "# ax1.axhline(0, c=\"k\")\n", + "ax1.axhline(pc.H.value, c=\"k\")\n", + "\n", + "xnew = np.linspace(np.amin(x), np.amax(x), 1000)\n", + "# ax1.plot(xnew, BSpline(*tck1)(xnew), '-', label='s=0')\n", + "ax1.plot(xnew, BSpline(*tck2)(xnew), \"-\", label=f\"s={len(x)}\")\n", + "ax1.plot(xnew, BSpline(*tck3)(xnew), \"-\", label=\"weighted smooth\")\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"data - model\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1bc11ca3", + "metadata": {}, + "outputs": [], + "source": [ + "# # use a running mean?\n", + "\n", + "# def movingaverage(interval, window_size):\n", + "# window = np.ones(int(window_size))/float(window_size)\n", + "# return np.convolve(interval, window, 'same')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a02fb173", + "metadata": {}, + "outputs": [], + "source": [ + "# # plot the data - model residuals\n", + "# x_plot = \"phaseAngle\"\n", + "\n", + "# x = getattr(obs, x_plot)\n", + "# yerr = obs.magErr\n", + "\n", + "# fig = plt.figure()\n", + "# gs = gridspec.GridSpec(1, 1)\n", + "# ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# # ax1.errorbar(x, res, yerr, fmt=\"o\")\n", + "# ax1.errorbar(x, abs_mag, yerr, fmt=\"o\")\n", + "\n", + "# # ax1.axhline(0, c=\"k\")\n", + "# ax1.axhline(pc.H.value, c=\"k\")\n", + "\n", + "# xnew = np.linspace(np.amin(x),np.amax(x),1000)\n", + "# # ax1.plot(xnew, BSpline(*tck1)(xnew), '-', label='s=0')\n", + "# ax1.plot(xnew, BSpline(*tck2)(xnew), '-', label=f's={len(x)}')\n", + "# ax1.plot(xnew, BSpline(*tck3)(xnew), '-', label='weighted smooth')\n", + "# ax1.scatter(x, BSpline(*tck3)(x), c = \"C2\")\n", + "\n", + "\n", + "# # x_plot = \"phaseAngle\"\n", + "# # x = getattr(obs, x_plot)\n", + "# # # y = res\n", + "# # y = abs_mag\n", + "# # sort_mask = np.argsort(x)\n", + "# # x = x[sort_mask]\n", + "# # y = y[sort_mask]\n", + "# # y_av = movingaverage(y, 5)\n", + "# # ax1.plot(x, y_av,label = \"moving avg\")\n", + "\n", + "\n", + "# ax1.invert_yaxis()\n", + "# ax1.set_xlabel(x_plot)\n", + "# ax1.set_ylabel(\"data - model\")\n", + "# ax1.legend()\n", + "\n", + "# plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6efb822b", + "metadata": {}, + "outputs": [], + "source": [ + "# generate the spline points at all points of a given filter set\n", + "# find the difference at these points\n", + "# do we consider differences where there is no data?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb5deb4b", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the data - model residuals\n", + "# x_plot = \"phaseAngle\"\n", + "x_plot = \"midPointMjdTai\"\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# # get min,max phase angle\n", + "# phase_min = []\n", + "# phase_max = []\n", + "# for filt in [\"g\",\"r\"]:\n", + "# obs = planetoid.observations_in_filter(filt)\n", + "# phase_min.append(np.amin(obs.phaseAngle))\n", + "# phase_max.append(np.amax(obs.phaseAngle))\n", + "# # phase_min = np.amin(np.array(phase_min))\n", + "# # phase_max = np.amax(np.array(phase_max))\n", + "\n", + "# # min and max that encloses both data sets\n", + "# phase_min = np.amax(np.array(phase_min))\n", + "# phase_max = np.amin(np.array(phase_max))\n", + "# print(phase_min,phase_max)\n", + "\n", + "# get min,max time\n", + "time_min = []\n", + "time_max = []\n", + "for filt in [\"g\", \"r\"]:\n", + " obs = planetoid.observations_in_filter(filt)\n", + " time_min.append(np.amin(obs.midPointMjdTai))\n", + " time_max.append(np.amax(obs.midPointMjdTai))\n", + "# min and max that encloses both data sets\n", + "time_min = np.amax(np.array(time_min))\n", + "time_max = np.amin(np.array(time_max))\n", + "print(time_min, time_max)\n", + "\n", + "splines = []\n", + "xdata = np.array([])\n", + "for filt in [\"g\", \"r\"]:\n", + " # get observations in filter\n", + " obs = planetoid.observations_in_filter(filt)\n", + " # get the stored AdlerData parameters for this filter\n", + " ad = adler_cols.get_phase_parameters_in_filter(filt, \"HG12_Pen16\")\n", + " # make the PhaseCurve object from AdlerData\n", + " pc = PhaseCurve().InitModelDict(ad.__dict__)\n", + " print(pc.__dict__)\n", + "\n", + " # calculate the model absolute magnitude\n", + " abs_mag = pc.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value\n", + "\n", + " x = getattr(obs, x_plot)\n", + " xdata = np.concatenate([xdata, x])\n", + " yerr = obs.magErr\n", + "\n", + " # ax1.errorbar(x, res, yerr, fmt=\"o\")\n", + " ax1.errorbar(x, abs_mag, yerr, fmt=\"o\", label=filt)\n", + "\n", + " # ax1.axhline(0, c=\"k\")\n", + " ax1.axhline(pc.H.value, c=\"k\", label=\"H_{} = {:.2f} mag\".format(filt, pc.H.value))\n", + "\n", + " y = abs_mag\n", + " sort_mask = np.argsort(x)\n", + " x = x[sort_mask]\n", + " y = y[sort_mask]\n", + " tck = splrep(x, y, w=1.0 / yerr)\n", + " splines.append(tck)\n", + " # xnew = np.linspace(np.amin(x),np.amax(x),1000)\n", + " # xnew = np.linspace(phase_min,phase_max,1000)\n", + " xnew = np.linspace(time_min, time_max, 1000)\n", + " ax1.plot(xnew, BSpline(*tck)(xnew), \"-\", label=\"weighted smooth\")\n", + "\n", + "# break\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"abs_mag\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88c61235", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# xnew = np.linspace(phase_min,phase_max,1000)\n", + "xnew = np.linspace(time_min, time_max, 1000)\n", + "spline_diff = BSpline(*splines[0])(xnew) - BSpline(*splines[1])(xnew)\n", + "ax1.plot(xnew, spline_diff, \"-\", label=\"weighted smooth\", c=\"k\", zorder=1)\n", + "\n", + "phase = np.array([])\n", + "spline_diff = np.array([])\n", + "for i, filt in enumerate([\"g\", \"r\"]):\n", + " obs = planetoid.observations_in_filter(filt)\n", + " # x = getattr(obs, \"phaseAngle\")\n", + " x = getattr(obs, \"midPointMjdTai\")\n", + " yerr = obs.magErr\n", + " # xdata_mask = (x>=phase_min) & (x<=phase_max)\n", + " xdata_mask = (x >= time_min) & (x <= time_max)\n", + " _spline_diff = BSpline(*splines[0])(x[xdata_mask]) - BSpline(*splines[1])(x[xdata_mask])\n", + " # ax1.scatter(x[xdata_mask], _spline_diff,\n", + " # label = filt, facecolor = \"none\", edgecolor = \"C{}\".format(i))\n", + " ax1.errorbar(x[xdata_mask], _spline_diff, yerr[xdata_mask], fmt=\"o\", label=filt)\n", + "\n", + " phase = np.concatenate([phase, x])\n", + " spline_diff = np.concatenate([spline_diff, _spline_diff])\n", + "\n", + "ad_g = adler_cols.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\")\n", + "ad_r = adler_cols.get_phase_parameters_in_filter(\"r\", \"HG12_Pen16\")\n", + "ax1.axhline(ad_g.H.value - ad_r.H.value, label=\"H_g - H_r\")\n", + "\n", + "spline_diff_med = np.median(spline_diff)\n", + "spline_diff_std = np.std(spline_diff)\n", + "ax1.axhline(spline_diff_med, label=\"median spline diff\", c=\"k\")\n", + "for i in [3, 2, 1]:\n", + " ax1.axhspan(\n", + " spline_diff_med - (i * spline_diff_std),\n", + " spline_diff_med + (i * spline_diff_std),\n", + " zorder=0,\n", + " color=\"k\",\n", + " alpha=0.15,\n", + " )\n", + "\n", + "ax1.legend()\n", + "ax1.set_xlabel(\"phaseAngle\")\n", + "ax1.set_ylabel(\"colour\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21d8adc0", + "metadata": {}, + "outputs": [], + "source": [ + "# account for photometric uncertainty!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8089e839", + "metadata": {}, + "outputs": [], + "source": [ + "# should we fit the spline in phase angle or time space?\n", + "# if time, should we consider apparitions separately?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1322e453", + "metadata": {}, + "outputs": [], + "source": [ + "# find apparitions, use time/solar elongation?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc139a24", + "metadata": {}, + "outputs": [], + "source": [ + "df_obs = pd.DataFrame(obs.__dict__)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2a5d219", + "metadata": {}, + "outputs": [], + "source": [ + "Time(np.amin(df_obs[\"midPointMjdTai\"]), format=\"mjd\").iso" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "532e01f6", + "metadata": {}, + "outputs": [], + "source": [ + "Time(np.amax(df_obs[\"midPointMjdTai\"]), format=\"mjd\").iso" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ac0599d", + "metadata": {}, + "outputs": [], + "source": [ + "R = df_obs[\"heliocentricDist\"]\n", + "Delta = df_obs[\"topocentricDist\"]\n", + "alpha = np.radians(df_obs[\"phaseAngle\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d88606c8", + "metadata": {}, + "outputs": [], + "source": [ + "# eps1 = (R * np.tan(alpha))\n", + "# eps2 = (Delta/np.cos(alpha)) - R\n", + "# df_obs[\"elong\"] = np.degrees(np.arctan(eps1/eps2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1521b6de", + "metadata": {}, + "outputs": [], + "source": [ + "# R_E = np.sqrt((R*R) + (Delta*Delta) - (2.0*R*Delta*np.cos(alpha)))\n", + "# df_obs[\"elong\"] = np.degrees(np.arcsin((R/R_E) * np.sin(alpha)))\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, + "id": "9c8f8630", + "metadata": {}, + "outputs": [], + "source": [ + "# _R = np.sqrt((R_E*R_E) + (Delta*Delta) - (2.0*R_E*Delta*np.cos(np.radians(df_obs[\"elong\"]))))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e73ba4bf", + "metadata": {}, + "outputs": [], + "source": [ + "df_obs = df_obs.sort_values(\"midPointMjdTai\")\n", + "df_obs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81314d44", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.hist(np.diff(df_obs[\"midPointMjdTai\"]), bins=\"auto\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "182dfdda", + "metadata": {}, + "outputs": [], + "source": [ + "_R = np.sqrt(\n", + " (df_obs[\"heliocentricX\"] ** 2.0) + (df_obs[\"heliocentricY\"] ** 2.0) + (df_obs[\"heliocentricZ\"] ** 2.0)\n", + ")\n", + "_Delta = np.sqrt(\n", + " (df_obs[\"topocentricX\"] ** 2.0) + (df_obs[\"topocentricY\"] ** 2.0) + (df_obs[\"topocentricZ\"] ** 2.0)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9505d8e5", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_obs[\"midPointMjdTai\"], df_obs[\"phaseAngle\"])\n", + "\n", + "ax1.set_xlabel(\"mjd\")\n", + "ax1.set_ylabel(\"phaseAngle\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ed95c44", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_obs[\"midPointMjdTai\"], df_obs[\"heliocentricDist\"])\n", + "ax1.axhline(q, label=\"perihelion\", ls=\":\")\n", + "ax1.axhline(Q, label=\"aphelion\", ls=\"--\")\n", + "\n", + "ax1.legend()\n", + "ax1.set_xlabel(\"mjd\")\n", + "ax1.set_ylabel(\"AU\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f63f1ef", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_obs[\"midPointMjdTai\"], df_obs[\"elong\"])\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba6258f6", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_obs[\"phaseAngle\"], df_obs[\"elong\"], c=df_obs[\"midPointMjdTai\"])\n", + "ax1.axhline(90)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2776289d", + "metadata": {}, + "outputs": [], + "source": [ + "np.amin(df_obs[\"elong\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15fce9da", + "metadata": {}, + "outputs": [], + "source": [ + "# astroquery the object\n", + "# target = \"2007 YP19\"\n", + "# target = \"2014 QL433\"\n", + "target = planetoid.MPCORB.fullDesignation\n", + "\n", + "site = \"X05\"\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, + "id": "f722a436", + "metadata": {}, + "outputs": [], + "source": [ + "target" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb1a5c67", + "metadata": {}, + "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, + "id": "44761d8e", + "metadata": {}, + "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, + "id": "3fecb8f2", + "metadata": {}, + "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, + "id": "3d0c1ab2", + "metadata": {}, + "outputs": [], + "source": [ + "df_eph[\"datetime_mjd\"] = Time(df_eph[\"datetime_jd\"], format=\"jd\").mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b9f978a", + "metadata": {}, + "outputs": [], + "source": [ + "df_vec[\"datetime_mjd\"] = Time(df_vec[\"datetime_jd\"], format=\"jd\").mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee69e367", + "metadata": {}, + "outputs": [], + "source": [ + "eph[\"datetime_jd\"].unit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1be1d08b", + "metadata": {}, + "outputs": [], + "source": [ + "df_eph.columns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "801ad4fd", + "metadata": {}, + "outputs": [], + "source": [ + "# plot only night time data points, replace day time with nans (will leave blank in plot)\n", + "# night_mask = (~np.isin(df_eph[\"solar_presence\"],[\"*\",\"N\",\"C\"])) & (df_eph[\"EL\"]>0)\n", + "night_mask = (~np.isin(df_eph[\"solar_presence\"], [\"*\", \"N\", \"C\"])) & (df_eph[\"EL\"] > 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43228dfa", + "metadata": {}, + "outputs": [], + "source": [ + "df_eph[[\"solar_presence\", \"EL\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e7ec662", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"],\n", + " df_obs[\"heliocentricDist\"],\n", + " edgecolor=\"C0\",\n", + " facecolor=\"none\",\n", + " label=\"heliocentricDist\",\n", + ")\n", + "# ax1.scatter(df_obs[\"midPointMjdTai\"],_R, marker = \"x\", c = \"C0\")\n", + "\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"],\n", + " df_obs[\"topocentricDist\"],\n", + " edgecolor=\"C1\",\n", + " facecolor=\"none\",\n", + " label=\"topocentricDist\",\n", + ")\n", + "# ax1.scatter(df_obs[\"midPointMjdTai\"],_Delta, marker = \"x\", c = \"C1\")\n", + "\n", + "ax1.scatter(df_obs[\"midPointMjdTai\"], R_E, edgecolor=\"C2\", facecolor=\"none\", label=\"heliocentricDist_Earth\")\n", + "# ax1.scatter(df_obs[\"midPointMjdTai\"],_R_E, marker = \"x\", c = \"C2\")\n", + "\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"r\"], c=\"C0\")\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"delta\"], c=\"C1\")\n", + "\n", + "ax1.legend()\n", + "ax1.set_xlabel(\"mjd\")\n", + "ax1.set_ylabel(\"AU\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51b6cf34", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"], df_obs[\"heliocentricX\"], edgecolor=\"C0\", facecolor=\"none\", label=\"heliocentricX\"\n", + ")\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"], df_obs[\"heliocentricY\"], edgecolor=\"C1\", facecolor=\"none\", label=\"heliocentricY\"\n", + ")\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"], df_obs[\"heliocentricZ\"], edgecolor=\"C2\", facecolor=\"none\", label=\"heliocentricZ\"\n", + ")\n", + "\n", + "ax1.plot(df_vec[\"datetime_mjd\"], df_vec[\"heliocentricX\"], c=\"C0\")\n", + "ax1.plot(df_vec[\"datetime_mjd\"], df_vec[\"heliocentricY\"], c=\"C1\")\n", + "ax1.plot(df_vec[\"datetime_mjd\"], df_vec[\"heliocentricZ\"], c=\"C2\")\n", + "\n", + "ax1.legend()\n", + "ax1.set_xlabel(\"mjd\")\n", + "ax1.set_ylabel(\"AU\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8208fe5f", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"], df_obs[\"topocentricX\"], edgecolor=\"C0\", facecolor=\"none\", label=\"topocentricX\"\n", + ")\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"], df_obs[\"topocentricY\"], edgecolor=\"C1\", facecolor=\"none\", label=\"topocentricY\"\n", + ")\n", + "ax1.scatter(\n", + " df_obs[\"midPointMjdTai\"], df_obs[\"topocentricZ\"], edgecolor=\"C2\", facecolor=\"none\", label=\"topocentricZ\"\n", + ")\n", + "\n", + "ax1.plot(df_vec[\"datetime_mjd\"], df_vec[\"topocentricX\"], c=\"C0\")\n", + "ax1.plot(df_vec[\"datetime_mjd\"], df_vec[\"topocentricY\"], c=\"C1\")\n", + "ax1.plot(df_vec[\"datetime_mjd\"], df_vec[\"topocentricZ\"], c=\"C2\")\n", + "\n", + "ax1.legend()\n", + "ax1.set_xlabel(\"mjd\")\n", + "ax1.set_ylabel(\"AU\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff0a02cb", + "metadata": {}, + "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, + "id": "fa4b8ab2", + "metadata": {}, + "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\", zorder=3)\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"alpha\"], c=\"k\", label=\"JPL\", zorder=0)\n", + "ax1.scatter(df_eph[night_mask][\"datetime_mjd\"], df_eph[night_mask][\"alpha\"], s=10)\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad8585b4", + "metadata": {}, + "outputs": [], + "source": [ + "df_eph[df_eph[\"alpha\"] > 120][[\"datetime_str\", \"alpha\", \"EL\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2eb7812", + "metadata": {}, + "outputs": [], + "source": [ + "Time(np.array(df_obs[df_obs[\"phaseAngle\"] > 120][[\"midPointMjdTai\"]]), format=\"mjd\").iso" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15c02874", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"alpha\"\n", + "y_plot = \"elong\"\n", + "c_plot = \"datetime_jd\"\n", + "df_plot = df_eph[night_mask]\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], c=df_plot[c_plot])\n", + "ax1.axhline(90)\n", + "\n", + "ax1.scatter(df_obs[\"phaseAngle\"], df_obs[\"elong\"], c=\"r\", marker=\"x\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33c4770e", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"r\"\n", + "y_plot = \"delta\"\n", + "c_plot = \"datetime_jd\"\n", + "df_plot = df_eph[night_mask]\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], c=df_plot[c_plot])\n", + "\n", + "ax1.scatter(df_obs[\"heliocentricDist\"], df_obs[\"topocentricDist\"], c=\"r\", marker=\"x\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a2a6a4d", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"r\"\n", + "y_plot = \"alpha\"\n", + "c_plot = \"datetime_jd\"\n", + "df_plot = df_eph[night_mask]\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], c=df_plot[c_plot])\n", + "\n", + "ax1.scatter(df_obs[\"heliocentricDist\"], df_obs[\"phaseAngle\"], c=\"r\", marker=\"x\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d6a9d9f", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"delta\"\n", + "y_plot = \"alpha\"\n", + "c_plot = \"datetime_jd\"\n", + "df_plot = df_eph[night_mask]\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], c=df_plot[c_plot])\n", + "\n", + "ax1.scatter(df_obs[\"topocentricDist\"], df_obs[\"phaseAngle\"], c=\"r\", marker=\"x\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c82e624", + "metadata": {}, + "outputs": [], + "source": [ + "# look up rsp notebooks to calculate x, y, z positions and to check orbits and angles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6004b00e", + "metadata": {}, + "outputs": [], + "source": [ + "df_obs_all = pd.DataFrame()\n", + "for filt in [\"r\", \"g\"]:\n", + " obs = planetoid.observations_in_filter(filt)\n", + " _df_obs = pd.DataFrame(obs.__dict__)\n", + " df_obs_all = pd.concat([df_obs_all, _df_obs])\n", + "df_obs_all = df_obs_all.sort_values(\"midPointMjdTai\")\n", + "\n", + "t_app = apparition_gap_finder(np.array(df_obs_all[\"midPointMjdTai\"]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8edcebd7", + "metadata": {}, + "outputs": [], + "source": [ + "df_obs_all" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9bcc201", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"heliocentricX\"\n", + "y_plot = \"heliocentricY\"\n", + "z_plot = \"heliocentricZ\"\n", + "df_plot = df_obs_all\n", + "\n", + "fig = plt.figure()\n", + "ax1 = fig.add_subplot(projection=\"3d\")\n", + "\n", + "ax1.scatter(df_plot[x_plot], df_plot[y_plot], df_plot[z_plot])\n", + "\n", + "ax1.set_aspect(\"equal\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97431295", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"topocentricX\"\n", + "y_plot = \"topocentricY\"\n", + "z_plot = \"topocentricZ\"\n", + "df_plot = df_obs_all\n", + "\n", + "fig = plt.figure()\n", + "ax1 = fig.add_subplot(projection=\"3d\")\n", + "\n", + "ax1.scatter(df_plot[x_plot], df_plot[y_plot], df_plot[z_plot])\n", + "\n", + "ax1.set_aspect(\"equal\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74f1c351", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"heliocentricX\"\n", + "y_plot = \"heliocentricY\"\n", + "c_plot = \"midPointMjdTai\"\n", + "df_plot = df_obs_all\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, + "id": "523e777a", + "metadata": {}, + "outputs": [], + "source": [ + "# use simple time diff for now" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66e70f52", + "metadata": {}, + "outputs": [], + "source": [ + "x = np.array(df_obs_all[\"midPointMjdTai\"])\n", + "x[0], x[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd09058a", + "metadata": {}, + "outputs": [], + "source": [ + "t_app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa2675d6", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"reduced_mag\"\n", + "df_plot = df_obs_all\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# for filt in [\"r\",\"g\"]:\n", + "for filt in [\"g\", \"r\"]:\n", + " _df_plot = df_plot[df_plot[\"filter_name\"] == filt]\n", + " ax1.scatter(_df_plot[x_plot], _df_plot[y_plot], label=filt)\n", + "\n", + "for x in t_app:\n", + " ax1.axvline(x, c=\"k\")\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6042ab2f", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "for i in range(len(t_app) - 1):\n", + " if i != 3:\n", + " continue\n", + "\n", + " time_min = t_app[i]\n", + " time_max = t_app[i + 1]\n", + " print(time_min, time_max)\n", + "\n", + " _df_obs_all = df_obs_all[\n", + " (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\n", + " print(i, len(_df_obs_all))\n", + " _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n", + "\n", + " splines = []\n", + " xdata = np.array([])\n", + " for filt in [\"g\", \"r\"]:\n", + " # get observations in filter\n", + " obs = planetoid.observations_in_filter(filt)\n", + " # get the stored AdlerData parameters for this filter\n", + " ad = adler_cols.get_phase_parameters_in_filter(filt, \"HG12_Pen16\")\n", + " # make the PhaseCurve object from AdlerData\n", + " pc = PhaseCurve().InitModelDict(ad.__dict__)\n", + " # print(pc.__dict__)\n", + "\n", + " # calculate the model absolute magnitude\n", + " abs_mag = pc.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value\n", + "\n", + " x = getattr(obs, x_plot)\n", + " xdata = np.concatenate([xdata, x])\n", + "\n", + " y = abs_mag\n", + " sort_mask = np.argsort(x)\n", + " x = x[sort_mask]\n", + " y = y[sort_mask]\n", + " yerr = obs.magErr[sort_mask]\n", + " # x_mask = (x>=time_min) & (x= time_min) & (x <= _time_max)\n", + " x = x[x_mask]\n", + " y = y[x_mask]\n", + " yerr = yerr[x_mask]\n", + " # print(len(x))\n", + "\n", + " ax1.errorbar(x, y, yerr, fmt=\"o\", label=filt)\n", + "\n", + " # ax1.axhline(0, c=\"k\")\n", + " ax1.axhline(pc.H.value, c=\"k\", label=\"H_{} = {:.2f} mag\".format(filt, pc.H.value))\n", + "\n", + " if len(x) < 5:\n", + " continue\n", + "\n", + " tck = splrep(x, y, w=1.0 / yerr)\n", + " splines.append(tck)\n", + "\n", + " # xnew = np.linspace(time_min,time_max,1000)\n", + " xnew = np.linspace(time_min, _time_max, 1000)\n", + " ax1.plot(xnew, BSpline(*tck)(xnew), \"-\", label=\"weighted smooth\", c=\"k\", zorder=3)\n", + "\n", + "# ax1.legend()\n", + "ax1.set_ylim(pc.H.value - 1, pc.H.value + 1)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d23ecf23", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "for i in range(len(t_app) - 1):\n", + " if i != 3:\n", + " continue\n", + "\n", + " time_min = t_app[i]\n", + " time_max = t_app[i + 1]\n", + " print(time_min, time_max)\n", + "\n", + " _df_obs_all = df_obs_all[\n", + " (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\n", + " print(i, len(_df_obs_all))\n", + " _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n", + "\n", + " splines = []\n", + " xdata = np.array([])\n", + " for filt in [\"g\", \"r\"]:\n", + " # get observations in filter\n", + " obs = planetoid.observations_in_filter(filt)\n", + " # get the stored AdlerData parameters for this filter\n", + " ad = adler_cols.get_phase_parameters_in_filter(filt, \"HG12_Pen16\")\n", + " # make the PhaseCurve object from AdlerData\n", + " pc = PhaseCurve().InitModelDict(ad.__dict__)\n", + " # print(pc.__dict__)\n", + "\n", + " # calculate the model absolute magnitude\n", + " abs_mag = pc.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value\n", + "\n", + " x = getattr(obs, x_plot)\n", + " xdata = np.concatenate([xdata, x])\n", + "\n", + " y = abs_mag\n", + " sort_mask = np.argsort(x)\n", + " x = x[sort_mask]\n", + " y = y[sort_mask]\n", + " yerr = obs.magErr[sort_mask]\n", + " # x_mask = (x>=time_min) & (x= time_min) & (x <= _time_max)\n", + " x = x[x_mask]\n", + " y = y[x_mask]\n", + " yerr = yerr[x_mask]\n", + " # print(len(x))\n", + "\n", + " if len(x) < 5:\n", + " continue\n", + "\n", + " tck = splrep(x, y, w=1.0 / yerr)\n", + " splines.append(tck)\n", + "\n", + " if len(splines) < 2:\n", + " continue\n", + "\n", + " # xnew = np.linspace(time_min,time_max,1000)\n", + " xnew = np.linspace(time_min, _time_max, 1000)\n", + " spline_diff = BSpline(*splines[0])(xnew) - BSpline(*splines[1])(xnew)\n", + " ax1.plot(xnew, spline_diff, \"-\", label=\"weighted smooth\", c=\"k\", zorder=1)\n", + "\n", + " phase = np.array([])\n", + " spline_diff = np.array([])\n", + " for i, filt in enumerate([\"g\", \"r\"]):\n", + " obs = planetoid.observations_in_filter(filt)\n", + " x = getattr(obs, \"midPointMjdTai\")\n", + " yerr = obs.magErr\n", + " # xdata_mask = (x>=time_min) & (x= time_min) & (x <= _time_max)\n", + " _spline_diff = BSpline(*splines[0])(x[xdata_mask]) - BSpline(*splines[1])(x[xdata_mask])\n", + " ax1.errorbar(x[xdata_mask], _spline_diff, yerr[xdata_mask], fmt=\"o\", label=filt)\n", + "\n", + " phase = np.concatenate([phase, x])\n", + " spline_diff = np.concatenate([spline_diff, _spline_diff])\n", + "\n", + " ad_g = adler_cols.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\")\n", + " ad_r = adler_cols.get_phase_parameters_in_filter(\"r\", \"HG12_Pen16\")\n", + " ax1.axhline(ad_g.H.value - ad_r.H.value, label=\"H_g - H_r\")\n", + "\n", + "# spline_diff_med = np.median(spline_diff)\n", + "# spline_diff_std = np.std(spline_diff)\n", + "# ax1.axhline(spline_diff_med,label = \"median spline diff\", c = \"k\")\n", + "# for i in [3,2,1]:\n", + "# ax1.axhspan(spline_diff_med - (i*spline_diff_std), spline_diff_med + (i*spline_diff_std),\n", + "# zorder = 0, color = \"k\", alpha = 0.15)\n", + "\n", + "# ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb06fef2", + "metadata": {}, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "for i in range(len(t_app) - 1):\n", + " if i != 3:\n", + " continue\n", + "\n", + " time_min = t_app[i]\n", + " time_max = t_app[i + 1]\n", + " print(time_min, time_max)\n", + "\n", + " _df_obs_all = df_obs_all[\n", + " (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\n", + " print(i, len(_df_obs_all))\n", + " _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n", + "\n", + " for filt in [\"g\", \"r\"]:\n", + " # get observations in filter\n", + " obs = planetoid.observations_in_filter(filt)\n", + "\n", + " y = getattr(obs, y_plot)\n", + " yerr = obs.magErr\n", + " x = getattr(obs, x_plot)\n", + "\n", + " # y = abs_mag\n", + " sort_mask = np.argsort(x)\n", + " x = x[sort_mask]\n", + " y = y[sort_mask]\n", + " yerr = obs.magErr[sort_mask]\n", + " # # x_mask = (x>=time_min) & (x= time_min) & (x <= _time_max)\n", + " x = x[x_mask]\n", + " y = y[x_mask]\n", + " yerr = yerr[x_mask]\n", + " # # print(len(x))\n", + "\n", + " ax1.errorbar(x, y, yerr, fmt=\"o\", label=filt)\n", + "\n", + " tck = splrep(x, y, w=1.0 / yerr)\n", + " splines.append(tck)\n", + "\n", + " xnew = np.linspace(time_min, _time_max, 1000)\n", + " ax1.plot(xnew, BSpline(*tck)(xnew), \"-\", label=\"weighted smooth\", c=\"k\", zorder=3)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "be147dff", + "metadata": {}, + "source": [ + "# test colour finding functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e16f9604", + "metadata": {}, + "outputs": [], + "source": [ + "# # set number of reference observations to use for colour estimate\n", + "# # N_mean = 5\n", + "# N_mean = 3\n", + "# # N_mean = 1\n", + "# # N_mean = None\n", + "\n", + "# # define observation column names to be analysed for colour\n", + "# x_plot = \"midPointMjdTai\"\n", + "# y_plot = \"AbsMag\"\n", + "# # y_plot = \"reduced_mag\"\n", + "# yerr_plot = \"magErr\"\n", + "# filt_obs = \"g\"\n", + "# filt_ref = \"r\"\n", + "# colour = \"{}-{}\".format(filt_obs,filt_ref)\n", + "# colErr = \"{}-{}Err\".format(filt_obs,filt_ref)\n", + "# delta_t_col = \"delta_t_{}\".format(colour)\n", + "\n", + "# # define new columns to be added to the observations dataframe\n", + "# new_obs_cols = [colour,colErr,delta_t_col]\n", + "\n", + "# fig = plt.figure()\n", + "# gs = gridspec.GridSpec(1, 1)\n", + "# ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# for i in range(len(t_app)-1):\n", + "\n", + "# if i!=3:\n", + "# continue\n", + "\n", + "# time_min = t_app[i]\n", + "# time_max = t_app[i+1]\n", + "# print(time_min,time_max)\n", + "\n", + "# _df_obs_all = df_obs_all[(df_obs_all[\"midPointMjdTai\"]>=time_min) & (df_obs_all[\"midPointMjdTai\"]= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\n", + " print(app_i, len(_df_obs_all))\n", + " _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n", + " print(time_min, _time_max)\n", + "\n", + " # get the phase curve model and observations for each filter\n", + " # get the stored AdlerData parameters for this filter\n", + " ad_g = adler_cols.get_phase_parameters_in_filter(filt_obs, \"HG12_Pen16\")\n", + " # make the PhaseCurve object from AdlerData\n", + " pc_g = PhaseCurve().InitModelDict(ad_g.__dict__)\n", + " ad_r = adler_cols.get_phase_parameters_in_filter(filt_ref, \"HG12_Pen16\")\n", + " pc_r = PhaseCurve().InitModelDict(ad_r.__dict__)\n", + " df_obs = get_df_obs_filt(\n", + " planetoid,\n", + " filt_obs,\n", + " x1=time_min,\n", + " x2=_time_max,\n", + " col_list=[obsId_col, y_plot, yerr_plot],\n", + " pc_model=pc_g,\n", + " )\n", + " df_obs_ref = get_df_obs_filt(\n", + " planetoid, filt_ref, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_r\n", + " )\n", + "\n", + " ax1.errorbar(df_obs[x_plot], df_obs[y_plot], df_obs[yerr_plot], fmt=\"o\", label=filt_obs)\n", + " ax1.errorbar(df_obs_ref[x_plot], df_obs_ref[y_plot], df_obs_ref[yerr_plot], fmt=\"o\", label=filt_ref)\n", + "\n", + " # simulate stepping through each filt_obs observation\n", + " x1 = time_min\n", + " for xi in range(len(df_obs)):\n", + " x2 = df_obs.iloc[xi][x_plot]\n", + " print(xi, x1, x2)\n", + "\n", + " # do the colour finding function here\n", + " col_dict = col_obs_ref(\n", + " planetoid,\n", + " adler_cols,\n", + " filt_obs=filt_obs,\n", + " filt_ref=filt_ref,\n", + " N_ref=N_ref,\n", + " x_col=x_plot,\n", + " y_col=y_plot,\n", + " yerr_col=yerr_plot,\n", + " x1=x1,\n", + " x2=x2,\n", + " )\n", + " col_dict_list.append(col_dict)\n", + "\n", + " # plot some lines to show the colour and mean reference\n", + " ax1.vlines(\n", + " df_obs.iloc[xi][x_plot], df_obs.iloc[xi][y_plot], col_dict[y_ref_col], color=\"k\", ls=\":\"\n", + " )\n", + " ax1.hlines(col_dict[y_ref_col], col_dict[x1_ref_col], col_dict[x2_ref_col], color=\"k\", ls=\"--\")\n", + "\n", + " df_col = pd.DataFrame(col_dict_list)\n", + " print(list(df_col))\n", + " print(list(df_obs))\n", + " # df_col = pd.concat([df_obs,df_col], axis =1)\n", + " df_col = df_col.merge(df_obs[[x for x in list(df_obs) if x != x_plot]], on=obsId_col)\n", + " df_col[\"diaSourceId\"] = df_col[\"diaSourceId\"].astype(int)\n", + "\n", + " ax1.set_xlabel(x_plot)\n", + " ax1.set_ylabel(y_plot)\n", + " ax1.legend()\n", + " ax1.invert_yaxis()\n", + "\n", + " plt.show()\n", + "\n", + " df_col_fname = \"df_{}_{}_{}_app_{}_N_ref_{}.csv\".format(ssoid, filt_obs, filt_ref, test_app_i, N_ref)\n", + " df_col.to_csv(df_col_fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca5b9204", + "metadata": {}, + "outputs": [], + "source": [ + "df_obs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29b9b35f", + "metadata": {}, + "outputs": [], + "source": [ + "col_dict" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a4d24f2", + "metadata": {}, + "outputs": [], + "source": [ + "df_col" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bacf848e", + "metadata": {}, + "outputs": [], + "source": [ + "# df_col_fname = \"df_{}_{}_{}_app_{}_N_ref_{}.csv\".format(ssoid,filt_obs,filt_ref,test_app_i,N_ref)\n", + "# df_col_fname" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "299bba1a", + "metadata": {}, + "outputs": [], + "source": [ + "# df_col.to_csv(df_col_fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad570220", + "metadata": {}, + "outputs": [], + "source": [ + "# find filt_obs - filt_ref of newest filt_obs observation to the mean of the previous N_ref filt_ref observations\n", + "# colour code by time diff between obs and most recent obs_ref\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = colour\n", + "y_plot_err = colErr\n", + "c_plot = delta_t_col\n", + "df_plot = df_col\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "\n", + "s1 = ax1.scatter(df_plot[x_plot], df_plot[y_plot], c=df_plot[c_plot], zorder=3)\n", + "cbar1 = plt.colorbar(s1)\n", + "ax1.errorbar(df_plot[x_plot], df_plot[y_plot], df_plot[yerr_plot], fmt=\".\", zorder=1)\n", + "\n", + "# TODO: leave the most recent obs out of this mean?\n", + "# loop through new obs and calculate current mean and std?\n", + "obs_ref_mean = np.mean(df_plot[y_plot])\n", + "obs_ref_median = np.median(df_plot[y_plot])\n", + "obs_ref_std = np.std(df_plot[y_plot])\n", + "print(obs_ref_median, obs_ref_std)\n", + "\n", + "ax1.axhline(obs_ref_mean, c=\"k\")\n", + "ax1.axhspan(obs_ref_mean - obs_ref_std, obs_ref_mean + obs_ref_std, zorder=0, color=\"k\", alpha=0.2)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "cbar1.set_label(c_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7f5253e", + "metadata": {}, + "outputs": [], + "source": [ + "# The colours can then be run through the previously written outlier detection functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6cd89e91", + "metadata": {}, + "outputs": [], + "source": [ + "planetoid.ssObjectId" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2246f74", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "adler-dev", + "language": "python", + "name": "adler-dev" + }, + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}