From 49fbc4ffef3f57615ed4cb66c03be2200b57ebf4 Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Wed, 31 Oct 2018 15:55:00 -0700 Subject: [PATCH 1/8] Moon Position --- src/lib.rs | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index bb5bcba..3dda66c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -148,3 +148,77 @@ fn test_pos(){ assert_eq!(0.6412750628729547, pos.azimuth); assert_eq!(-0.7000406838781611, pos.altitude); } + + +// general moon calculations + +fn lunar_mean_anomaly(d: f64) -> f64 { + (134.963 + 13.064993 * d).to_radians() +} + +fn lunar_ecliptic_longitude(d: f64) -> f64 { + (218.316 + 13.176396 * d).to_radians() +} + +fn lunar_mean_distance(d: f64) -> f64 { + (93.272 + 13.229350 * d).to_radians() +} + +fn astro_refraction(h: f64) -> f64 { + let hh = if h < 0.0 { + 0.0 + } else { + h + }; + + 0.0002967 / (hh + 0.00312536 / (hh + 0.08901179)).tan() +} + +#[derive(Debug)] +struct Coords { + pub right_ascension: f64, + pub declination: f64, + pub distance: f64, +} + +fn moon_coords(d: f64) -> Coords { + let l = lunar_ecliptic_longitude(d); + let m = lunar_mean_anomaly(d); + let f = lunar_mean_distance(d); + + let lng = l + TO_RAD * 6.289 * m.sin(); + let lat = TO_RAD * 5.128 * f.sin(); + let distance = 385001.0 - 20905.0 * m.cos(); + + Coords { + right_ascension: right_ascension(lng, lat), + declination: declination(lng, lat), + distance: distance, + } +} + +pub fn moon_pos(unixtime_in_ms: i64, lat: f64, lon: f64) -> Position { + let lw = TO_RAD * -lon; + let phi = TO_RAD * lat; + let d = to_days(unixtime_in_ms); + + let c = moon_coords(d); + + let h = sidereal_time(d, lw) - c.right_ascension; + let mut alt = altitude(h, phi, c.declination); + alt = alt + astro_refraction(alt); + + Position { + azimuth: azimuth(h, phi, c.declination), + altitude: alt + } +} + +#[test] +fn test_moon_pos() { + let date = 1362441600000; + let pos = moon_pos(date, 50.5, 30.5); + + assert_eq!(-0.9783999522438225, pos.azimuth - PI); + assert_eq!(0.014551482243892251, pos.altitude); +} From dfec544a669172d1cac9b4663fd466b2ea338e21 Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Thu, 1 Nov 2018 21:18:59 -0700 Subject: [PATCH 2/8] Remove distance field from Coords --- src/lib.rs | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 3dda66c..190244e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -178,7 +178,6 @@ fn astro_refraction(h: f64) -> f64 { struct Coords { pub right_ascension: f64, pub declination: f64, - pub distance: f64, } fn moon_coords(d: f64) -> Coords { @@ -188,12 +187,10 @@ fn moon_coords(d: f64) -> Coords { let lng = l + TO_RAD * 6.289 * m.sin(); let lat = TO_RAD * 5.128 * f.sin(); - let distance = 385001.0 - 20905.0 * m.cos(); Coords { right_ascension: right_ascension(lng, lat), declination: declination(lng, lat), - distance: distance, } } From 0d464fb24fc8cc55611072ec1ccab5e538027b3b Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Thu, 1 Nov 2018 21:20:29 -0700 Subject: [PATCH 3/8] Extract sun_coords fn --- src/lib.rs | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 190244e..eb0b598 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -128,18 +128,25 @@ pub fn pos(unixtime_in_ms: i64, lat: f64, lon: f64) -> Position { let lw = -lon.to_radians(); let phi = lat.to_radians(); let d = to_days(unixtime_in_ms); - let m = solar_mean_anomaly(d); - let l = ecliptic_longitude(m); - let dec = declination(l, 0.0); - let ra = right_ascension(l, 0.0); - let h = sidereal_time(d, lw) - ra; + let coords = sun_coords(d); + let h = sidereal_time(d, lw) - coords.right_ascension; Position { - azimuth : azimuth(h, phi, dec), - altitude : altitude(h, phi, dec) + azimuth : azimuth(h, phi, coords.declination), + altitude : altitude(h, phi, coords.declination) } } +fn sun_coords (d: f64) -> Coords { + let m = solar_mean_anomaly(d); + let l = ecliptic_longitude(m); + + Coords { + right_ascension: right_ascension(l, 0.0), + declination: declination(l, 0.0), + } +} + #[test] fn test_pos(){ // 2013-03-05 UTC From 48df545a57c5ebdd2d4710cb0486573f88203b3a Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Thu, 1 Nov 2018 21:37:00 -0700 Subject: [PATCH 4/8] Moon Illumination --- src/lib.rs | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index eb0b598..bda8be3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -226,3 +226,49 @@ fn test_moon_pos() { assert_eq!(-0.9783999522438225, pos.azimuth - PI); assert_eq!(0.014551482243892251, pos.altitude); } + +#[derive(Debug)] +pub struct Moon { + pub fraction: f64, + pub phase: f64, + pub angle: f64, +} + +pub fn moon_illumination(unixtime_in_ms: i64) -> Moon { + let d = to_days(unixtime_in_ms); + let s = sun_coords(d); + let m = moon_coords(d); + let a = lunar_mean_anomaly(d); + + let distance = 385001.0 - 20905.0 * a.cos(); // distance to the moon in km + + let sdist = 149598000 as f64; + + let phi = (s.declination.sin() * m.declination.sin() + s.declination.cos() * m.declination.cos() * (s.right_ascension - m.right_ascension).cos()).acos(); + + let inc = (sdist * phi.sin()).atan2(distance - sdist * phi.cos()); + let angle = (s.declination.cos() * (s.right_ascension - m.right_ascension).sin()).atan2(s.declination.sin() * m.declination.cos() - s.declination.cos() * (m.declination).sin() * (s.right_ascension - m.right_ascension).cos()); + + + let sign = if angle < 0.0 { + -1.0 + } else { + 1.0 + }; + + Moon { + fraction: (1.0 + inc.cos()) / 2.0, + phase: 0.5 + 0.5 * inc * sign / PI, + angle: angle, + } +} + +#[test] +fn test_moon_illumination() { + let date = 1362441600000; + let moon_illum = moon_illumination(date); + + assert_eq!(moon_illum.fraction, 0.4848068202456373); + assert_eq!(moon_illum.phase, 0.7548368838538762); + assert_eq!(moon_illum.angle, 1.6732942678578346); +} From bf39043d171cdc0ddb931bd10d2e69ad9a4129a1 Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Thu, 1 Nov 2018 21:46:59 -0700 Subject: [PATCH 5/8] Organizing --- src/lib.rs | 67 ++++++++++++++++++++++++++++-------------------------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index bda8be3..3e1573b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -40,6 +40,12 @@ pub struct Position { pub altitude : f64 } +#[derive(Debug)] +struct Coords { + pub right_ascension: f64, + pub declination: f64, +} + fn to_julian(unixtime_in_ms: i64) -> f64 { unixtime_in_ms as f64 / (MILLISECONDS_PER_DAY as f64) - 0.5 + J1970 as f64 @@ -116,6 +122,16 @@ fn ecliptic_longitude(m:f64) -> f64 { m + equation_of_center(m) + PERIHELION_OF_EARTH + PI } +fn sun_coords (d: f64) -> Coords { + let m = solar_mean_anomaly(d); + let l = ecliptic_longitude(m); + + Coords { + right_ascension: right_ascension(l, 0.0), + declination: declination(l, 0.0), + } +} + /// Calculates the sun position for a given date and latitude/longitude. /// The angles are calculated as [radians](https://en.wikipedia.org/wiki/Radian). /// @@ -137,16 +153,6 @@ pub fn pos(unixtime_in_ms: i64, lat: f64, lon: f64) -> Position { } } -fn sun_coords (d: f64) -> Coords { - let m = solar_mean_anomaly(d); - let l = ecliptic_longitude(m); - - Coords { - right_ascension: right_ascension(l, 0.0), - declination: declination(l, 0.0), - } -} - #[test] fn test_pos(){ // 2013-03-05 UTC @@ -156,20 +162,16 @@ fn test_pos(){ assert_eq!(-0.7000406838781611, pos.altitude); } +// Moon -// general moon calculations - -fn lunar_mean_anomaly(d: f64) -> f64 { - (134.963 + 13.064993 * d).to_radians() -} - -fn lunar_ecliptic_longitude(d: f64) -> f64 { - (218.316 + 13.176396 * d).to_radians() +#[derive(Debug)] +pub struct Moon { + pub fraction: f64, + pub phase: f64, + pub angle: f64, } -fn lunar_mean_distance(d: f64) -> f64 { - (93.272 + 13.229350 * d).to_radians() -} +// general moon calculations fn astro_refraction(h: f64) -> f64 { let hh = if h < 0.0 { @@ -181,10 +183,16 @@ fn astro_refraction(h: f64) -> f64 { 0.0002967 / (hh + 0.00312536 / (hh + 0.08901179)).tan() } -#[derive(Debug)] -struct Coords { - pub right_ascension: f64, - pub declination: f64, +fn lunar_mean_anomaly(d: f64) -> f64 { + (134.963 + 13.064993 * d).to_radians() +} + +fn lunar_ecliptic_longitude(d: f64) -> f64 { + (218.316 + 13.176396 * d).to_radians() +} + +fn lunar_mean_distance(d: f64) -> f64 { + (93.272 + 13.229350 * d).to_radians() } fn moon_coords(d: f64) -> Coords { @@ -201,6 +209,7 @@ fn moon_coords(d: f64) -> Coords { } } +/// calculates the moon position for a given date and latitude/longitude pub fn moon_pos(unixtime_in_ms: i64, lat: f64, lon: f64) -> Position { let lw = TO_RAD * -lon; let phi = TO_RAD * lat; @@ -227,13 +236,7 @@ fn test_moon_pos() { assert_eq!(0.014551482243892251, pos.altitude); } -#[derive(Debug)] -pub struct Moon { - pub fraction: f64, - pub phase: f64, - pub angle: f64, -} - +/// calculates the moon illumination, phase, and angle for a given date pub fn moon_illumination(unixtime_in_ms: i64) -> Moon { let d = to_days(unixtime_in_ms); let s = sun_coords(d); From 0e905d5e2d4332a9ed64b90f538bcbbdb66736b8 Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Thu, 1 Nov 2018 22:30:49 -0700 Subject: [PATCH 6/8] Rename struct. Moon -> Illumination --- src/lib.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 3e1573b..ff60eac 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -165,7 +165,7 @@ fn test_pos(){ // Moon #[derive(Debug)] -pub struct Moon { +pub struct Illumination { pub fraction: f64, pub phase: f64, pub angle: f64, @@ -237,7 +237,7 @@ fn test_moon_pos() { } /// calculates the moon illumination, phase, and angle for a given date -pub fn moon_illumination(unixtime_in_ms: i64) -> Moon { +pub fn moon_illumination(unixtime_in_ms: i64) -> Illumination { let d = to_days(unixtime_in_ms); let s = sun_coords(d); let m = moon_coords(d); @@ -259,7 +259,7 @@ pub fn moon_illumination(unixtime_in_ms: i64) -> Moon { 1.0 }; - Moon { + Illumination { fraction: (1.0 + inc.cos()) / 2.0, phase: 0.5 + 0.5 * inc * sign / PI, angle: angle, From 893d1647e12f2d2b7a525fa261ed0987cb1ab113 Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Thu, 1 Nov 2018 22:31:29 -0700 Subject: [PATCH 7/8] Improve doc comments --- src/lib.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index ff60eac..470eb66 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -164,6 +164,7 @@ fn test_pos(){ // Moon +/// Holds illuminated _fraction_ of the moon, the _phase_, and the _angle_. #[derive(Debug)] pub struct Illumination { pub fraction: f64, @@ -209,7 +210,7 @@ fn moon_coords(d: f64) -> Coords { } } -/// calculates the moon position for a given date and latitude/longitude +/// Calculates the moon position for a given date and latitude/longitude pub fn moon_pos(unixtime_in_ms: i64, lat: f64, lon: f64) -> Position { let lw = TO_RAD * -lon; let phi = TO_RAD * lat; @@ -236,7 +237,7 @@ fn test_moon_pos() { assert_eq!(0.014551482243892251, pos.altitude); } -/// calculates the moon illumination, phase, and angle for a given date +/// Calculates the moon illumination, phase, and angle for a given date pub fn moon_illumination(unixtime_in_ms: i64) -> Illumination { let d = to_days(unixtime_in_ms); let s = sun_coords(d); From 410a401964892f425dd5502e888847df09f06359 Mon Sep 17 00:00:00 2001 From: Travis LaDuke Date: Mon, 5 Nov 2018 20:57:36 -0800 Subject: [PATCH 8/8] Make float comparisons approximate --- src/lib.rs | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 470eb66..b96a7ce 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -232,9 +232,10 @@ pub fn moon_pos(unixtime_in_ms: i64, lat: f64, lon: f64) -> Position { fn test_moon_pos() { let date = 1362441600000; let pos = moon_pos(date, 50.5, 30.5); + const EPSILON : f64 = 1.0e-15; - assert_eq!(-0.9783999522438225, pos.azimuth - PI); - assert_eq!(0.014551482243892251, pos.altitude); + assert!((-0.9783999522438226 - pos.azimuth + PI).abs() < EPSILON); + assert!((0.014551482243892251 - pos.altitude).abs() < EPSILON); } /// Calculates the moon illumination, phase, and angle for a given date @@ -271,8 +272,9 @@ pub fn moon_illumination(unixtime_in_ms: i64) -> Illumination { fn test_moon_illumination() { let date = 1362441600000; let moon_illum = moon_illumination(date); + const EPSILON : f64 = 1.0e-15; - assert_eq!(moon_illum.fraction, 0.4848068202456373); - assert_eq!(moon_illum.phase, 0.7548368838538762); - assert_eq!(moon_illum.angle, 1.6732942678578346); + assert!((moon_illum.fraction - 0.4848068202456373).abs() < EPSILON); + assert!((moon_illum.phase - 0.7548368838538762).abs() < EPSILON); + assert!((moon_illum.angle - 1.6732942678578346).abs() < EPSILON); }