"""When does the Arctic route pay? Voyage-cost comparison of the Northern Sea Route (NSR), Suez and the Cape.
Run: python3 arctic-route-economics.py   (Python 3.9 or later, standard library only, fixed seed)
Writes arctic-route-economics-results.json beside this file.
Source: https://secopsnomad.si/blog/arctic-ice-melt-geopolitics/

Every input is a range, not a point estimate. Anchors from published sources are cited next to each range;
ranges marked "assumption" are SecOpsNomad's judgement. Costs are per one-way laden voyage in US dollars.
The model leaves out things that favour Suez (extra port calls on the way, liner schedule value), so it
flatters the Arctic route for container ships; read its container results as an upper bound.
"""
import json
import random
import statistics as st
from pathlib import Path

OUT = Path(__file__).resolve().parent / "arctic-route-economics-results.json"
SEED = 20261028
DRAWS = 50_000


def u(rng, lo, hi):
    return rng.uniform(lo, hi)


def tri(rng, lo, mode, hi):
    return rng.triangular(lo, hi, mode)


BUNKER = (650, 900)  # VLSFO, $/t, October 2026: Rotterdam about $683, Singapore $890-904 (Ship & Bunker, ENGINE)


# ---------- container ship, Shanghai -> Rotterdam, about 4,000 TEU ----------
def container(rng, year, waive_arctic_fees=False):
    nm_suez, nm_cape = 10_600, 13_800             # Port Economics, Management and Policy (no port calls)
    nm_nsr = u(rng, 7_600, 8_700)                 # 3,000 nm shorter than Suez (trade press, 2026); 18.1% shorter (Shin & Kim, 2026)
    ice_leg = u(rng, 2_400, 3_000)                # Kara Gate to the Bering Strait, nm (assumption)
    v_open = u(rng, 13, 15)                       # knots (assumption)
    v_ice = u(rng, 8, 12) if year == 2026 else u(rng, 10, 13)          # knots in ice (assumption)
    delay = tri(rng, 0, 1, 6) if year == 2026 else tri(rng, 0, 0.5, 3)  # ice and convoy delays, days (assumption)
    day_cost = u(rng, 35_000, 65_000)             # hire and running cost, $/day (assumption; 2026 feeder fixtures $26,000-62,500)
    ice_premium = u(rng, 1.10, 1.30)              # ice-class ship costs more per day (assumption)
    burn = u(rng, 45, 65)                         # t/day at service speed (assumption)
    ice_burn = burn * u(rng, 1.0, 1.3)            # t/day in ice (assumption)
    bunker = u(rng, *BUNKER)
    suez_toll = u(rng, 250_000, 450_000)          # assumption; Shin & Kim used $300,000 for a Panamax ship
    hull = u(rng, 40e6, 60e6)                     # assumption
    war_risk = hull * (u(rng, 0.002, 0.007) if year == 2026 else u(rng, 0, 0.0005))  # Red Sea additional premium (assumption)
    nsr_fees = 0 if waive_arctic_fees else (u(rng, 150_000, 600_000) if year == 2026 else u(rng, 100_000, 400_000))  # escort, permit, pilots (assumption)
    nsr_insurance = 0 if waive_arctic_fees else u(rng, 50_000, 250_000)  # Arctic hull and P&I loading (assumption)

    sea = lambda nm, v: nm / v / 24
    suez_days = sea(nm_suez, v_open) + u(rng, 0.5, 1.5)
    cape_days = sea(nm_cape, v_open)
    ice_days = sea(ice_leg, v_ice) + delay
    nsr_days = sea(nm_nsr - ice_leg, v_open) + ice_days
    cost = {
        "suez": suez_days * day_cost + sea(nm_suez, v_open) * burn * bunker + suez_toll + war_risk,
        "cape": cape_days * (day_cost + burn * bunker),
        "nsr": nsr_days * day_cost * ice_premium + (nsr_days - ice_days) * burn * bunker + ice_days * ice_burn * bunker
               + nsr_fees + nsr_insurance,
    }
    return cost, {"suez": suez_days, "cape": cape_days, "nsr": nsr_days}


# ---------- Arc7 LNG carrier, Sabetta (Yamal LNG) -> east China ----------
def lng(rng, year, war=True):
    nsr_days = u(rng, 17, 23) + tri(rng, 0, 1, 5)  # Novatek: 19 days via the NSR in July 2018
    suez_days = u(rng, 33, 38)                     # Novatek: 35 days via Suez and Malacca
    cape_days = suez_days + u(rng, 8, 11)          # the Cape adds about 8 days at 16 knots (Port Economics)
    day_cost = u(rng, 70_000, 130_000)             # UBS estimate for an Arc7 charter $70,000-73,000/day, plus running cost (assumption)
    burn = u(rng, 90, 130)                         # t/day fuel equivalent (assumption)
    bunker = u(rng, *BUNKER)
    hull = u(rng, 250e6, 320e6)                    # Arc7 carriers cost about $300 million (assumption)
    war_risk = hull * u(rng, 0.002, 0.007) if (war and year == 2026) else 0
    cost = {
        "suez": suez_days * (day_cost + burn * bunker) + u(rng, 400_000, 800_000) + war_risk,  # toll: assumption
        "cape": cape_days * (day_cost + burn * bunker),
        "nsr": nsr_days * (day_cost + burn * 1.1 * bunker) + u(rng, 50_000, 300_000),  # Arc7 often sails without escort in summer
    }
    return cost, {"suez": suez_days, "cape": cape_days, "nsr": nsr_days}


# ---------- ice-class Aframax crude tanker, Murmansk or an Arctic terminal -> Ningbo / Shandong ----------
def crude(rng, year, war=True):
    nsr_days = u(rng, 19, 24) + tri(rng, 0, 2, 22)  # Rosatomflot: 19.9 days Murmansk-Ningbo; 2023 ice delays pushed one voyage to 46-47 days
    suez_days = u(rng, 35, 40)                     # Rosatomflot: 35.8 days; Ruksha: 40 days to Shanghai
    cape_days = suez_days + u(rng, 8, 12)          # assumption, after Port Economics
    day_cost = u(rng, 40_000, 80_000)              # assumption; one-year VLCC hire near $200,000/day in late September 2026
    burn = u(rng, 35, 50)                          # t/day (assumption)
    bunker = u(rng, *BUNKER)
    hull = u(rng, 50e6, 90e6)                      # assumption
    war_risk = hull * u(rng, 0.002, 0.007) if (war and year == 2026) else 0
    cost = {
        "suez": suez_days * (day_cost + burn * bunker) + u(rng, 250_000, 400_000) + war_risk,  # toll: assumption
        "cape": cape_days * (day_cost + burn * bunker),
        "nsr": nsr_days * (day_cost + burn * 1.1 * bunker) + u(rng, 100_000, 500_000) + u(rng, 30_000, 150_000),
    }
    return cost, {"suez": suez_days, "cape": cape_days, "nsr": nsr_days}


def pct(sorted_values, q):
    return sorted_values[min(len(sorted_values) - 1, max(0, round(q * (len(sorted_values) - 1))))]


def run(rng, fn, *args, **kw):
    """Draw the case DRAWS times; report how often the NSR is cheaper and by how much."""
    save_suez, save_cape, days = [], [], {"suez": [], "cape": [], "nsr": []}
    for _ in range(DRAWS):
        cost, d = fn(rng, *args, **kw)
        save_suez.append(cost["suez"] - cost["nsr"])
        save_cape.append(cost["cape"] - cost["nsr"])
        for k in days:
            days[k].append(d[k])
    out = {}
    for name, s in (("vsSuez", save_suez), ("vsCape", save_cape)):
        s.sort()
        out[name] = {"pNsrCheaper": round(sum(1 for v in s if v > 0) / DRAWS, 3),
                     "medianSaving": round(pct(s, 0.5), -3), "p10": round(pct(s, 0.1), -3), "p90": round(pct(s, 0.9), -3)}
    out["days"] = {k: round(st.median(v), 1) for k, v in days.items()}
    return out


def main():
    rng = random.Random(SEED)
    result = {
        "seed": SEED, "draws": DRAWS, "bunkerUsdPerTonne": list(BUNKER),
        "container": {
            "summer2026": run(rng, container, 2026),
            "summer2026FeesWaived": run(rng, container, 2026, waive_arctic_fees=True),
            "summer2040": run(rng, container, 2040),
        },
        "lng": {"summer2026": run(rng, lng, 2026), "summer2026NoWarRisk": run(rng, lng, 2026, war=False)},
        "crude": {"summer2026": run(rng, crude, 2026), "summer2026NoWarRisk": run(rng, crude, 2026, war=False)},
        # Independent check: Shin & Kim (Future Transportation, 2026) found the NSR cheaper than Suez for a container
        # ship in at most 23.5% of their simulations (easing tensions, summer) and 1.4% in winter.
        "benchmark": {"study": "Shin & Kim 2026, Future Transportation 6(4):168", "maxPNsrCheaperThanSuez": 0.235},
    }
    OUT.write_text(json.dumps(result, ensure_ascii=False, indent=1) + "\n", encoding="utf-8")
    print(json.dumps(result, ensure_ascii=False, indent=1))


if __name__ == "__main__":
    main()