// 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