Moon

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).
 ///
@@ -128,15 +144,12 @@ 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)
   }
 }
 
@@ -148,3 +161,120 @@ fn test_pos(){
   assert_eq!(0.6412750628729547, pos.azimuth);
   assert_eq!(-0.7000406838781611, pos.altitude);
 }
+
+// Moon
+
+/// Holds illuminated _fraction_ of the moon, the _phase_, and the _angle_.
+#[derive(Debug)]
+pub struct Illumination {
+    pub fraction: f64,
+    pub phase: f64,
+    pub angle: f64,
+}
+
+// general moon calculations
+
+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()
+}
+
+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 {
+    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();
+
+    Coords {
+        right_ascension: right_ascension(lng, lat),
+        declination: declination(lng, lat),
+    }
+}
+
+/// 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;
+    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);
+    const EPSILON : f64 = 1.0e-15;
+
+    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
+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);
+    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
+    };
+
+    Illumination {
+        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);
+    const EPSILON : f64 = 1.0e-15;
+
+    assert!((moon_illum.fraction - 0.4848068202456373).abs() < EPSILON);
+    assert!((moon_illum.phase - 0.7548368838538762).abs() < EPSILON);
+    assert!((moon_illum.angle - 1.6732942678578346).abs() < EPSILON);
+}