"""Checks every number and identity quoted in apse-line-rotation.qmd (apse-line rotation). Prints FAIL on a mismatch."""
import numpy as np
mu=398600.;a=12000.;e=.4;p=a*(1-e*e);k=np.sqrt(mu/p);ok=lambda n,x,y,tol:print(('ok  ' if abs(x-y)<=tol else 'FAIL'),n,round(x,5),y)
ok('p',p,10080,1e-9);ok('k',k,6.288,5e-4);ok('dv30',2*k*e*np.sin(np.radians(15)),1.302,5e-4);ok('rounded dv30',2*6.288*.4*np.sin(np.radians(15)),1.302,5e-4)
for fd,R,vr,vt,v,g in((15,7271,.651,8.718,8.742,4.27),(195,16427,-.651,3.859,3.913,-9.58)):
    f=np.radians(fd);VR=k*e*np.sin(f);VT=k*(1+e*np.cos(f))
    ok('r',p/(1+e*np.cos(f)),R,.5);ok('vr',VR,vr,5e-4);ok('vt',VT,vt,5e-4);ok('v',np.hypot(VR,VT),v,5e-4);ok('g',np.degrees(np.arctan2(VR,VT)),g,5e-3);ok('2|vr|',2*abs(VR),1.302,5e-4)
ra=a*1.4;rp=a*.6;vc=lambda r:np.sqrt(mu/r);vv=lambda r:np.sqrt(mu*(2/r-1/a))
ok('ra',ra,16800,0);ok('vc',vc(ra),4.871,5e-4);ok('va',vv(ra),3.773,5e-4);ok('each',vc(ra)-vv(ra),1.098,5e-4);T=2*(vc(ra)-vv(ra));ok('two',T,2.196,5e-4)
ok('closed form',2*np.sqrt(mu/a)*(1-np.sqrt(1-e))/np.sqrt(1+e),T,1e-12)
tc=np.pi/2*np.sqrt(ra**3/mu);ok('coast s',tc,5418,.5);ok('coast h',tc/3600,1.50,5e-3)
S90=2*k*e*np.sin(np.pi/4);ok('single90',S90,3.557,5e-4);ok('rounded',2*6.288*.4*np.sin(np.pi/4),3.557,5e-4);ok('saving',S90-T,1.361,5e-4)
ok('rp',rp,7200,0);ok('vp',vv(rp),8.804,5e-4);ok('vcp',vc(rp),7.441,5e-4);P=2*(vv(rp)-vc(rp));ok('each',P/2,1.363,5e-4);ok('peri',P,2.726,5e-4)
tp=np.pi/2*np.sqrt(rp**3/mu);ok('coast s',tp,1520,.5);ok('coast h',tp/3600,.42,5e-3);ok('more',P-T,.53,5e-3);ok('ratio',P/T,1.24,5e-3);ok('tratio',tc/tp,3.56,5e-3)
s,t=np.sqrt(1+e),np.sqrt(1-e);ok('ratio formula',s/t*(s-1)/(1-t),P/T,1e-12);ok('p=a s2 t2',a*s*s*t*t,p,1e-9)
sc=t*(1-t)/e;ok('t',t,.7746,5e-5);ok('sin',sc,.4365,5e-5);ok('eta_c',2*np.degrees(np.arcsin(sc)),51.8,.05);ok('exact crossing',2*k*e*sc,T,1e-9);ok('dv180',2*k*e,5.03,5e-3)
E=1e-3;S_=np.sqrt(1+E);T_=np.sqrt(1-E);ok('small e ratio',S_/T_*(S_-1)/(1-T_),1+E/2,1e-5)
# the crossings really are on both orbits with velocities as claimed, for many eta
good=True
for eta in np.radians(np.arange(5,355,5)):
    for f in (eta/2,np.pi+eta/2):
        g=f-eta;good&=abs(np.cos(f)-np.cos(g))<1e-12 and abs(np.sin(f)+np.sin(g))<1e-12 and abs(abs(np.sin(f))-np.sin(eta/2))<1e-12
print('ok  ' if good else 'FAIL','crossing identities for eta in 5..350 deg')
for E,R in ((.01,1.005),(.4,1.24),(.99,6.44)):
    S_=np.sqrt(1+E);T_=np.sqrt(1-E);ok(f'ratio e={E}',S_/T_*(S_-1)/(1-T_),R,5e-3 if R>1.1 else 5e-4)
ok('simplified sin(eta_c/2)',t/(1+t),sc,1e-12);ok('t/(1+t) value',t/(1+t),.4365,5e-5)
g=all(np.sqrt(1-x)/(1+np.sqrt(1-x))<.5 for x in np.linspace(1e-6,.999,999))
print('ok  ' if g else 'FAIL','sin(eta_c/2) < 1/2 for all e, so eta_c < 60 deg')
ok('eta_c as e->0',2*np.degrees(np.arcsin(np.sqrt(1-1e-8)/(1+np.sqrt(1-1e-8)))),60,1e-3)
for E,D in ((.2,56.3),(.8,36.0)):
    T_=np.sqrt(1-E);ok(f'eta_c e={E}',2*np.degrees(np.arcsin(T_/(1+T_))),D,.05)
# the trade-off table (eta = 90 deg, I_sp = 320 s)
ve=9.80665e-3*320
for n,rp_,ra_,E_,A_,P_,ta_,tp_,ma_,mp_,tol in (('LEO',6930.,7070.,.01,.0753,.0757,24.65,23.92,.976,.976,5e-5),('example',7200.,16800.,.40,2.196,2.726,90.3,25.3,.497,.419,5e-4),('GTO',6678.,42164.,.73,2.934,4.852,359,22.6,.393,.213,5e-4)):
    a_=(rp_+ra_)/2;vis_=lambda r:np.sqrt(mu*(2/r-1/a_));vc_=lambda r:np.sqrt(mu/r)
    A=2*(vc_(ra_)-vis_(ra_));P=2*(vis_(rp_)-vc_(rp_));ta=np.pi/2*np.sqrt(ra_**3/mu)/60;tp=np.pi/2*np.sqrt(rp_**3/mu)/60
    ok(n+' e',(ra_-rp_)/(ra_+rp_),E_,5e-3);ok(n+' dv apo',A,A_,tol);ok(n+' dv peri',P,P_,tol);ok(n+' coast apo',ta,ta_,.5 if ta_>100 else .05);ok(n+' coast peri',tp,tp_,.05)
    ok(n+' mass apo',np.exp(-A/ve),ma_,5e-4);ok(n+' mass peri',np.exp(-P/ve),mp_,5e-4)
    if n=='LEO': ok('LEO extra m/s',(P-A)*1e3,.4,.05);ok('LEO saved s',(ta-tp)*60,44,.5);ok('LEO half period min',np.pi*np.sqrt(a_**3/mu)/60,49,.5)
    if n=='GTO': ok('GTO saved h',(ta-tp)/60,5.6,.05);ok('GTO extra',P-A,1.92,5e-3)
print('ok  ' if (.213/.393)<.56 and (.213/.393)>.5 else 'FAIL','GTO gives up almost half of what remained')
ok('30 deg burn uses about a third of the mass',1-np.exp(-1.302/ve),.34,5e-3)
ok('critical inclination',np.degrees(np.arccos(1/np.sqrt(5))),63.4,.05)