"""engine.py — the computations behind
"We Are Holding Our Own — The Wreck of the Edmund Fitzgerald, Read by Jyotisa"
https://ayurastro.com/writings/we-are-holding-our-own-the-wreck-of-the-edmund-fitzgerald-read-by-jyotisha.html

Requires:  pip install pyswisseph      (Swiss Ephemeris; Moshier theory, no data files)
Run:       python engine.py

Sidereal zodiac, Lahiri ayanamsa, true node, whole-sign houses.
All clock times are Eastern Standard Time (UT - 5 h), as in the
U.S. Coast Guard Marine Board report.
"""
import swisseph as swe, json, math
swe.set_sid_mode(swe.SIDM_LAHIRI)
FL = swe.FLG_MOSEPH | swe.FLG_SIDEREAL | swe.FLG_SPEED
G = [('Su',swe.SUN),('Mo',swe.MOON),('Ma',swe.MARS),('Me',swe.MERCURY),('Ju',swe.JUPITER),('Ve',swe.VENUS),('Sa',swe.SATURN),('Ra',swe.TRUE_NODE)]
SIGNS=['Mesha','Vrsabha','Mithuna','Karka','Simha','Kanya','Tula','Vrscika','Dhanus','Makara','Kumbha','Mina']
NAK=['Asvini','Bharani','Krttika','Rohini','Mrgasirsa','Ardra','Punarvasu','Pusya','Aslesa','Magha','P.Phalguni','U.Phalguni','Hasta','Citra','Svati','Visakha','Anuradha','Jyestha','Mula','P.Asadha','U.Asadha','Sravana','Dhanistha','Satabhisaj','P.Bhadrapada','U.Bhadrapada','Revati']
def dms(x):
    d=int(x); m=int((x-d)*60); return f"{d:02d}°{m:02d}′"
def jd(y,mo,d,h,mi,s=0,tz=0):
    return swe.julday(y,mo,d,h+mi/60+s/3600-tz)
def chart(jdut,lat,lon,label=''):
    out={}
    for n,p in G:
        r,_=swe.calc_ut(jdut,p,FL); out[n]=(r[0]%360,r[3])
    out['Ke']=((out['Ra'][0]+180)%360,out['Ra'][1])
    cusps,asc=swe.houses_ex(jdut,lat,lon,b'W',swe.FLG_SIDEREAL)
    out['La']=(asc[0],0)
    out['ayan']=swe.get_ayanamsa_ut(jdut)
    return out
def show(c,label):
    print('==',label,' ayan',dms(c['ayan']))
    ls=int(c['La'][0]//30)
    for k in ['La','Su','Mo','Ma','Me','Ju','Ve','Sa','Ra','Ke']:
        lon,sp=c[k]; s=int(lon//30); n=int(lon//(360/27)); pada=int((lon%(360/27))//(360/108))+1
        print(f"{k} {SIGNS[s]:8s} {dms(lon%30)} {'R' if sp<0 and k not in('Ra','Ke','La') else ' '} H{(s-ls)%12+1:2d} {NAK[n]} p{pada}  ({lon:.3f})")
    t=((c['Mo'][0]-c['Su'][0])%360)/12
    print('tithi index',int(t)+1,'frac',round(t%1,3),'yoga',int(((c['Mo'][0]+c['Su'][0])%360)/(360/27))+1, 'karana idx', int(((c['Mo'][0]-c['Su'][0])%360)/6)+1)

# ---------- almanac limbs, stations, eclipses, risings ----------
EST=-5
def cal(jdut,tz=EST):
    y,m,d,h=swe.revjul(jdut+tz/24.0)
    hh=int(h); mm=int(round((h-hh)*60)); 
    if mm==60: hh+=1; mm=0
    return f"{y}-{m:02d}-{d:02d} {hh:02d}:{mm:02d}"
def lon(jdut,p):
    return swe.calc_ut(jdut,p,FL)[0][0]%360
def spd(jdut,p): return swe.calc_ut(jdut,p,FL)[0][3]
def solve(f,a,b,tol=1e-6):
    fa=f(a)
    for _ in range(80):
        m=(a+b)/2; fm=f(m)
        if (fa<=0)==(fm<=0): a,fa=m,fm
        else: b=m
        if b-a<tol: break
    return (a+b)/2
def ang(x): return (x+180)%360-180
def elong(j): return (lon(j,swe.MOON)-lon(j,swe.SUN))%360
def ysum(j): return (lon(j,swe.MOON)+lon(j,swe.SUN))%360
def crossings(fn,step_deg,j0,j1,dj=0.02):
    out=[]; j=j0; prev=int(fn(j)//step_deg)
    while j<j1:
        j2=j+dj; cur=int(fn(j2)//step_deg)
        if cur!=prev:
            target=(cur*step_deg) if (cur-prev)%int(360/step_deg) in (1,) else cur*step_deg
            t=solve(lambda x: ang(fn(x)-cur*step_deg),j,j2)
            out.append((t,cur)); prev=cur
        j=j2
    return out
KAR=['Bava','Bālava','Kaulava','Taitila','Gara','Vaṇij','Viṣṭi']
def karana_name(k): # k 0..59
    if k==0: return 'Kiṁstughna'
    if k>=57: return ['Śakuni','Catuṣpada','Nāga'][k-57]
    return KAR[(k-1)%7]
YOG=['Viṣkambha','Prīti','Āyuṣmān','Saubhāgya','Śobhana','Atigaṇḍa','Sukarmā','Dhṛti','Śūla','Gaṇḍa','Vṛddhi','Dhruva','Vyāghāta','Harṣaṇa','Vajra','Siddhi','Vyatīpāta','Varīyān','Parigha','Śiva','Siddha','Sādhya','Śubha','Śukla','Brahma','Indra','Vaidhṛti']

# ---------- Vimsottari ----------
ORDER=[('Ke',7),('Ve',20),('Su',6),('Mo',10),('Ma',7),('Ra',18),('Ju',16),('Sa',19),('Me',17)]
def vim(jd0, moonlon, yr=365.25):
    n=moonlon/(360/27); i=int(n); frac=n-i
    k=i%9
    start=jd0-frac*ORDER[k][1]*yr
    out=[]
    t=start
    for a in range(9):
        L,y=ORDER[(k+a)%9]
        subs=[]; ts=t
        for b in range(9):
            L2,y2=ORDER[(k+a+b)%9]
            d=y*y2/120*yr
            subs.append((L2,ts,ts+d)); ts+=d
        out.append((L,t,t+y*yr,subs)); t+=y*yr
    return out

if __name__=='__main__':
    RR=(42.275,-83.125)                 # Great Lakes Engineering Works, River Rouge
    SUP=(46.74,-92.08)                  # Burlington Northern Dock No. 1, Superior
    WRK=(46+59.9/60,-(85+6.6/60))       # the wreck
    show(chart(jd(1958,6,7,12,34,tz=-5),*RR),'launch 12:34 EST, 7 June 1958')
    show(chart(jd(1958,6,7,12,38,tz=-5),*RR),'launch 12:38 EST (the other reported minute)')
    show(chart(jd(1975,11,9,14,15,tz=-5),*SUP),'departure 14:15 EST, 9 Nov 1975')
    show(chart(jd(1975,11,10,19,10,tz=-5),*WRK),'last word 19:10 EST, 10 Nov 1975')

    for label,(a,b) in {'launch':(jd(1958,6,6,0,0,tz=EST),jd(1958,6,9,0,0,tz=EST)),'voyage':(jd(1975,11,8,0,0,tz=EST),jd(1975,11,12,0,0,tz=EST))}.items():
        print('=====',label)
        print(' tithi:',[(cal(t),i+1) for t,i in crossings(elong,12,a,b)])
        print(' karana:',[(cal(t),karana_name(i)) for t,i in crossings(elong,6,a,b)])
        print(' naksh:',[(cal(t),NAK[i]) for t,i in crossings(lambda j:lon(j,swe.MOON),360/27,a,b)])
        print(' yoga:',[(cal(t),YOG[i]) for t,i in crossings(ysum,360/27,a,b)])
    # Mars station 1975
    j=jd(1975,10,25,0,0); 
    t=solve(lambda x: spd(x,swe.MARS), j+30, j)  # speed + to -
    t=solve(lambda x: -spd(x,swe.MARS), jd(1975,10,20,0,0), jd(1975,11,20,0,0))
    print('Mars station R', cal(t), dms(lon(t,swe.MARS)%30), lon(t,swe.MARS))
    t2=solve(lambda x: spd(x,swe.MARS), jd(1975,12,20,0,0), jd(1976,2,20,0,0)); print('Mars direct',cal(t2),lon(t2,swe.MARS))
    # Mars in Ardra interval (66.667-80)
    def find_cross(p,deg,a,b,step=0.5):
        res=[]; j=a; prev=ang(lon(j,p)-deg)
        while j<b:
            cur=ang(lon(j+step,p)-deg)
            if (prev<0)!=(cur<0) and abs(cur)<20: res.append(cal(solve(lambda x: ang(lon(x,p)-deg) if prev<0 else -ang(lon(x,p)-deg), j,j+step)))
            prev=cur; j+=step
        return res
    print('Mars 66.667:',find_cross(swe.MARS,66.6667,jd(1975,8,1,0,0),jd(1976,6,1,0,0)))
    print('Mars 80:',find_cross(swe.MARS,80,jd(1975,8,1,0,0),jd(1976,6,1,0,0)))
    print('Mars 60:',find_cross(swe.MARS,60,jd(1975,8,1,0,0),jd(1976,6,1,0,0)))
    print('Saturn Pusya 93.333:',find_cross(swe.SATURN,93.3333,jd(1975,1,1,0,0),jd(1977,1,1,0,0)))
    print('Saturn 106.667:',find_cross(swe.SATURN,106.6667,jd(1975,1,1,0,0),jd(1977,6,1,0,0)))
    print('Saturn 90:',find_cross(swe.SATURN,90,jd(1974,6,1,0,0),jd(1977,1,1,0,0)))
    ts=solve(lambda x: -spd(x,swe.SATURN), jd(1975,11,1,0,0), jd(1975,12,30,0,0)); print('Saturn station R',cal(ts),lon(ts,swe.SATURN))
    print('Jupiter stations', cal(solve(lambda x: -spd(x,swe.JUPITER), jd(1975,7,20,0,0), jd(1975,9,10,0,0))), cal(solve(lambda x: spd(x,swe.JUPITER), jd(1975,11,20,0,0), jd(1976,1,20,0,0))))
    # eclipses
    r=swe.sol_eclipse_when_glob(jd(1975,10,20,0,0),swe.FLG_MOSEPH,0); print('solar ecl',r[0],cal(r[1][0],0),'UT', 'sunlon',lon(r[1][0],swe.SUN), 'moon', lon(r[1][0],swe.MOON))
    try:
        w=swe.sol_eclipse_where(r[1][0],swe.FLG_MOSEPH); print(' where',w[1][:2])
        h=swe.sol_eclipse_how(r[1][0],(-85.11,46.998,183),swe.FLG_MOSEPH); print(' how at wreck',h[0],h[1][:3])
    except Exception as e: print(e)
    r=swe.lun_eclipse_when(jd(1975,11,5,0,0),swe.FLG_MOSEPH,0); print('lunar ecl',r[0],cal(r[1][0],0),'UT','moon',lon(r[1][0],swe.MOON),'sun',lon(r[1][0],swe.SUN),[cal(x) if x else None for x in r[1][:8]])
    # keel day moon
    for h in (0,7,12,17,24): print('keel moon',h,lon(jd(1957,8,7,h,0,tz=EST),swe.MOON))
    # sunset at wreck Nov 10 1975, sunrise
    geo=(-85.11,46.998,183)
    rs=swe.rise_trans(jd(1975,11,10,5,0,tz=EST),swe.SUN,swe.CALC_SET,geo,flags=swe.FLG_MOSEPH); print('sunset wreck',cal(rs[1][0]))
    rs=swe.rise_trans(jd(1975,11,10,0,0,tz=EST),swe.SUN,swe.CALC_RISE,geo,flags=swe.FLG_MOSEPH); print('sunrise wreck',cal(rs[1][0]))
    rs=swe.rise_trans(jd(1975,11,10,12,0,tz=EST),swe.MOON,swe.CALC_SET,geo,flags=swe.FLG_MOSEPH); print('moonset wreck',cal(rs[1][0]))
    rs=swe.rise_trans(jd(1975,11,9,0,0,tz=EST),swe.SUN,swe.CALC_SET,(-92.08,46.74,183),flags=swe.FLG_MOSEPH); print('sunset superior nov9',cal(rs[1][0]))
    rs=swe.rise_trans(jd(1958,6,7,0,0,tz=EST),swe.SUN,swe.CALC_RISE,(-83.125,42.275,175),flags=swe.FLG_MOSEPH); print('sunrise RR',cal(rs[1][0]))
    # lagna sign windows at wreck evening Nov 10
    def asc(j,lat,lo): return swe.houses_ex(j,lat,lo,b'W',swe.FLG_SIDEREAL)[1][0]
    for (nm,la,lo,a,b) in [('wreck',46.998,-85.11,jd(1975,11,10,16,0,tz=EST),jd(1975,11,10,22,0,tz=EST)),('superior',46.74,-92.08,jd(1975,11,9,11,0,tz=EST),jd(1975,11,9,17,0,tz=EST)),('RR',42.275,-83.125,jd(1958,6,7,9,0,tz=EST),jd(1958,6,7,15,0,tz=EST))]:
        print(nm,[(cal(t),SIGNS[i]) for t,i in crossings(lambda j:asc(j,la,lo),30,a,b,dj=0.002)])
    # navamsa boundaries at launch
    for mm in range(28,44,2):
        a_=asc(jd(1958,6,7,12,mm,tz=EST),42.275,-83.125); print('launch',mm,dms(a_%30), 'nav', SIGNS[int((a_*9)//30)%12], 'drek', int((a_%30)//10)+1)
    for mm in (0,15,30,45): 
        a_=asc(jd(1975,11,9,14,mm,tz=EST),46.74,-92.08); print('dep 14:%02d'%mm,SIGNS[int(a_//30)],dms(a_%30),'nav',SIGNS[int((a_*9)//30)%12])
    a_=asc(jd(1975,11,9,15,15,tz=EST),46.74,-92.08); print('dep 15:15 EST (if CST)',SIGNS[int(a_//30)],dms(a_%30))

    # ----- Vimsottari -----

    for mm in (34,38):
        j=jd(1958,6,7,12,mm,tz=EST); m=lon(j,swe.MOON)
        for yr in (365.25,360):
            print('--- launch 12:%d'%mm,'yr',yr)
            for L,a,b,subs in vim(j,m,yr)[:3]:
                print(L,cal(a)[:10],cal(b)[:10])
                if L=='Ra':
                    for s in subs: print('    Ra-'+s[0],cal(s[1])[:10],cal(s[2])[:10])
    jsink=jd(1975,11,10,19,15,tz=EST)
    j=jd(1958,6,7,12,34,tz=EST); print('age yrs',(jsink-j)/365.25, 'days',jsink-j)
    # pratyantar at sinking
    for L,a,b,subs in vim(j,lon(j,swe.MOON)):
        for s in subs:
            if s[1]<=jsink<s[2]:
                print('at sinking:',L,s[0],cal(s[1]),cal(s[2]))
                # pratyantara
                k=[x[0] for x in ORDER].index(s[0]); t=s[1]; tot=s[2]-s[1]
                for c in range(9):
                    L3,y3=ORDER[(k+c)%9]; d=tot*y3/120
                    if t<=jsink<t+d: print('  pratyantara',L3,cal(t),cal(t+d))
                    t+=d
    # how many minutes of launch time shift moves Ra-Ma start past sinking?
    for mm in range(0,60,5):
        j2=jd(1958,6,7,12,mm,tz=EST)
        for L,a,b,subs in vim(j2,lon(j2,swe.MOON)):
            if L=='Ra': print(mm,[ (s[0],cal(s[1])[:10]) for s in subs if s[0] in('Ma',)])