Comet phase correction

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"

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<f64>,
+    phase_corr: [f64; 2],
 ) -> (Option<f64>, Option<f64>) {
     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),
     )
 }
 

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"

src/kete_core/src/flux/comets.rs

@@ -5,7 +5,23 @@ use serde::{Deserialize, Serialize};
 ///
 /// <https://en.wikipedia.org/wiki/Absolute_magnitude#Cometary_magnitudes>
 ///
+/// 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<f64>,
         sun2obj: &Vector3<f64>,
@@ -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<f64>,
         sun2obj: &Vector3<f64>,
@@ -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);
+    }
+}