import math,random
H=8760
EXPO={"pelton":0.85,"pompes":0.95,"caverne":0.7,"puits":0.75,"bassins":0.8,"retenue":0.8,"raccord":0.75,"instrum":0.3}
BASE={"pelton":210.0,"pompes":210.0,"pv":57.0,"raccord":25.0,"instrum":0.7,"caverne":60.0,"puits":58.0,"bassins":23.0,"retenue":3.0}
MACHINES=("pelton","pompes","pv","raccord","instrum")
def simule(k,chute=1500.0,surplus=400.0):
    random.seed(2026)
    volume=800000.0*k
    pel=300.0*k
    pom=300.0*k
    pvc=100.0*k
    sur=surplus*k
    stmax=volume*chute*9.81*1000*0.90/3.6e9
    def forme(h):
        d=23.45*math.sin(2*math.pi*(284+h//24)/365)
        ha=15*(h%24+0.5-12)
        sa=math.sin(math.radians(35.0))*math.sin(math.radians(d))+math.cos(math.radians(35.0))*math.cos(math.radians(d))*math.cos(math.radians(ha))
        return max(sa,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=[forme(h)*nuage[h//24] for h in range(H)]
    somme=sum(brut)
    pv=[b/somme*pvc*1700 for b in brut]
    for h in range(H):
        random.gauss(0,1)
    dispo=[]
    for h in range(H):
        heure=h%24
        v=0.5+0.5*random.random()
        if 10<=heure<16: d=sur*nuage[h//24]
        elif heure>=23 or heure<6: d=sur*v
        else: d=0.0
        dispo.append(d)
    stock=stmax*0.5
    ppv=pres=turb=pvv=0.0
    for h in range(H):
        heure=h%24
        pointe=(17<=heure<23) or (6<=heure<9)
        pp=0.0
        if not pointe and stock<stmax:
            p1=min(pv[h],pom,(stmax-stock)/0.87)
            p2=min(dispo[h],pom-p1,(stmax-stock)/0.87-p1)
            ppv+=p1;pres+=p2;pp=p1+p2
            pvv+=pv[h]-p1
        else:
            pvv+=pv[h]
        stock+=pp*0.87
        if pointe and stock>0:
            t=min(pel,stock*0.90)
            stock-=t/0.90
            turb+=t
    return dict(stock=stmax,turbine=turb/1000,pompe_reseau=pres/1000,pv_vendu=pvv/1000,pv_pompe=ppv/1000)
def centrale(k,chute=1500.0,surplus=400.0):
    s=simule(k,chute,surplus)
    p={c:BASE[c]*(k if c=="pv" else k**EXPO[c]) for c in BASE}
    mach=sum(p[c] for c in MACHINES)
    civ=sum(p[c] for c in BASE if c not in MACHINES)
    inv=mach+civ
    opex=0.015*inv
    serv_prive=(mach*0.0578+civ*0.0442+opex)*1e6/(s["turbine"]*1000)
    serv_public=(mach*0.0465+civ*0.0361+opex)*1e6/(s["turbine"]*1000)
    rec=s["turbine"]*120/1000+s["pv_vendu"]*70/1000+12*k**1+2*k
    net=rec-s["pompe_reseau"]*20/1000-opex
    return dict(k=k,postes=p,invest=inv,opex=opex,service_prive=serv_prive,service_public=serv_public,livre_prive=serv_prive+s["pompe_reseau"]*20/s["turbine"],net=net,retour=inv/net,**s)
if __name__=="__main__":
    for k,lab in [(1,"V15"),(2,"MAX"),(4,"MAX 2"),(6,"MAX 3")]:
        r=centrale(k)
        print(f"{lab}: stock {r['stock']:.0f} MWh, turbine {r['turbine']:.0f} GWh, invest {r['invest']:.1f} M, {r['invest']/r['stock']*1000:.0f} k/MWh, service {r['service_prive']:.1f}/{r['service_public']:.1f}, livre {r['livre_prive']:.1f}, net {r['net']:.1f}, retour {r['retour']:.1f}")