# Requires Python, pyswisseph and scipy. Swiss Ephemeris AGPL or commercial licensing applies.
import gzip,json,math,datetime as dt
from pathlib import Path
import swisseph as s
from scipy.optimize import brentq,minimize_scalar
O=Path(__file__).resolve().parent; data=json.load(gzip.open(O/'usa-pair.json.gz','rt'));start=s.julday(1776,1,1,12);events=[]
def calc(j):
 p=[]
 for b in [s.JUPITER,s.VENUS]:
  x,f=s.calc_ut(j,b,s.FLG_MOSEPH|s.FLG_SPEED);assert f&s.FLG_MOSEPH;p.append(x)
 return p
def sep(j):
 a,b=calc(j);r=math.pi/180
 return math.degrees(math.acos(max(-1,min(1,math.sin(a[1]*r)*math.sin(b[1]*r)+math.cos(a[1]*r)*math.cos(b[1]*r)*math.cos((a[0]-b[0])*r)))))
def dif(j):
 a,b=calc(j);return (a[0]-b[0]+180)%360-180
def iso(j):
 y,m,d,h=s.revjul(j);t=dt.datetime(y,m,d,tzinfo=dt.timezone.utc)+dt.timedelta(hours=h);return t.isoformat(timespec='seconds').replace('+00:00','Z')
for i,(a,b) in enumerate(zip(data['rows'],data['rows'][1:])):
 x=(a[0]-a[3]+180)%360-180;y=(b[0]-b[3]+180)%360-180
 if x*y<0 and abs(y-x)<180:
  j=brentq(dif,start+i,start+i+1,xtol=1e-9);opt=minimize_scalar(lambda t:sep(j+t),bounds=(-3,3),method='bounded',options={'xatol':1e-9});p=calc(j)
  events.append({'utc':iso(j),'longitude':p[0][0],'fullSeparation':sep(j),'minimumUtc':iso(j+opt.x),'minimumSeparation':float(opt.fun)})
json.dump({'meta':{'method':'Every consecutive noon pair brackets wrapped equal-longitude roots; reject 180° discontinuities, refine with Brent. Minimize spherical separation separately within ±3 days of each root. Swiss Ephemeris explicit Moshier. UTC dates; historical event hours not asserted.','start':'1776-01-01','end':'2027-12-31'},'events':events},open(O/'usa-conjunctions.json','w'),indent=2)
print(len(events),'searched conjunctions');print([e for e in events if e['utc'].startswith('1776')])