#!/usr/bin/env python3
# Serpent Mound summer-solstice sunset azimuth check used on truthinstone.com
# Published by Truth in Stone (https://truthinstone.com/claims/serpent-mound-solstice-alignment/).
# License: CC0 1.0 Universal (public domain dedication). https://creativecommons.org/publicdomain/zero/1.0/
# Requires: Python 3, numpy, pyerfa.
"""
Obliquity of the ecliptic from the IAU long-term precession model of Vondrak, Capitaine & Wallace (2011)
(angle between the mean pole of date from erfa.ltp and the ecliptic pole of date from erfa.ltecm).
At the June solstice the Sun's declination equals the obliquity (mean; nutation +/-0.003 deg ignored).
Sunset azimuth for a Sun of declination d seen at altitude h from latitude phi:
    cos(Az_from_north_west_side) -> Az = 360 - arccos((sin d - sin phi sin h)/(cos phi cos h))
h = apparent horizon altitude minus refraction/semidiameter. We report three cases:
   h = -0.833 deg (upper limb on a flat sea-level horizon, standard refraction)
   h = +0.0 and +1.0 deg (Sun centre on a raised local horizon; the true skyline at the site is not surveyed here)
Serpent Mound: 39.0253 N, 83.4303 W.
"""
import math, numpy as np, erfa
LAT = 39.0253
def obliquity(y):
    pole = erfa.ltp(float(y))[2]          # mean celestial pole of date, in J2000 frame
    epole = erfa.ltecm(float(y))[2]       # ecliptic pole of date, in J2000 frame
    return math.degrees(math.acos(np.dot(pole, epole)))
def sunset_az(dec, h, lat=LAT):
    d, p, hh = map(math.radians, (dec, lat, h))
    c = (math.sin(d) - math.sin(p)*math.sin(hh))/(math.cos(p)*math.cos(hh))
    return 360 - math.degrees(math.acos(c))
print("year   obliquity  az(h=-0.833)  az(h=0)  az(h=+1)")
for lab, y in [("13,000 yr ago (~11,000 BC)", -10999), ("Adena, 300 BC", -299), ("AD 1070 (charcoal dates)", 1070), ("AD 1988 (survey)", 1988), ("AD 2026", 2026)]:
    e = obliquity(y)
    print(f"{lab:28s} {e:8.4f}  {sunset_az(e,-0.833):9.2f}  {sunset_az(e,0):8.2f}  {sunset_az(e,1):8.2f}")
# Which epoch would put sunset at 302.0 deg (Hancock's figure) or 300.08 deg (Fletcher & Cameron 1988)?
ys = np.arange(-30000, 2027, 10)
for target in (302.0, 300.0833):
    for h in (-0.833, 0.0, 1.0):
        az = np.array([sunset_az(obliquity(y), h) for y in ys])
        i = int(np.argmin(abs(az - target)))
        print(f"target {target:.2f} h={h:+.3f}: closest az {az[i]:.2f} at year {ys[i]} (max az over 30 kyr {az.max():.2f}, min {az.min():.2f})")
e = obliquity(-10999) - obliquity(2026)
print(f"Change in obliquity 13,000 yr ago vs today: {e:+.3f} deg -> azimuth shift {sunset_az(obliquity(-10999),-0.833)-sunset_az(obliquity(2026),-0.833):+.2f} deg")