diff --git a/notebooks/adler_examples.ipynb b/notebooks/adler_examples.ipynb new file mode 100644 index 0000000..e594f29 --- /dev/null +++ b/notebooks/adler_examples.ipynb @@ -0,0 +1,489 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "d591f5d8-9148-46ff-a62b-0f2a29eb806c", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:11.533501Z", + "iopub.status.busy": "2024-02-29T16:09:11.533247Z", + "iopub.status.idle": "2024-02-29T16:09:13.062886Z", + "shell.execute_reply": "2024-02-29T16:09:13.061510Z", + "shell.execute_reply.started": "2024-02-29T16:09:11.533477Z" + } + }, + "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": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:13.065333Z", + "iopub.status.busy": "2024-02-29T16:09:13.064482Z", + "iopub.status.idle": "2024-02-29T16:09:13.070452Z", + "shell.execute_reply": "2024-02-29T16:09:13.069235Z", + "shell.execute_reply.started": "2024-02-29T16:09:13.065297Z" + } + }, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"8268570668335894776\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10b36aab-b322-49b8-8ff3-49bef68d7416", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:13.072314Z", + "iopub.status.busy": "2024-02-29T16:09:13.071967Z", + "iopub.status.idle": "2024-02-29T16:09:14.279820Z", + "shell.execute_reply": "2024-02-29T16:09:14.278772Z", + "shell.execute_reply.started": "2024-02-29T16:09:13.072275Z" + } + }, + "outputs": [], + "source": [ + "# retrieve the object data via adler\n", + "planetoid = AdlerPlanetoid(ssoid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9a0623d-0dc7-49c1-99dd-a76ef970a3ff", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.281553Z", + "iopub.status.busy": "2024-02-29T16:09:14.281143Z", + "iopub.status.idle": "2024-02-29T16:09:14.287399Z", + "shell.execute_reply": "2024-02-29T16:09:14.286612Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.281511Z" + } + }, + "outputs": [], + "source": [ + "# inspect the object\n", + "planetoid.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d360360-025b-4a77-acf5-325b2f2d1873", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.290960Z", + "iopub.status.busy": "2024-02-29T16:09:14.290631Z", + "iopub.status.idle": "2024-02-29T16:09:14.307210Z", + "shell.execute_reply": "2024-02-29T16:09:14.306215Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.290937Z" + } + }, + "outputs": [], + "source": [ + "planetoid.SSObject.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0bcad0b-32ea-4c0f-b489-00d7491ff3b6", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.308956Z", + "iopub.status.busy": "2024-02-29T16:09:14.308489Z", + "iopub.status.idle": "2024-02-29T16:09:14.333001Z", + "shell.execute_reply": "2024-02-29T16:09:14.332167Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.308917Z" + } + }, + "outputs": [], + "source": [ + "# dir(planetoid.Observations)\n", + "planetoid.Observations.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b226fefa-a252-40d2-925e-549eee16858e", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.334599Z", + "iopub.status.busy": "2024-02-29T16:09:14.334226Z", + "iopub.status.idle": "2024-02-29T16:09:14.351732Z", + "shell.execute_reply": "2024-02-29T16:09:14.350910Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.334562Z" + } + }, + "outputs": [], + "source": [ + "getattr(planetoid.Observations, \"phaseAngle\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d7dc125-06c1-49ad-8854-17d8c8b6954f", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.353538Z", + "iopub.status.busy": "2024-02-29T16:09:14.352944Z", + "iopub.status.idle": "2024-02-29T16:09:14.544075Z", + "shell.execute_reply": "2024-02-29T16:09:14.543137Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.353500Z" + } + }, + "outputs": [], + "source": [ + "# plot the observations\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(planetoid.Observations, x_plot)\n", + "y = getattr(planetoid.Observations, y_plot)\n", + "xerr = planetoid.Observations.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": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.545627Z", + "iopub.status.busy": "2024-02-29T16:09:14.545342Z", + "iopub.status.idle": "2024-02-29T16:09:14.553323Z", + "shell.execute_reply": "2024-02-29T16:09:14.552540Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.545603Z" + } + }, + "outputs": [], + "source": [ + "# define the phase curve\n", + "pc = PhaseCurve(\n", + " abs_mag=planetoid.SSObject.r_H * u.mag, phase_param=planetoid.SSObject.r_G12, model_name=\"HG12_Pen16\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80f552f1-8907-4cc9-b57c-2e667eab459c", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.554644Z", + "iopub.status.busy": "2024-02-29T16:09:14.554319Z", + "iopub.status.idle": "2024-02-29T16:09:14.571588Z", + "shell.execute_reply": "2024-02-29T16:09:14.570921Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.554620Z" + } + }, + "outputs": [], + "source": [ + "pc.model_function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24c1955e-95cd-4d77-ad05-aa5b8d18620a", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.572987Z", + "iopub.status.busy": "2024-02-29T16:09:14.572491Z", + "iopub.status.idle": "2024-02-29T16:09:14.593610Z", + "shell.execute_reply": "2024-02-29T16:09:14.592732Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.572960Z" + } + }, + "outputs": [], + "source": [ + "alpha = np.linspace(0, np.amax(planetoid.Observations.phaseAngle)) * u.deg\n", + "alpha" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3f30fe0-0d89-4ffa-8237-9c71181d44ee", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.594977Z", + "iopub.status.busy": "2024-02-29T16:09:14.594621Z", + "iopub.status.idle": "2024-02-29T16:09:14.614140Z", + "shell.execute_reply": "2024-02-29T16:09:14.613401Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.594943Z" + } + }, + "outputs": [], + "source": [ + "red_mag = pc.ReducedMag(alpha)\n", + "red_mag" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04be98a1-e4dc-4216-bcd9-ef777f6053fb", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.615656Z", + "iopub.status.busy": "2024-02-29T16:09:14.615050Z", + "iopub.status.idle": "2024-02-29T16:09:14.783515Z", + "shell.execute_reply": "2024-02-29T16:09:14.782664Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.615623Z" + } + }, + "outputs": [], + "source": [ + "# %matplotlib widget\n", + "\n", + "# plot the observations with the LSST phase curve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(planetoid.Observations, x_plot)\n", + "y = getattr(planetoid.Observations, y_plot)\n", + "xerr = planetoid.Observations.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": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.784754Z", + "iopub.status.busy": "2024-02-29T16:09:14.784498Z", + "iopub.status.idle": "2024-02-29T16:09:14.953022Z", + "shell.execute_reply": "2024-02-29T16:09:14.952237Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.784731Z" + } + }, + "outputs": [], + "source": [ + "# plot the observations\n", + "x_plot = \"mjd\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(planetoid.Observations, x_plot)\n", + "y = getattr(planetoid.Observations, y_plot)\n", + "xerr = planetoid.Observations.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": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.954181Z", + "iopub.status.busy": "2024-02-29T16:09:14.953935Z", + "iopub.status.idle": "2024-02-29T16:09:14.957556Z", + "shell.execute_reply": "2024-02-29T16:09:14.956913Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.954159Z" + } + }, + "outputs": [], + "source": [ + "# do a different phase curve fit to the data\n", + "# adler should be able to fit different models, and perform more sophisticated fits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f92891c9-6ccf-4dac-8887-9545f633ba90", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.963464Z", + "iopub.status.busy": "2024-02-29T16:09:14.958584Z", + "iopub.status.idle": "2024-02-29T16:09:14.975015Z", + "shell.execute_reply": "2024-02-29T16:09:14.974187Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.963434Z" + } + }, + "outputs": [], + "source": [ + "pc_fit = PhaseCurve(abs_mag=pc.abs_mag, model_name=\"HG\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db24432b-6d05-4ff2-9d98-e52d8c2e4342", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.976283Z", + "iopub.status.busy": "2024-02-29T16:09:14.975943Z", + "iopub.status.idle": "2024-02-29T16:09:14.992854Z", + "shell.execute_reply": "2024-02-29T16:09:14.992071Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.976260Z" + } + }, + "outputs": [], + "source": [ + "pc_fit.model_function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9039e2e2-27d9-4d21-b2f6-9504a5b85ce4", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:14.994511Z", + "iopub.status.busy": "2024-02-29T16:09:14.993863Z", + "iopub.status.idle": "2024-02-29T16:09:15.023285Z", + "shell.execute_reply": "2024-02-29T16:09:15.022583Z", + "shell.execute_reply.started": "2024-02-29T16:09:14.994486Z" + } + }, + "outputs": [], + "source": [ + "pc_fit.FitModel(\n", + " phase_angle=planetoid.Observations.phaseAngle * u.deg,\n", + " reduced_mag=planetoid.Observations.reduced_mag * u.mag,\n", + " mag_err=planetoid.Observations.magErr * u.mag,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d4b7144-ee72-45e0-9606-c40f83c443c6", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:15.024555Z", + "iopub.status.busy": "2024-02-29T16:09:15.024229Z", + "iopub.status.idle": "2024-02-29T16:09:15.217263Z", + "shell.execute_reply": "2024-02-29T16:09:15.216519Z", + "shell.execute_reply.started": "2024-02-29T16:09:15.024533Z" + } + }, + "outputs": [], + "source": [ + "# %matplotlib widget\n", + "\n", + "# plot the observations with the LSST phase curve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(planetoid.Observations, x_plot)\n", + "y = getattr(planetoid.Observations, y_plot)\n", + "xerr = planetoid.Observations.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": { + "execution": { + "iopub.execute_input": "2024-02-29T16:09:15.221317Z", + "iopub.status.busy": "2024-02-29T16:09:15.220742Z", + "iopub.status.idle": "2024-02-29T16:09:15.224414Z", + "shell.execute_reply": "2024-02-29T16:09:15.223697Z", + "shell.execute_reply.started": "2024-02-29T16:09:15.221290Z" + } + }, + "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" + ] + } + ], + "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.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/adler/adler.py b/src/adler/adler.py index e73ab8c..c327aa3 100644 --- a/src/adler/adler.py +++ b/src/adler/adler.py @@ -1,6 +1,8 @@ import argparse +import astropy.units as u from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid +from adler.science.PhaseCurve import PhaseCurve def runAdler(args): @@ -8,6 +10,27 @@ def runAdler(args): planetoid.do_pretend_science() + # now let's do some phase curves! + + # get the RSP r filter model + pc = PhaseCurve( + abs_mag=planetoid.SSObject.H[2] * u.mag, + phase_param=planetoid.SSObject.G12[2], + model_name="HG12_Pen16", + ) + print(pc) + print(pc.abs_mag, pc.phase_param) + + # get the r filter observations + obs_r = planetoid.observations_in_filter("r") + alpha = obs_r.phaseAngle * u.deg + red_mag = obs_r.reduced_mag * u.mag + mag_err = obs_r.magErr * u.mag + + # do a simple fit to all data + pc_fit = pc.FitModel(alpha, red_mag, mag_err) + print(pc_fit) + def main(): parser = argparse.ArgumentParser(description="Runs Adler for a select planetoid and given user input.") diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py new file mode 100644 index 0000000..cbc328a --- /dev/null +++ b/src/adler/science/PhaseCurve.py @@ -0,0 +1,184 @@ +from sbpy.photometry import HG, HG1G2, HG12, HG12_Pen16, LinearPhaseFunc +import astropy.units as u +from astropy.modeling.fitting import LevMarLSQFitter + + +class PhaseCurve: + """A class to define the phasecurve model and associated functions. + Units - by default no units are set but astropy units can be passed. + It is up to the user to ensure that units are correct for the relevant phasecurve model. + + Attibutes + ----------- + abs_mag : float + Absolute magnitude, H, of the phasecurve model (often units of mag). + + phase_param: float + The first phase parameter of the phasecurve model, e.g. G from HG + (often dimensionless units unless S from LinearPhaseFunc, which has units mag/deg or mag/rad). + + phase_param2: float + The second phase parameter, only used for the 3 parameter HG1G2 phasecurve model. + + model_name: str + Label for the phasecurve model to be used. + Choice of: "HG", "HG1G2", "HG12", "HG12_Pen16", "LinearPhaseFunc" + + """ + + def __init__(self, abs_mag=18, phase_param=0.2, phase_param2=None, model_name="HG"): + self.abs_mag = abs_mag + self.phase_param = phase_param + self.phase_param2 = phase_param2 + self.model_name = model_name + + if model_name == "HG": + self.model_function = HG(H=abs_mag, G=self.phase_param) + elif model_name == "HG1G2": + self.model_function = HG1G2(H=abs_mag, G1=self.phase_param, G2=self.phase_param) + elif model_name == "HG12": + self.model_function = HG12(H=abs_mag, G12=self.phase_param) + elif model_name == "HG12_Pen16": + self.model_function = HG12_Pen16(H=abs_mag, G12=self.phase_param) + elif model_name == "LinearPhaseFunc": + self.model_function = LinearPhaseFunc(H=abs_mag, S=self.phase_param) + else: + print("no model selected") + + def ReturnModelDict(self): + """Return the values for the PhaseCurve class as a dict + + Returns + ---------- + + self.__dict__ : dict + The dict of PhaseCurve object parameters. + + """ + + return self.__dict__ + + def InitModelDict(self, model_dict): + """Set up a new PhaseCurve model object from a dictionary. + This could be written by the user or generated from another PhaseCurve object using ReturnModelDict + + Parameters + ----------- + model_dict : dict + Dictionary containing the PhaseCurve parameters you wish to set, e.g. abs_mag, phase_param + + Returns + ---------- + + model : object + The new PhaseCurve class object + + """ + + model = PhaseCurve() + for key, value in model_dict.items(): + setattr(model, key, value) + return model + + def InitModelSbpy(self, model_sbpy): + """Set up a new PhaseCurve model object from an existing sbpy model + ### or create dict from sbpy model and then use InitModelDict? + + Parameters + ----------- + model_sbpy : object + The sbpy model object, e.g. HG() + + Returns + ---------- + + model : object + The new PhaseCurve class object + + """ + # get model name from the sbpy model object + model_name = model_sbpy.__class__.name + + # get the sbpy model parameters + param_names = model_sbpy.param_names + parameters = [] + for p in param_names: + # try get the quantity (value with units) + x = getattr(model_sbpy, p).quantity + # if there are no units get just the value + if x is None: + x = getattr(model_sbpy, p).value + parameters.append(x) + print(param_names, parameters) + + # create a PhaseCurve object with the extracted parameters + model = PhaseCurve(*parameters, model_name=model_name) + + return model + + def ReducedMag(self, phase_angle): + """Return the reduced magnitude of the phasecurve model for a given phase angle(s) + + phase_angle - value or array, must have astropy units of degrees + + Parameters + ----------- + phase_angle : float or array + value or array of phase angles at which to evaluate the phasecurve model, must have astropy units of degrees. + + Returns + ---------- + + return_value : float or array + The phasecurve model reduced magnitude at the given phase angle(s) + + """ + + return self.model_function(phase_angle) + + def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): + """Fit the phasecurve model parameters to observations. + starts with a phase curve model as an initial guess for parameters. + fits model to phase angle and reduced magnitude. + + phase_angle - phase angle of each observations + reduced_mag - distance corrected reduced magnitudes + mag_err - photometric uncertainties to weight the measurements + fitter - can pass a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter + + Parameters + ----------- + phase_angle : float or array + The Sun-object-observer phase angles of the observations. + + reduced_mag : float or array + The observed reduced magnitudes at the corresponding phase angles. + + mag_err : float or array + Uncertainty on the reduced magnitude, used to weight the fit. + + fitter : object + Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter + + Returns + ---------- + + return_value : float or array + The phasecurve model reduced magnitude at the given phase angle(s) + + """ + + # use the LevMarLSQFitter by default + if fitter is None: + fitter = LevMarLSQFitter() + print(fitter) + + if mag_err is not None: # fit weighted by photometric uncertainty + model_fit = fitter(self.model_function, phase_angle, reduced_mag, weights=1.0 / mag_err) + else: # unweighted fit + model_fit = fitter(self.model_function, phase_angle, reduced_mag) + + ### if overwrite_model: # add an overwrite option? + # redo __init__ with the new fitted parameters + + return model_fit diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py new file mode 100644 index 0000000..f95b10d --- /dev/null +++ b/tests/adler/science/test_PhaseCurve.py @@ -0,0 +1,57 @@ +from numpy.testing import assert_array_equal +import pytest + + +def test_PhaseCurve_init(): + import numpy as np + import astropy.units as u + from adler.science.PhaseCurve import PhaseCurve + + # test creating a model + H = 18.9 + G = 0.12 + pc = PhaseCurve(abs_mag=H * u.mag, phase_param=G, model_name="HG") + + assert pc.abs_mag.value == 18.9 + assert pc.abs_mag.unit == u.mag + assert pc.phase_param == 0.12 + assert pc.model_name == "HG" + + +def test_PhaseCurve_ReducedMag(): + import numpy as np + import astropy.units as u + from adler.science.PhaseCurve import PhaseCurve + + # define the phase angles + alpha = np.array([0, 10]) * u.deg + + # linear phase curve model + pc_lin = PhaseCurve(model_name="LinearPhaseFunc", abs_mag=18 * u.mag, phase_param=0.1 * (u.mag / u.deg)) + + # find the reduced mag + red_mag = pc_lin.ReducedMag(alpha) + + assert red_mag.unit == u.mag + assert_array_equal(red_mag.value, np.array([18.0, 19.0])) + + +def test_PhaseCurve_FitModel(): + import numpy as np + import astropy.units as u + from adler.science.PhaseCurve import PhaseCurve + + # define the observations + alpha = np.array([0, 10]) * u.deg + red_mag = np.array([18.0, 19.0]) * u.mag + + # empty linear phase curve model + pc_lin = PhaseCurve(model_name="LinearPhaseFunc") + + # fit the model to the data + pc_fit = pc_lin.FitModel(alpha, red_mag) + + assert pc_fit.H.unit == u.mag + assert pc_fit.H.value == 18.0 + assert pc_fit.S.unit == u.mag / u.deg + assert pc_fit.S.value == 0.1