#!/usr/bin/env python3
# Teotihuacan Sun Pyramid sunset-date check used on truthinstone.com
# Published by Truth in Stone (https://truthinstone.com/claims/teotihuacan-astronomical-layout/).
# License: CC0 1.0 Universal. Requires: Python 3, numpy, pyerfa.
"""
For the Sun Pyramid's east-west axis (azimuth of the western direction about 285.5 deg, i.e. 15.5 deg north of west;
Sprajc 2000) we find the Sun declination at which the Sun's upper limb sets on that azimuth, then the dates in AD 200
(proleptic Gregorian) when the Sun has that declination at local sunset. Sun position: the standard low-precision
solar theory (Meeus, Astronomical Algorithms, ch. 25) with the long-term obliquity of Vondrak et al. (2011) from erfa.
Its declination error (a few hundredths of a degree) is far below the azimuth uncertainty of the structure.
A check run for AD 2000 gives 29 April and 13 August, the dates usually quoted today.
Latitude 19.692 N, longitude 98.844 W. Sun-centre altitudes -0.833 (upper limb on a flat horizon), 0, +0.5 and +1.0 deg (western skyline not surveyed here).
"""
import math, numpy as np, erfa, datetime
LAT = math.radians(19.692)
def obliq(y):
    return math.acos(np.dot(erfa.ltp(float(y))[2], erfa.ltecm(float(y))[2]))
def sun_dec_at(jd):
    # low-precision solar longitude (Meeus ch. 25), referred to mean equinox of date
    T = (jd - 2451545.0)/36525
    L0 = 280.46646 + 36000.76983*T + 0.0003032*T*T
    M = math.radians(357.52911 + 35999.05029*T - 0.0001537*T*T)
    C = (1.914602 - 0.004817*T)*math.sin(M) + (0.019993 - 0.000101*T)*math.sin(2*M) + 0.000289*math.sin(3*M)
    lam = math.radians((L0 + C) % 360)
    y = 2000 + T*100
    return math.asin(math.sin(obliq(y))*math.sin(lam))
def needed_dec(az_deg, h_deg):
    A, h = math.radians(az_deg), math.radians(h_deg)  # h = altitude of the Sun's centre at the moment of sunset
    return math.asin(math.sin(LAT)*math.sin(h) + math.cos(LAT)*math.cos(h)*math.cos(A))
def gregorian_from_jd(jd):
    Z = int(jd + 0.5); a = int((Z - 1867216.25)/36524.25); A = Z + 1 + a - a//4
    B = A + 1524; C = int((B - 122.1)/365.25); D = int(365.25*C); E = int((B - D)/30.6001)
    day = B - D - int(30.6001*E); month = E - 1 if E < 14 else E - 13; year = C - 4716 if month > 2 else C - 4715
    return year, month, day
def jd_greg(y, m, d):
    a = (14 - m)//12; yy = y + 4800 - a; mm = m + 12*a - 3
    return d + (153*mm + 2)//5 + 365*yy + yy//4 - yy//100 + yy//400 - 32045 - 0.5
for AZ in (285.47, 285.5, 286.5):
  for h in (-0.833, 0.0, 0.5, 1.0):
    target = needed_dec(AZ, h)
    base = jd_greg(200, 1, 1) + 1.0  # local sunset ~ 0h UT next day at 98.8 W; use +1.0 d (~ 24h UT)
    days = [base + k for k in range(366)]
    decs = [sun_dec_at(j) for j in days]
    hits = []
    for i in range(len(days)-1):
        if (decs[i]-target)*(decs[i+1]-target) <= 0:
            f = (target-decs[i])/(decs[i+1]-decs[i]); hits.append(days[i] + f)
    ds = [gregorian_from_jd(j) for j in hits]
    gap = hits[1]-hits[0] if len(hits) > 1 else float('nan')
    print(f"az {AZ:.2f} skyline {h:+.3f}: Sun dec {math.degrees(target):6.2f} -> dates AD 200 (Gregorian) {ds}; interval spring->autumn {gap:.1f} d, autumn->spring {365.2422-gap:.1f} d")