import math,random,sys
def arg(i,d):
    return float(sys.argv[i]) if len(sys.argv)>i else d
N_EOL=int(arg(1,0))
VENT_MOYEN=arg(2,6.0)
PV_TRACKERS=arg(3,1500.0)
H_TRACKERS=arg(4,2375.0)
POSTE=arg(5,2080.0)
GWH_BARRAGE=arg(6,3400.0)
CHUTE=arg(7,170.0)
P_POMPES=arg(8,600.0)
P_TURBINES=arg(9,2080.0)
LONGUEUR_CONDUITE=arg(10,2.0)
PRIX_PV=arg(11,0.57)
LIGNE_KM=arg(12,20.0)
LAT=36.0
H=8760
random.seed(2026)
P_EOL=6.0
FACTEUR_ROTOR=1.09
PERTES_EOL=0.08
SILLAGE=1.018
RENDEMENT_POMPE=0.87
RENDEMENT_TURBINE=0.90
PRIX=[60]*6+[70]*3+[30]*7+[80]*2+[130]*4+[80]*2
HEURES_BARRAGE=GWH_BARRAGE*1000/365/POSTE
def weibull_quantile(u,mean,k=2.0):
    c=mean/math.gamma(1+1/k)
    return c*(-math.log(1-u))**(1/k)
def normal_cdf(x):
    return 0.5*(1+math.erf(x/math.sqrt(2)))
vent=[]
z=0.0
rho=0.96
for h in range(H):
    z=rho*z+math.sqrt(1-rho*rho)*random.gauss(0,1)
    saison=1+0.15*math.cos(2*math.pi*(h/24-15)/365)
    u=min(max(normal_cdf(z),1e-6),1-1e-6)
    vent.append(weibull_quantile(u,VENT_MOYEN*saison))
COURBE=[(3,0.0),(4,0.05),(5,0.13),(6,0.24),(7,0.39),(8,0.56),(9,0.74),(10,0.89),(11,0.97),(11.5,1.0)]
def p_eolienne(v):
    if v<3 or v>=25: return 0.0
    if v>=11.5: return P_EOL
    for (v1,f1),(v2,f2) in zip(COURBE,COURBE[1:]):
        if v1<=v<v2: return P_EOL*(f1+(f2-f1)*(v-v1)/(v2-v1))
    return P_EOL
eol=[N_EOL*p_eolienne(v*FACTEUR_ROTOR)*(1-PERTES_EOL)*SILLAGE for v in vent]
def pv_forme(h):
    jour=h//24
    heure=h%24
    decl=23.45*math.sin(2*math.pi*(284+jour)/365)
    ha=15*(heure+0.5-12)
    sinalt=math.sin(math.radians(LAT))*math.sin(math.radians(decl))+math.cos(math.radians(LAT))*math.cos(math.radians(decl))*math.cos(math.radians(ha))
    return max(sinalt,0.0)
nuage=[]
c=0.8
for j in range(366):
    c=0.6*c+0.4*random.uniform(0.55,1.0)
    nuage.append(c)
brut=[pv_forme(h)*nuage[h//24] for h in range(H)]
somme=sum(brut)
pv=[min(b/somme*PV_TRACKERS*H_TRACKERS,0.85*PV_TRACKERS) for b in brut]
def barrage(heure):
    debut=17.0
    return POSTE if debut<=heure<debut+HEURES_BARRAGE else (POSTE*(debut+HEURES_BARRAGE-heure) if heure<debut+HEURES_BARRAGE+1 and heure>=debut+math.floor(HEURES_BARRAGE) else 0.0)
KWH_M3=9.81*CHUTE/3.6e6*1000
stock=0.0
vendu=turbine=pompe_e=ecrete=prod_pv=prod_eol=recette=eau=0.0
heures_turbine=0
STOCK_MAX=P_POMPES*RENDEMENT_POMPE*8
for h in range(H):
    heure=h%24
    prod=eol[h]+pv[h]
    prod_pv+=pv[h];prod_eol+=eol[h]
    libre=max(POSTE-barrage(heure),0.0)
    export=min(prod,libre)
    surplus=prod-export
    pompage=0.0
    if surplus>0 and stock<STOCK_MAX:
        pompage=min(surplus,P_POMPES,(STOCK_MAX-stock)/RENDEMENT_POMPE)
        surplus-=pompage
    creux=(heure>=9 and heure<16)
    if creux and stock<STOCK_MAX*0.95:
        dispo=min(export,P_POMPES-pompage,(STOCK_MAX-stock)/RENDEMENT_POMPE)
        if dispo>0:
            pompage+=dispo
            export-=dispo
    stock+=pompage*RENDEMENT_POMPE
    pompe_e+=pompage
    ecrete+=surplus
    turb=0.0
    if (heure>=17 or heure<5) and stock>0:
        turb=min(P_TURBINES,stock*RENDEMENT_TURBINE,libre-export)
        if turb>0:
            stock-=turb/RENDEMENT_TURBINE
            heures_turbine+=1
    turbine+=turb
    vendu+=export+turb
    recette+=(export+turb)*PRIX[heure]
eau=pompe_e*1000/KWH_M3*RENDEMENT_POMPE/1e6
eolien_m=N_EOL*P_EOL*229/240
pv_m=PV_TRACKERS*PRIX_PV
pompes=67.8*P_POMPES/112*0.55
conduite=13.3*P_POMPES/64*LONGUEUR_CONDUITE/1.5
prise=40.0+0.05*P_POMPES
ligne=LIGNE_KM*1.0
poste=15.0+0.01*(PV_TRACKERS+N_EOL*P_EOL)
machines=eolien_m+pv_m+pompes+poste+0.5
civil=conduite+prise+ligne
invest=machines+civil
opex=0.006*PV_TRACKERS+2.8*N_EOL*P_EOL/240+0.004*P_POMPES+1.0
annuite_prive=machines*0.0578+civil*0.0442+opex
annuite_public=machines*0.0465+civil*0.0361+opex
print(f"Barrage existant : poste {POSTE:.0f} MW, {GWH_BARRAGE:.0f} GWh par an par l'eau du fleuve ({HEURES_BARRAGE:.1f} h par jour a pleine puissance), chute {CHUTE:.0f} m, turbines existantes {P_TURBINES:.0f} MW")
print(f"Hydroloop ajoutee : eolien {N_EOL} x {P_EOL:.0f} MW a {VENT_MOYEN:.1f} m/s, PV trackers {PV_TRACKERS:.0f} MWc a {H_TRACKERS:.0f} h, pompes {P_POMPES:.0f} MW, conduite {LONGUEUR_CONDUITE:.1f} km, ligne {LIGNE_KM:.0f} km")
print(f"Eolien net : {prod_eol/1000:.1f} GWh, PV : {prod_pv/1000:.1f} GWh, pompe : {pompe_e/1000:.1f} GWh, turbine : {turbine/1000:.1f} GWh, ecrete : {ecrete/1000:.1f} GWh, eau remontee : {eau:.0f} hm3 par an")
print(f"Electricite vendue en plus : {vendu/1000:.1f} GWh, dont pilotable le soir : {100*turbine/vendu:.1f} %, turbines existantes en marche pour l'Hydroloop : {heures_turbine} h")
print(f"Investissement : {invest:.1f} M euros (eolien {eolien_m:.1f}, PV {pv_m:.1f}, pompes {pompes:.1f}, conduite et prise {conduite+prise:.1f}, ligne et poste {ligne+poste:.1f}), soit {invest/vendu*1000:.3f} M euros par GWh annuel")
print(f"LCOE prive / public : {annuite_prive*1e6/vendu:.1f} / {annuite_public*1e6/vendu:.1f} euros/MWh")
print(f"Valeur au profil de prix ensoleille (moyenne {sum(PRIX)/24:.1f}) : {recette/vendu:.1f} euros/MWh, recette {recette/1e6:.1f} M euros, temps de retour {invest/(recette/1e6-opex):.1f} ans")