import math,random,sys
random.seed(2026)
H=8760
VENT_MOYEN=float(sys.argv[1]) if len(sys.argv)>1 else 6.5
CHARGE_MAX=float(sys.argv[2]) if len(sys.argv)>2 else 0.0
CHARGE_MIN=CHARGE_MAX/3
N_EOL=10
P_EOL=6.0
PV_TRACKERS=15.0
PV_FLOTTANT=6.0
PV_DIGUES=1.0
H_TRACKERS=1550
H_FLOTTANT=1170
H_DIGUES=1100
RACCORDEMENT=50.0
P_POMPES=16.0
P_PELTON=12.0
STOCK_MAX=147.0
RENDEMENT_POMPE=0.87
RENDEMENT_TURBINE=0.86
PLUIE_GWH=1.9
PERTES_EOL=0.08
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)*(1-PERTES_EOL) 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)
    lat=45.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.7
for j in range(366):
    c=0.6*c+0.4*random.uniform(0.15,1.0)
    nuage.append(c)
brut=[pv_forme(h)*nuage[h//24] for h in range(H)]
somme=sum(brut)
def pv_serie(puissance,heures):
    return [b/somme*puissance*heures for b in brut]
pv_t=pv_serie(PV_TRACKERS,H_TRACKERS)
pv_f=pv_serie(PV_FLOTTANT,H_FLOTTANT)
pv_d=pv_serie(PV_DIGUES,H_DIGUES)
pv=[pv_t[h]+pv_f[h]+pv_d[h] for h in range(H)]
pluie_h=PLUIE_GWH*1000/H
stock=STOCK_MAX*0.5
vendu=0.0;turbine=0.0;pompe_e=0.0;ecrete=0.0;prod_pv=0.0;prod_eol=0.0;charge_surplus=0.0;charge_base=0.0
heures_pelton=0
for h in range(H):
    heure=h%24
    prod=eol[h]+pv[h]
    prod_pv+=pv[h];prod_eol+=eol[h]
    stock=min(STOCK_MAX,stock+pluie_h)
    export=min(prod,RACCORDEMENT)
    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
    nuit=(heure>=23 or heure<6)
    creux=(heure>=5 and heure<17)
    if (nuit or creux) and stock<STOCK_MAX*0.95:
        marge=P_POMPES-pompage
        dispo=min(export,marge,(STOCK_MAX-stock)/RENDEMENT_POMPE)
        if dispo>0:
            pompage+=dispo
            export-=dispo
    stock+=pompage*RENDEMENT_POMPE
    pompe_e+=pompage
    absorbe=min(surplus,CHARGE_MAX)
    charge_surplus+=absorbe
    surplus-=absorbe
    base=max(0.0,CHARGE_MIN-absorbe)
    base=min(base,export)
    export-=base
    charge_base+=base
    ecrete+=surplus
    turb=0.0
    pointe=(heure>=17 or heure<5)
    if pointe and stock>0:
        turb=min(P_PELTON,stock*RENDEMENT_TURBINE,RACCORDEMENT-export)
        stock-=turb/RENDEMENT_TURBINE
        if turb>0: heures_pelton+=1
    turbine+=turb
    vendu+=export+turb
print(f"Vent moyen a 120 m : {VENT_MOYEN:.1f} m/s, heures equivalentes eolien : {prod_eol/(N_EOL*P_EOL):.0f} h")
print(f"Eolien net : {prod_eol/1000:.1f} GWh, PV : {prod_pv/1000:.1f} GWh, pluie : {PLUIE_GWH:.1f} GWh")
print(f"Energie pompee : {pompe_e/1000:.1f} GWh, turbinee : {turbine/1000:.1f} GWh, ecretee : {ecrete/1000:.1f} GWh")
if CHARGE_MAX>0: print(f"Charge pilotable {CHARGE_MAX:.0f} MW : {charge_surplus/1000:.1f} GWh pris sur le surplus, {charge_base/1000:.1f} GWh pris sur la vente")
print(f"Electricite vendue : {vendu/1000:.1f} GWh, dont pilotable : {100*turbine/vendu:.1f} %, Pelton en marche : {heures_pelton} h")