"""Monte Carlo: first Starship ship-to-ship propellant transfer by 2027-06-30 (register entry 2026-10-038).
All durations in months from 2026-10-23. Inputs are stated assumptions, not data.

    python3 scripts/forecasts/2026-10-038_mc.py

Needs numpy. Seed 20261023, 1,000,000 trials: the output is the same on every run."""
import numpy as np
rng=np.random.default_rng(20261023)
N=1_000_000
DEADLINE=(np.datetime64('2027-06-30')-np.datetime64('2026-10-23')).astype(int)/30.44  # ~8.25 mo

def lognorm(median, p90, n):
    s=np.log(p90/median)/1.2816
    return median*np.exp(s*rng.standard_normal(n))

# ---------- Model A: bottom-up ----------
# A1 first demo vehicle: "in build and integration" today -> launch-ready
a1=lognorm(4.0, 7.5, N)
# A2 second vehicle ready and launched within the demo window (after the first is ready)
a2=lognorm(1.0, 3.0, N)
# A3 intervening-flight mishaps: flights before the demo ~3 (Flight 15 + others), each mishap p=0.25,
#     each causing a stand-down of 1.5 mo median (p90 4 mo)
flights=3; p_mishap=0.25
k=rng.binomial(flights,p_mishap,N)
standdown=np.zeros(N)
for i in range(1,flights+1):
    m=k>=i
    standdown[m]+=lognorm(1.5,4.0,m.sum())
# A4 first attempt achieves a measurable transfer: p=0.70; else retry after median 2.5 mo
ok1=rng.random(N)<0.70
retry=np.where(ok1,0.0,lognorm(2.5,5.0,N))
ok2=rng.random(N)<0.80   # second attempt success
T_A=a1+a2+standdown+retry
succ_A=ok1|ok2           # beyond two attempts treated as past deadline
pA=np.mean((T_A<=DEADLINE)&succ_A)

# ---------- Model B: slip factor on the stated target ----------
# Stated/reported target "early 2027": uniform Jan 15 - Apr 15 2027
t0=rng.uniform((np.datetime64('2027-01-15')-np.datetime64('2026-10-23')).astype(int)/30.44,
               (np.datetime64('2027-04-15')-np.datetime64('2026-10-23')).astype(int)/30.44,N)
# Remaining-duration slip factor, from this demo's own history: two prior targets each slipped
# ~12 months (NASA OIG IG-26-004). Lognormal median 2.0, p90 3.1, floor 1.0. A generic SpaceX
# milestone slip (median 1.7) gives a blend of ~0.59 instead.
slip=np.maximum(1.0,lognorm(2.0,3.1,N))
T_B=t0*slip
pB=np.mean(T_B<=DEADLINE)

blend=0.5*pA+0.5*pB
def q(x): return np.percentile(x,[10,50,90])
print(f"deadline months from today: {DEADLINE:.2f}")
print(f"Model A bottom-up:   P={pA:.3f}  T p10/p50/p90 = {np.round(q(T_A),1)}")
print(f"Model B slip factor: P={pB:.3f}  T p10/p50/p90 = {np.round(q(T_B),1)}")
print(f"Blend 50/50: {blend:.3f}  -> display {round(blend*20)/20:.2f}")
print(f"correlated_mean={blend:.4f}")   # read by scripts/check_forecast_models.py
# sensitivities
for name,val in [('first-attempt success 0.55',None),('slip median 1.4',None),('slip median 2.2',None)]:
    pass