1 files changed, 182 insertions, 0 deletions
diff --git a/src/astro.cc b/src/astro.cc
new file mode 100644
index 0000000..17cd463
--- /dev/null
+++ b/src/astro.cc
@@ -0,0 +1,182 @@
+// Copyright 2022 Benjamin Barenblat
+//
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not
+// use this file except in compliance with the License. You may obtain a copy of
+// the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+// License for the specific language governing permissions and limitations under
+// the License.
+
+#include "src/astro.h"
+
+#include <math.h>
+
+#include <chrono>
+#include <utility>
+
+#include "third_party/date/include/date/tz.h"
+
+namespace glplanet {
+
+namespace {
+
+// The J2000 epoch--January 1, 2000, at 12:00:00 Terrestrial Time. This is
+// January 1, 2000, at 11:59:27.816 TAI (because TAI lags TT(TAI) by 32.184 s),
+// or January 1, 2000, at 11:58:55.816 UTC (because TAI-UTC was 32 seconds on
+// January 1, 2000).
+constexpr auto kJ2000Epoch =
+    std::chrono::time_point<date::tai_clock, std::chrono::milliseconds>() +
+    std::chrono::days(42 * 365 + 10 /* leap days */) + std::chrono::hours(11) +
+    std::chrono::minutes(59) + std::chrono::milliseconds(27'816);
+
+constexpr double kMillisecondsPerEphemerisCentury =
+    uint64_t{36'525} * 24 * 60 * 60 * 1'000;
+
+constexpr double kRadiansPerDegree = M_PI / 180.0;
+
+struct CommonValues {
+  // The time, measured in ephemeris centuries from the J2000 epoch.
+  double t;
+
+  double sun_mean_longitude_degrees;
+  double omega;
+  double ecliptic_obliquity;
+};
+
+CommonValues ComputeCommonValues(
+    std::chrono::time_point<date::tai_clock, std::chrono::milliseconds>
+        now) noexcept {
+  const double t =
+      static_cast<double>(std::chrono::duration_cast<std::chrono::milliseconds>(
+                              now - kJ2000Epoch)
+                              .count()) /
+      kMillisecondsPerEphemerisCentury;
+
+  // Digit groupings in these magic constants match those in Meeus for easy
+  // auditing.
+  const double mean_longitude_degrees =
+      (0.000'3032 * t + 36'000.769'83) * t + 280.46646;
+  const double omega =
+      (((2.222222222222222e-6 * t + 0.002'0708) * t + -1934.136'261) * t +
+       125.04452) *
+      kRadiansPerDegree;
+
+  // From Meeus (22.2) and (25.8).
+  const double ecliptic_obliquity_degrees =
+      ((5.036111111111111e-7 * t + -1.6388888888888888e-7) * t +
+       -1.3004166666666667e-2) *
+          t +
+      0.00256 * cos(omega) +  // (25.8)
+      23.43929111111111;
+
+  return {
+      .t = t,
+      .sun_mean_longitude_degrees = mean_longitude_degrees,
+      .omega = omega,
+      .ecliptic_obliquity = ecliptic_obliquity_degrees * kRadiansPerDegree,
+  };
+}
+
+std::pair<double, double> SunEquatorialPositionWithCommonValues(
+    const CommonValues& common) noexcept {
+  // Reference: Jean Meeus, _Astronomical Algorithms_, 2nd edition,
+  // Willmann-Bell, 1998 (ISBN 0-943396-61-1), chapter 25.
+
+  const auto& [t, mean_longitude_degrees, omega, ecliptic_obliquity] = common;
+
+  // Again, digit groupings match Meeus.
+  const double mean_anomaly_degrees =
+      (0.000'1537 * t + 35'999.050'29) * t + 357.52911;
+  const double equation_of_center_degrees =
+      0.000'289 * sin(mean_anomaly_degrees * (M_PI / 60.0)) +
+      (-0.000'101 * t + 0.019'993) * sin(mean_anomaly_degrees * (M_PI / 90.0)) +
+      ((-0.000'014 * t - 0.004'817) * t + 1.914'602) *
+          sin(mean_anomaly_degrees * (M_PI / 180.0));
+  const double apparent_longitude_degrees = -0.00478 * sin(omega) - 0.00569 +
+                                            equation_of_center_degrees +
+                                            mean_longitude_degrees;
+
+  const double apparent_longitude =
+      apparent_longitude_degrees * kRadiansPerDegree;
+
+  // Per Meeus, latitude never exceeds 1.2 arcseconds = 0.0003 degrees,
+  // within margin of error. We therefore have:
+  const double right_ascension_radians =
+      atan2(cos(ecliptic_obliquity) * sin(apparent_longitude),
+            cos(apparent_longitude));
+  const double declination_radians =
+      asin(sin(ecliptic_obliquity) * sin(apparent_longitude));
+
+  return {right_ascension_radians, declination_radians};
+}
+
+double ApparentSiderealTimeAtGreenwichWithCommonValues(
+    const CommonValues& common) noexcept {
+  // Reference: Meeus, page 88.
+
+  const auto& [t, mean_longitude_degrees, omega, ecliptic_obliquity] = common;
+
+  const double moon_mean_longitude_degrees = 481'267.8831 * t + 218.3165;
+  const double longitude_nutation_degrees =
+      5.833333333333333e-5 * sin(2 * omega) +
+      -6.38888888888888e-5 * sin(moon_mean_longitude_degrees * M_PI / 90.0) +
+      -3.6666666666666666e-4 * sin(mean_longitude_degrees * M_PI / 90.0) +
+      -4.7777777777777777e-3 * sin(omega);
+  double greenwich_apparent_sidereal_time_degrees =
+      ((-2.5833118057349522e-8 * t + 0.000'387'933) * t + 13185000.770053742) *
+          t +
+      longitude_nutation_degrees * cos(ecliptic_obliquity) +
+      280.460'618'37;  // (12.4)
+  greenwich_apparent_sidereal_time_degrees =
+      fmod(greenwich_apparent_sidereal_time_degrees, 360.0);
+  if (greenwich_apparent_sidereal_time_degrees < -180.0) {
+    greenwich_apparent_sidereal_time_degrees += 360.0;
+  }
+
+  return greenwich_apparent_sidereal_time_degrees * kRadiansPerDegree;
+}
+
+}  // namespace
+
+std::pair<double, double> SunEquatorialPosition(
+    const std::chrono::time_point<date::tai_clock, std::chrono::milliseconds>
+        now) noexcept {
+  return SunEquatorialPositionWithCommonValues(ComputeCommonValues(now));
+}
+
+double ApparentSiderealTimeAtGreenwich(
+    std::chrono::time_point<date::tai_clock, std::chrono::milliseconds>
+        now) noexcept {
+  return ApparentSiderealTimeAtGreenwichWithCommonValues(
+      ComputeCommonValues(now));
+}
+
+std::pair<double, double> HighNoonLocation(
+    std::chrono::time_point<date::tai_clock, std::chrono::milliseconds>
+        now) noexcept {
+  CommonValues common = ComputeCommonValues(now);
+  const auto& [right_ascension, declination] =
+      SunEquatorialPositionWithCommonValues(common);
+  double greenwich_apparent_sidereal_time =
+      ApparentSiderealTimeAtGreenwichWithCommonValues(common);
+
+  // The declination is the latitude.
+  const double latitude = declination;
+
+  // Computing the longitude is only a bit more complex. We're sitting right
+  // under the sun, so our local hour angle is 0, giving
+  //
+  //   east longitude = right ascension - apparent sidereal time at Greenwich
+  //
+  // (cf. Meeus, page 92).
+  const double longitude = right_ascension - greenwich_apparent_sidereal_time;
+
+  return {longitude, latitude};
+}
+
+}  // namespace glplanet