#!/usr/bin/env python3
"""Reproduce T61 solar context; plain Python, no external dependencies.

NOAA Solar Calculator (Meeus), geometric solar-centre altitude; azimuth
clockwise from true north. Sunrise: centre -0.833 degrees (standard
refraction plus solar radius); civil dawn: centre -6 degrees. Unobstructed
sea-level horizon assumed. Gregorian Julian day, valid before 1900.
https://gml.noaa.gov/grad/solcalc/main.js
https://gml.noaa.gov/grad/solcalc/calcdetails.html
"""
import json
import math
from datetime import datetime, timezone, timedelta
from pathlib import Path

LAT, LON = 41.34, -71.91  # Approximate Stonington context supplied in task brief.
DAY = datetime(1892, 5, 14, tzinfo=timezone.utc)
D = math.degrees
R = math.radians

def solar(utc_minutes):
    instant = DAY + timedelta(minutes=utc_minutes)
    jd = 2440587.5 + (instant - datetime(1970, 1, 1, tzinfo=timezone.utc)).total_seconds()/86400
    t = (jd - 2451545)/36525
    l0 = (280.46646 + t*(36000.76983 + t*0.0003032)) % 360
    m = 357.52911 + t*(35999.05029 - 0.0001537*t)
    e = 0.016708634 - t*(0.000042037 + 0.0000001267*t)
    c = (math.sin(R(m))*(1.914602-t*(0.004817+0.000014*t))
         + math.sin(R(2*m))*(0.019993-0.000101*t) + math.sin(R(3*m))*0.000289)
    omega = 125.04 - 1934.136*t
    apparent_long = l0+c-0.00569-0.00478*math.sin(R(omega))
    obliq = 23+(26+(21.448-t*(46.815+t*(0.00059-t*0.001813)))/60)/60
    obliq += 0.00256*math.cos(R(omega))
    decl = math.asin(math.sin(R(obliq))*math.sin(R(apparent_long)))
    y = math.tan(R(obliq)/2)**2
    eq = 4*D(y*math.sin(R(2*l0))-2*e*math.sin(R(m))
             +4*e*y*math.sin(R(m))*math.cos(R(2*l0))
             -0.5*y*y*math.sin(R(4*l0))-1.25*e*e*math.sin(R(2*m)))
    ha = R((utc_minutes+eq+4*LON)/4-180)
    alt = D(math.asin(math.sin(R(LAT))*math.sin(decl)+math.cos(R(LAT))*math.cos(decl)*math.cos(ha)))
    az = (D(math.atan2(math.sin(ha),math.cos(ha)*math.sin(R(LAT))-math.tan(decl)*math.cos(R(LAT))))+180)%360
    return {'altitude_deg':alt,'azimuth_deg':az,'utc':instant.isoformat()}

def crossing(altitude):
    lo,hi = 6*60,12*60
    for _ in range(50):
        mid=(lo+hi)/2
        if solar(mid)['altitude_deg']<altitude:lo=mid
        else:hi=mid
    return (lo+hi)/2

def clock(minutes):
    return (DAY+timedelta(minutes=minutes)).strftime('%H:%M:%S')

def distance(lat,lon):
    a=math.sin(R(lat-LAT)/2)**2+math.cos(R(LAT))*math.cos(R(lat))*math.sin(R(lon-LON)/2)**2
    return 6371.0088*2*math.asin(math.sqrt(a))

def calculate():
    result={'date':'1892-05-14','coordinate_assumption':{'latitude':LAT,'longitude':LON,'source':'User brief; approximate town location, not an easel position'},
            'method':'NOAA Solar Calculator / Meeus, Python transcription; no atmospheric refraction in reported altitude',
            'method_urls':['https://gml.noaa.gov/grad/solcalc/main.js','https://gml.noaa.gov/grad/solcalc/calcdetails.html'],
            'local_mean_time_ahead_of_EST_minutes':300+4*LON,'clocks':{}}
    dawn,rise=crossing(-6),crossing(-0.833)
    for name,offset in [('Eastern Standard Time',-300),('Stonington local mean time',4*LON)]:
        at_five=300-offset
        result['clocks'][name]={'offset_from_UT_minutes':offset,'at_0500':solar(at_five),
              'civil_dawn':clock(dawn+offset),'sunrise':clock(rise+offset),
              'minutes_after_sunrise_at_0500':at_five-rise,
              'map_0800_EST_on_this_clock':clock(780+offset),'map_delay_minutes':780-at_five}
    result['block_island_distance']={'km':distance(41.1667,-71.5667),'station_coordinate_proxy':[41.1667,-71.5667],
        'station_id':'USW00014799','source_url':'https://www.ncei.noaa.gov/pub/data/ghcn/daily/ghcnd-stations.txt',
        'method':'Great-circle haversine, mean Earth radius 6371.0088 km; later station coordinates used only as an island-location proxy, not the verified 1892 instrument position.'}
    return result

if __name__=='__main__':
    output=Path(__file__).with_name('solar.json')
    output.write_text(json.dumps(calculate(),indent=2)+'\n')
    print(output.read_text())