"""AGPL-3.0-or-later. Requires pyswisseph and scipy.
Modern samples use Moshier. Ancient calculations additionally require
DE441 compressed Swiss files seplm06.se1 and semom06.se1 in ephe/ beside this script,
from Astrodienst’s official ephemeris distribution. No network calls are made.
"""
from pathlib import Path
import swisseph as s,json,gzip,math,datetime as dt
from scipy.optimize import brentq,minimize_scalar
O=Path(__file__).resolve().parent;B=[('Sun',s.SUN),('Moon',s.MOON),('Mercury',s.MERCURY),('Venus',s.VENUS),('Mars',s.MARS),('Jupiter',s.JUPITER),('Saturn',s.SATURN)]
wrap=lambda x:x%360
signed=lambda x:(x+180)%360-180
start=dt.datetime(2025,12,30,tzinfo=dt.timezone.utc);end=dt.datetime(2028,1,3,tzinfo=dt.timezone.utc);j0=s.julday(2025,12,30,0);rows=[]
for hour in range(int((end-start).total_seconds()/3600)+1):
 j=j0+hour/24;r=[]
 for name,code in B:
  p,flag=s.calc_ut(j,code,s.FLG_MOSEPH|s.FLG_SPEED);assert flag&s.FLG_MOSEPH;r.extend(round(v,7) for v in [p[0],p[1],p[3]])
 r.append(round(s.pheno_ut(j,s.MOON,s.FLG_MOSEPH)[1]*100,5));rows.append(r)
def modern(j,code):return s.calc_ut(j,code,s.FLG_MOSEPH|s.FLG_SPEED)[0]
def sep(a,b):
 la,lb,dl=map(math.radians,[a[1],b[1],a[0]-b[0]]);return math.degrees(math.acos(max(-1,min(1,math.sin(la)*math.sin(lb)+math.cos(la)*math.cos(lb)*math.cos(dl)))))
def stamp(j):
 y,m,d,h=s.revjul(j,s.GREG_CAL);return (dt.datetime(y,m,d,tzinfo=dt.timezone.utc)+dt.timedelta(hours=h)).isoformat(timespec='seconds').replace('+00:00','Z')
events=[]
for angle in [0,60,90,120,180]:
 for target in ([angle] if angle in [0,180] else [angle,-angle]):
  def f(j):return signed(modern(j,s.VENUS)[0]-modern(j,s.JUPITER)[0]-target)
  a=s.julday(2026,1,1,0);fa=f(a)
  while a<s.julday(2028,1,1,0):
   b=a+.5;fb=f(b)
   if fa*fb<=0 and abs(fb-fa)<180:
    root=brentq(f,a,b,xtol=1e-9);p=modern(root,s.VENUS);q=modern(root,s.JUPITER);e={'angle':angle,'utc':stamp(root),'longitudeSeparation':abs(signed(p[0]-q[0])),'fullSeparation':round(sep(p,q),6),'VenusLongitude':p[0],'JupiterLongitude':q[0]}
    if angle==0:
     minimum=minimize_scalar(lambda offset:sep(modern(root+offset,s.VENUS),modern(root+offset,s.JUPITER)),bounds=(-3,3),method='bounded',options={'xatol':1e-10});e['minimumUtc']=stamp(root+minimum.x);e['minimumSeparation']=round(minimum.fun,6)
    events.append(e)
   a=b;fa=fb
meta={'start':start.isoformat().replace('+00:00','Z'),'end':end.isoformat().replace('+00:00','Z'),'bodies':[n for n,_ in B],'columns':'Each body: apparent geocentric tropical longitude, latitude, longitude speed (degrees/day). Final column: lunar illumination percent.','stepHours':1,'method':'Swiss Ephemeris 2.10.03, explicit Moshier analytical mode. Hourly data linearly interpolated with longitude wrap. UT input; built-in Delta T.','events':'Half-day brackets for signed Venus−Jupiter longitude minus 0, ±60, ±90, ±120, 180 degrees; reject wrap discontinuities; Brent roots. Great-circle conjunction minima separately minimized within ±3 days.','source':'https://www.astro.com/swisseph/swephprg.htm','license':'Swiss Ephemeris © Astrodienst AG; AGPL/professional dual license. This generator is distributed under AGPL-3.0-or-later.'}
with gzip.open(O/'sky.json.gz','wt') as f:json.dump({'meta':meta,'rows':rows,'events':sorted(events,key=lambda x:x['utc'])},f,separators=(',',':'))
# Ancient snapshots and the requested 536 BCE search use the existing DE441 compressed Swiss files.
s.set_ephe_path(str((O/'ephe').resolve()));s.set_tid_acc(s.TIDAL_DE441)
def ancient(j,code):
 p,flag=s.calc(j,code,s.FLG_SWIEPH|s.FLG_SPEED);assert flag&s.FLG_SWIEPH and not flag&s.FLG_MOSEPH;assert s.get_current_file_data(0)[3]==441;return p
def oldstamp(j):
 y,m,d,h=s.revjul(j,s.JUL_CAL);hh=int(h);mm=int((h-hh)*60);return f'{1-y} BCE-{m:02}-{d:02} {hh:02}:{mm:02} TT'
ancient_data={'meta':{'ephemeris':'DE441 compressed Swiss files; runtime file header asserted 441.','calendar':'Proleptic Julian; calculations in TT. Astronomical year −535 = 536 BCE, −567 = 568 BCE.','snapshot':'First day of each selected month, 12:00 TT; no claim of an observed Babylonian civil date.','deltaT':'Ancient Earth rotation is uncertain. Built-in ΔT is recorded; TT angular geometry is more secure than an exact ancient local clock or horizon reconstruction.','observer':'Babylon 32.54 N, 44.42 E; angular positions here are geocentric.','search':'0.5-day brackets of wrapped Venus−Jupiter longitude with discontinuity rejection, Brent refinement. Separately minimize spherical separation within ±3 days.','source':'https://www.astro.com/swisseph/swephprg.htm'},'snapshots':{},'conjunctions536':[]}
for bce in [568,536]:
 samples=[]
 for month in range(1,13):
  j=s.julday(1-bce,month,1,12,s.JUL_CAL);p={n:ancient(j,c) for n,c in B};samples.append({'label':oldstamp(j),'month':month,'positions':{n:{'longitude':v[0],'latitude':v[1],'speed':v[3]} for n,v in p.items()},'separation':sep(p['Venus'],p['Jupiter']),'venusElongation':sep(p['Venus'],p['Sun']),'jupiterElongation':sep(p['Jupiter'],p['Sun']),'deltaTSeconds':s.deltat_ex(j,s.FLG_SWIEPH)*86400})
 ancient_data['snapshots'][str(bce)]=samples
a=s.julday(-535,1,1,0,s.JUL_CAL);endj=s.julday(-534,1,1,0,s.JUL_CAL)
def f(j):return signed(ancient(j,s.VENUS)[0]-ancient(j,s.JUPITER)[0])
while a<endj:
 b=a+.5;fa,fb=f(a),f(b)
 if fa*fb<=0 and abs(fb-fa)<180:
  root=brentq(f,a,b,xtol=1e-9);mn=minimize_scalar(lambda offset:sep(ancient(root+offset,s.VENUS),ancient(root+offset,s.JUPITER)),bounds=(-3,3),method='bounded',options={'xatol':1e-10});ancient_data['conjunctions536'].append({'conjunctionTT':oldstamp(root),'conjunctionJDTT':root,'minimumTT':oldstamp(root+mn.x),'minimumSeparation':mn.fun,'venusElongation':sep(ancient(root,s.VENUS),ancient(root,s.SUN))})
 a=b
(O/'ancient-sky.json').write_text(json.dumps(ancient_data,indent=2))
print(len(rows),'modern hourly rows;',len(events),'exact modern aspect contacts; ancient 536 BCE conjunctions:',ancient_data['conjunctions536'])