From 73c600be7832640799b3dded4cc69468956e267a Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Sat, 23 Nov 2024 11:04:49 -0800 Subject: [PATCH] Comet phase correction --- Cargo.toml | 2 +- src/kete/rust/flux/common.rs | 22 +++++++------ src/kete_core/Cargo.toml | 2 +- src/kete_core/src/flux/comets.rs | 56 +++++++++++++++++++++++++++----- 4 files changed, 62 insertions(+), 20 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index ff7d1c3..6fecc81 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,7 +9,7 @@ keywords = ["physics", "simulation", "astronomy"] [dependencies] kete_core = { version = "*", path = "src/kete_core", features=["pyo3", "polars"]} -pyo3 = { version = "^0.22.1", features = ["extension-module"] } +pyo3 = { version = "^0.22.6", features = ["extension-module"] } serde = { version = "^1.0.203", features = ["derive"] } nalgebra = {version = "^0.33.0"} itertools = "^0.13.0" diff --git a/src/kete/rust/flux/common.rs b/src/kete/rust/flux/common.rs index 884244c..65802bf 100644 --- a/src/kete/rust/flux/common.rs +++ b/src/kete/rust/flux/common.rs @@ -333,8 +333,8 @@ pub fn frm_thermal_py( /// /// The model for apparent magnitudes are: /// -/// m1 + k1 * log10(sun2obj.r) + 5.0 * log10(obj2obs.r) + phase_mag_slope * phase -/// m2 + k2 * log10(sun2obj.r) + 5.0 * log10(obj2obs.r) + phase_mag_slope * phase +/// m1 + k1 * log10(sun2obj.r) + 5.0 * log10(obj2obs.r) + phase_mag_slope_1 * phase +/// m2 + k2 * log10(sun2obj.r) + 5.0 * log10(obj2obs.r) + phase_mag_slope_2 * phase /// /// Where m1/k1 are related to total magnitudes and m2/k2 are nucleus magnitudes. /// @@ -346,6 +346,9 @@ pub fn frm_thermal_py( /// which are present in the wikipedia definition, this matches the definitions used by /// JPL Horizons. /// +/// Note that this includes seperate phase magnitude slope for both the nucleus and +/// total fluxes. +/// /// This does a best effort to compute both magnitudes, if any values are missing this /// will return None in the respective calculation. /// @@ -365,29 +368,28 @@ pub fn frm_thermal_py( /// and K_2 nucleus magnitude slope as a function of heliocentric distance. /// phase_corr : /// Magnitude variation of the comet as a function of observing phase, units are -/// Mag/Deg of phase, this defaults to 0.035 Mag/Deg. +/// Mag/Deg of phase, this defaults to 0.0 for the total magnitude, and 0.035 +/// Mag/Deg for the nucleus. /// /// Returns /// ------- /// float /// (Total apparent magnitude, Magnitude of the nucleus) -#[allow(clippy::too_many_arguments)] #[pyfunction] -#[pyo3(name = "comet_apparent_mags", signature = (sun2obj, sun2obs, mk_1=None, mk_2=None, phase_corr=None))] +#[pyo3(name = "comet_apparent_mags", signature = (sun2obj, sun2obs, mk_1=None, mk_2=None, phase_corr=[0.0, 0.035]))] pub fn comet_mags_py( sun2obj: VectorLike, sun2obs: VectorLike, mk_1: Option<[f64; 2]>, mk_2: Option<[f64; 2]>, - phase_corr: Option, + phase_corr: [f64; 2], ) -> (Option, Option) { let sun2obj = sun2obj.into_vec(PyFrames::Ecliptic); let sun2obs = sun2obs.into_vec(PyFrames::Ecliptic); - let corr = phase_corr.unwrap_or(0.035); - let mk_params = CometMKParams::new("".into(), mk_1, mk_2, corr); + let mk_params = CometMKParams::new("".into(), mk_1, mk_2, phase_corr); ( - mk_params.apparent_total_flux(&sun2obs, &sun2obj), - mk_params.apparent_nuclear_flux(&sun2obs, &sun2obj), + mk_params.apparent_total_mag(&sun2obs, &sun2obj), + mk_params.apparent_nuclear_mag(&sun2obs, &sun2obj), ) } diff --git a/src/kete_core/Cargo.toml b/src/kete_core/Cargo.toml index def24b5..1c57c53 100644 --- a/src/kete_core/Cargo.toml +++ b/src/kete_core/Cargo.toml @@ -20,7 +20,7 @@ lazy_static = "^1.5.0" nalgebra = {version = "^0.33.0"} polars = {version = "0.43.1", optional=true, features=["parquet", "polars-io"]} pathfinding = "^4.10.0" -pyo3 = { version = "^0.22.1", optional=true} +pyo3 = { version = "^0.22.6", optional=true} rayon = "^1.10.0" serde = { version = "^1.0.203", features = ["derive"] } sgp4 = "^2.2.0" diff --git a/src/kete_core/src/flux/comets.rs b/src/kete_core/src/flux/comets.rs index 4f39364..89ca9f6 100644 --- a/src/kete_core/src/flux/comets.rs +++ b/src/kete_core/src/flux/comets.rs @@ -5,7 +5,23 @@ use serde::{Deserialize, Serialize}; /// /// /// +/// The model for apparent magnitudes are: +/// +/// m1 + k1 * log10(sun2obj.r) + 5.0 * log10(obj2obs.r) + phase_mag_slope_1 * phase +/// m2 + k2 * log10(sun2obj.r) + 5.0 * log10(obj2obs.r) + phase_mag_slope_2 * phase +/// +/// Where m1/k1 are related to total magnitudes and m2/k2 are nucleus magnitudes. +/// +/// This model is based off of these: +/// https://ssd.jpl.nasa.gov/horizons/manual.html#obsquan (see section 9) +/// https://en.wikipedia.org/wiki/Absolute_magnitude#Cometary_magnitudes +/// +/// Note that the above model does not include a 2.5x term attached to the K1/2 terms +/// which are present in the wikipedia definition, this matches the definitions used by +/// JPL Horizons. +/// /// This model additionally includes a correction for phase effects. +/// #[derive(Debug, Deserialize, Serialize)] pub struct CometMKParams { /// Designation (name) of the object. @@ -17,8 +33,8 @@ pub struct CometMKParams { /// M2 and K2 if defined. pub mk_2: Option<[f64; 2]>, - /// Phase correction coefficient in units of Mag/Deg - pub phase_corr_coef: f64, + /// Phase correction coefficients in units of Mag/Deg + pub phase_corr_coef: [f64; 2], } impl CometMKParams { @@ -27,19 +43,19 @@ impl CometMKParams { desig: String, mk_1: Option<[f64; 2]>, mk_2: Option<[f64; 2]>, - phase_corr: f64, + phase_corr_coef: [f64; 2], ) -> Self { Self { desig, mk_1, mk_2, - phase_corr_coef: phase_corr, + phase_corr_coef, } } /// Compute the apparent total flux including both coma and nucleus of the comet. /// This includes an additional 0.035 Mag/Deg phase correction. - pub fn apparent_total_flux( + pub fn apparent_total_mag( &self, sun2obs: &Vector3, sun2obj: &Vector3, @@ -48,13 +64,13 @@ impl CometMKParams { let obj2obs = -sun2obj + sun2obs; let obs_dist = obj2obs.norm(); let helio_dist = sun2obj.norm(); - let phase_corr = self.phase_corr_coef * obj2obs.angle(&-sun2obj).to_degrees(); + let phase_corr = self.phase_corr_coef[0] * obj2obs.angle(&-sun2obj).to_degrees(); Some(m1 + k1 * helio_dist.log10() + 5.0 * obs_dist.log10() + phase_corr) } /// Compute the apparent nuclear flux of the comet, not including the coma. /// This includes an additional 0.035 Mag/Deg phase correction. - pub fn apparent_nuclear_flux( + pub fn apparent_nuclear_mag( &self, sun2obs: &Vector3, sun2obj: &Vector3, @@ -63,7 +79,31 @@ impl CometMKParams { let obj2obs = -sun2obj + sun2obs; let obs_dist = obj2obs.norm(); let helio_dist = sun2obj.norm(); - let phase_corr = self.phase_corr_coef * obj2obs.angle(&-sun2obj).to_degrees(); + let phase_corr = self.phase_corr_coef[1] * obj2obs.angle(&-sun2obj).to_degrees(); Some(m2 + k2 * helio_dist.log10() + 5.0 * obs_dist.log10() + phase_corr) } } + +mod tests { + use super::*; + + #[test] + fn test_comet_mags() { + // Testing 12P against JPL horizons values. + + // Horizons values: 15.757 19.192 + let mk_1 = [5.0, 15.0]; + let mk_2 = [11.0, 10.0]; + + let pos = Vector3::from([-1.154216937776, -2.385684461737, -1.893975509337]); + let obs = Vector3::from([0.469511752038, 0.868580775506, -5.2896978e-5]); + + let comet_mag = CometMKParams::new("12P".into(), Some(mk_1), Some(mk_2), [0.0, 0.0]); + + let nulc_mags = comet_mag.apparent_nuclear_mag(&obs, &pos).unwrap(); + let total_mags = comet_mag.apparent_total_mag(&obs, &pos).unwrap(); + + assert!((total_mags - 15.757).abs() < 1e-3); + assert!((nulc_mags - 19.192).abs() < 1e-3); + } +}