import math,random,sys
random.seed(2026)
H=8760
def arg(i,d):
    return float(sys.argv[i]) if len(sys.argv)>i else d
CHUTE=arg(1,1500.0)
VOLUME=arg(2,800000.0)
P_PELTON=arg(3,300.0)
P_POMPES=arg(4,300.0)
PV=arg(5,100.0)
SURPLUS_MW=arg(6,400.0)
PRIX_SURPLUS=arg(7,20.0)
PRIX_POINTE=arg(8,120.0)
P_EOL=arg(9,0.0)
VENT_MOYEN=arg(10,6.5)
N_HYB=arg(11,0.0)
RENDEMENT_POMPE=0.87
RENDEMENT_TURBINE=0.90
STOCK_MAX=VOLUME*CHUTE*9.81*1000*RENDEMENT_TURBINE/3.6e9
H_PV=1700
def pv_forme(h):
    jour=h//24
    heure=h%24
    decl=23.45*math.sin(2*math.pi*(284+jour)/365)
    lat=35.0
    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.4,1.0)
    nuage.append(c)
brut=[pv_forme(h)*nuage[h//24] for h in range(H)]
somme=sum(brut)
pv=[b/somme*PV*H_PV for b in brut]
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
for h in range(H):
    z=0.96*z+math.sqrt(1-0.96*0.96)*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_eol(v):
    if v<3 or v>=25: return 0.0
    if v>=11.5: return 1.0
    for (v1,f1),(v2,f2) in zip(COURBE,COURBE[1:]):
        if v1<=v<v2: return f1+(f2-f1)*(v-v1)/(v2-v1)
    return 1.0
eol=[P_EOL*p_eol(v*1.09)*0.92 for v in vent]
hyb=[N_HYB*(0.062*b/somme*1400+0.010*p_eol(vent[h]*0.85)*0.75) for h,b in enumerate(brut)]
pv=[pv[h]+eol[h]+hyb[h] for h in range(H)]
prod_eol=sum(eol);prod_hyb=sum(hyb)
dispo=[]
for h in range(H):
    heure=h%24
    jour=h//24
    vent=0.5+0.5*random.random()
    soleil=nuage[jour]
    if 10<=heure<16: d=SURPLUS_MW*soleil
    elif heure>=23 or heure<6: d=SURPLUS_MW*vent
    else: d=0.0
    dispo.append(d)
stock=STOCK_MAX*0.5
pompe_pv=0.0;pompe_reseau=0.0;turbine=0.0;pv_vendu=0.0;heures_pelton=0;heures_pompe=0;cycles=0.0
for h in range(H):
    heure=h%24
    pointe=(17<=heure<23) or (6<=heure<9)
    pompage=0.0
    if not pointe and stock<STOCK_MAX:
        marge=P_POMPES
        p1=min(pv[h],marge,(STOCK_MAX-stock)/RENDEMENT_POMPE)
        pompe_pv+=p1
        marge-=p1
        p2=min(dispo[h],marge,(STOCK_MAX-stock)/RENDEMENT_POMPE-p1)
        pompe_reseau+=p2
        pompage=p1+p2
        pv_vendu+=pv[h]-p1
        if pompage>0: heures_pompe+=1
    else:
        pv_vendu+=pv[h]
    stock+=pompage*RENDEMENT_POMPE
    turb=0.0
    if pointe and stock>0:
        turb=min(P_PELTON,stock*RENDEMENT_TURBINE)
        stock-=turb/RENDEMENT_TURBINE
        if turb>0: heures_pelton+=1
    turbine+=turb
cycles=turbine/STOCK_MAX
achat=pompe_reseau*PRIX_SURPLUS/1e6
vente=(turbine*PRIX_POINTE+pv_vendu*70)/1e6
print(f"Chute {CHUTE:.0f} m, bassins 2 x {VOLUME:.0f} m3 : stock utile {STOCK_MAX:.0f} MWh ({STOCK_MAX/P_PELTON:.1f} h de Pelton a {P_PELTON:.0f} MW)")
print(f"Pompes {P_POMPES:.0f} MW, surplus reseau disponible jusqu'a {SURPLUS_MW:.0f} MW la nuit et a la mi-journee, PV {PV:.0f} MWc sur le plateau")
print(f"Production propre : {sum(pv)/1000:.0f} GWh (PV {sum(pv)/1000-prod_eol/1000-prod_hyb/1000:.0f}, eolien {prod_eol/1000:.0f}, panneaux hybrides {prod_hyb/1000:.0f}), dont pompe sur place : {pompe_pv/1000:.0f} GWh, vendu direct : {pv_vendu/1000:.0f} GWh")
print(f"Surplus reseau pompe : {pompe_reseau/1000:.0f} GWh ; energie turbinee en pointe : {turbine/1000:.0f} GWh ; cycles complets : {cycles:.0f} par an")
print(f"Pelton en marche : {heures_pelton} h, pompes en marche : {heures_pompe} h, rendement de boucle : {100*RENDEMENT_POMPE*RENDEMENT_TURBINE:.0f} %")
print(f"Achat du surplus a {PRIX_SURPLUS:.0f} euros/MWh : {achat:.1f} M euros ; vente en pointe a {PRIX_POINTE:.0f} et production directe a 70 : {vente:.1f} M euros par an")