diff --git a/docs/notebooks.rst b/docs/notebooks.rst index 7f7e544..61e3308 100644 --- a/docs/notebooks.rst +++ b/docs/notebooks.rst @@ -4,3 +4,4 @@ Notebooks .. toctree:: Introducing Jupyter Notebooks + Adler phasecurve models diff --git a/docs/notebooks/adler_phasecurve_example.ipynb b/docs/notebooks/adler_phasecurve_example.ipynb new file mode 100644 index 0000000..6e64763 --- /dev/null +++ b/docs/notebooks/adler_phasecurve_example.ipynb @@ -0,0 +1,359 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8738000d", + "metadata": {}, + "source": [ + "# Adler phasecurve models\n", + "This notebook demonstrates how Adler implements phasecurve models. An example object with photometric observations is loaded. We can create a phasecurve model object from the SSObject parameters associated with this object. We can also fit a phasecurve model of our choice to the observations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d591f5d8-9148-46ff-a62b-0f2a29eb806c", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import astropy.units as u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "645efb98-567d-481e-a79c-b1cfdc828726", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"8268570668335894776\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10b36aab-b322-49b8-8ff3-49bef68d7416", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the object data via adler\n", + "\n", + "# # here we use an offline SQL database which contains the observations of the sso\n", + "fname = \"../../tests/data/testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)\n", + "\n", + "# alternatively we can retrieve the object data directly from the RSP\n", + "# planetoid = AdlerPlanetoid.construct_from_RSP(ssoid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9a0623d-0dc7-49c1-99dd-a76ef970a3ff", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the whole object\n", + "planetoid.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d360360-025b-4a77-acf5-325b2f2d1873", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect just the ssObject table\n", + "planetoid.SSObject.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be8f8d63", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve all observations in the r filter\n", + "obs_r = planetoid.observations_in_filter(\"r\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da1e483d", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the fields available in the observations table\n", + "obs_r.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d7dc125-06c1-49ad-8854-17d8c8b6954f", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations as a phasecurve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6efe3b5a-09dd-4d5e-9f41-20ea6e1b43df", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the phase curve model parameters provided in the ssObject table\n", + "\n", + "sso_r = planetoid.SSObject_in_filter(\"r\")\n", + "\n", + "r_H = sso_r.H\n", + "r_G12 = sso_r.G12\n", + "\n", + "pc = PhaseCurve(H=r_H * u.mag, phase_parameter_1=r_G12, model_name=\"HG12_Pen16\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80f552f1-8907-4cc9-b57c-2e667eab459c", + "metadata": {}, + "outputs": [], + "source": [ + "# what sbpy model is being used?\n", + "pc.model_function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24c1955e-95cd-4d77-ad05-aa5b8d18620a", + "metadata": {}, + "outputs": [], + "source": [ + "# set up an array of phase angles to plot the model\n", + "alpha = np.linspace(0, np.amax(obs_r.phaseAngle)) * u.deg\n", + "alpha" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3f30fe0-0d89-4ffa-8237-9c71181d44ee", + "metadata": {}, + "outputs": [], + "source": [ + "# calculate the model reduced magnitude over these phase angles\n", + "red_mag = pc.ReducedMag(alpha)\n", + "red_mag" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04be98a1-e4dc-4216-bcd9-ef777f6053fb", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations with the model phase curve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.plot(alpha.value, red_mag.value)\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9815543d-6140-4bdb-8bad-8296994723f4", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations as a lightcurve\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de462b92-3914-4091-b0af-bddd9e9c1ef1", + "metadata": {}, + "outputs": [], + "source": [ + "# do a different phase curve fit to the data\n", + "# adler is able to fit different models, and perform more sophisticated fits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f92891c9-6ccf-4dac-8887-9545f633ba90", + "metadata": {}, + "outputs": [], + "source": [ + "# create a new PhaseCurve object with a different sbpy model\n", + "pc_fit = PhaseCurve(H=pc.H, model_name=\"HG\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db24432b-6d05-4ff2-9d98-e52d8c2e4342", + "metadata": {}, + "outputs": [], + "source": [ + "pc_fit.model_function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9039e2e2-27d9-4d21-b2f6-9504a5b85ce4", + "metadata": {}, + "outputs": [], + "source": [ + "# use adler to fit this new phase curve model to the data\n", + "pc_fit.FitModel(\n", + " phase_angle=obs_r.phaseAngle * u.deg,\n", + " reduced_mag=obs_r.reduced_mag * u.mag,\n", + " mag_err=obs_r.magErr * u.mag,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d4b7144-ee72-45e0-9606-c40f83c443c6", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations with both\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "xerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, xerr, fmt=\"o\")\n", + "\n", + "ax1.plot(alpha.value, pc.ReducedMag(alpha).value, label=pc.model_name)\n", + "ax1.plot(alpha.value, pc_fit.ReducedMag(alpha).value, label=pc_fit.model_name)\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7f39ed4-8334-4e10-a97c-a9471105225b", + "metadata": {}, + "outputs": [], + "source": [ + "# # now we would add our calculated values back into planetoid\n", + "# planetoid.AdlerSchema.r_H = pc_fit.abs_mag\n", + "# planetoid.AdlerSchema.r_G = pc_fit.phase_param" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ca4bbfd-1954-469f-8608-40c52838d300", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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 +} diff --git a/docs/requirements.txt b/docs/requirements.txt index 3979f83..a02ca08 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -8,3 +8,6 @@ jupyter astropy sbpy matplotlib +numpy +ipykernel +scipy \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 7482ed1..7874b80 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,6 +21,7 @@ dependencies = [ "sbpy", "matplotlib", # for plotting "pandas", + "scipy", ] [project.urls]