Varenius/more coordinates (#39)

* Replace 'now' with 'when' internally

* Added calculations for Sun position in ecliptic, equatorial and horizontal coordinates + test

* Added ecliptic-->equatorial conversion

common/src/coords.rs

@@ -2,21 +2,21 @@ use crate::{Direction, Location};
 use chrono::prelude::*;
 use std::f64::consts::PI;
 
-fn julian_day(now: DateTime<Utc>) -> f64 {
-    // Calculate decimal julian day for current date. We can simplify
+fn julian_day(when: DateTime<Utc>) -> f64 {
+    // Calculate decimal julian day for specified date. We can simplify
     // since we do not need to cover dates in the past, only the future!
     // From https://aa.usno.navy.mil/data/JulianDate we get that for
     // A.D. 2000 January 1 	12:00:00.0 correspond to julian day 2451545.0.
     // Calculate difference to this date
     let jdref = Utc.with_ymd_and_hms(2000, 1, 1, 12, 0, 0).unwrap();
-    let diff = now.signed_duration_since(jdref);
+    let diff = when.signed_duration_since(jdref);
     // Need f64 for precision
     2451545.0 + (diff.num_milliseconds() as f64 / (24.0 * 60.0 * 60.0 * 1000.0))
 }
 
-fn gmst(now: DateTime<Utc>) -> f64 {
+fn gmst(when: DateTime<Utc>) -> f64 {
     // Algoritm from https://aa.usno.navy.mil/faq/GAST
-    let jd = julian_day(now);
+    let jd = julian_day(when);
     let jd0 = jd.floor() + 0.5;
     let h = (jd - jd0) * 24.0;
     let dtt = jd - 2451545.0;
@@ -90,6 +90,50 @@ pub fn horizontal_from_galactic(
     horizontal_from_equatorial(location, when, ra, dec)
 }
 
+fn equatorial_from_ecliptic(l: f64, b: f64) -> (f64, f64) {
+    // From https://astrophysicsandpython.com/2021/04/19/celestial-coordinate-conversions/#Equatorial_Ecliptic
+    let _ecliptic = 23.43927944_f64.to_radians();
+    let dec = (b.sin() * _ecliptic.cos() + b.cos() * _ecliptic.sin() * l.sin()).asin();
+    let ra =
+        (b.cos() * l.sin() * _ecliptic.cos() - b.sin() * _ecliptic.sin()).atan2(b.cos() * l.cos());
+    (ra, dec)
+}
+
+fn ecliptic_from_sun(when: DateTime<Utc>) -> (f64, f64) {
+    // Algorithm from https://aa.usno.navy.mil/faq/sun_approx
+    // for computing the Sun's angular coordinates to an accuracy of about 1 arcminute within two centuries of 2000
+    let d = julian_day(when) - 2451545.0;
+    //Mean anomaly of the Sun:
+    let g = (357.529 + 0.98560028 * d) % 360.0;
+    //Mean longitude of the Sun:
+    let q = (280.459 + 0.98564736 * d) % 360.0;
+    // Geocentric apparent ecliptic longitude of the Sun (adjusted for aberration):
+    let l = (q + 1.915 * g.to_radians().sin() + 0.020 * (2.0 * g).to_radians().sin()) % 360.0;
+    // where all the constants (therefore g, q, and L) are in degrees.
+    // It may be necessary or desirable to reduce g, q, and L to the range 0° to 360°.
+    //The Sun's ecliptic latitude, b, can be approximated by b=0.
+    let b = 0.0;
+    // return in radians
+    (l.to_radians(), b)
+}
+
+fn equatorial_from_sun(when: DateTime<Utc>) -> (f64, f64) {
+    // Algorithm from https://aa.usno.navy.mil/faq/sun_approx
+    // for computing the Sun's angular coordinates to an accuracy of about 1 arcminute within two centuries of 2000
+    let (l, _b) = ecliptic_from_sun(when);
+    let d = julian_day(when) - 2451545.0;
+    //First compute the mean obliquity of the ecliptic:
+    let e = (23.439 - 0.00000036 * d).to_radians();
+    let ra = (e.cos() * l.sin()).atan2(l.cos());
+    let dec = (e.sin() * l.sin()).asin();
+    (ra, dec)
+}
+
+pub fn horizontal_from_sun(location: Location, when: DateTime<Utc>) -> Direction {
+    let (ra, dec) = equatorial_from_sun(when);
+    horizontal_from_equatorial(location, when, ra, dec)
+}
+
 #[cfg(test)]
 mod test {
     use chrono::Duration;
@@ -124,4 +168,21 @@ mod test {
             1e-6
         );
     }
+    #[test]
+    fn test_horizontal_from_sun() {
+        // Test that we get the correct horizontal position for the Sun
+        // given specific location and time
+        let jdref = Utc.with_ymd_and_hms(2023, 4, 4, 12, 0, 0).unwrap();
+        // Use SALSA Onsala location
+        let locref = Location {
+            longitude: 0.20802143022,
+            latitude: 1.00170457462,
+        };
+        let dir = horizontal_from_sun(locref, jdref);
+        // Expected horizontal coordinates in radians
+        let expected_az = 3.386904823113701;
+        let expected_alt = 0.6557470215389855;
+        assert_similar!(dir.azimuth, expected_az, 1e-6);
+        assert_similar!(dir.altitude, expected_alt, 1e-6);
+    }
 }