import math,random,sys
H=8760
LAT=25.0
P_EOL=6.0
PRIX=[60]*6+[70]*3+[30]*7+[80]*2+[130]*4+[80]*2
EXPO={"bassins":0.8,"conduites":0.75,"machines":0.9,"acces":0.4,"raccord":0.75,"ligne":0.6,"instrum":0.3,"liaison":0.6,"surface":2/3}
def profils(vent_moyen,graine=2026):
    random.seed(graine)
    def wq(u,m,k=2.0):
        c=m/math.gamma(1+1/k)
        return c*(-math.log(1-u))**(1/k)
    def ncdf(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**2)*random.gauss(0,1)
        s=1+0.15*math.cos(2*math.pi*(h/24-15)/365)
        u=min(max(ncdf(z),1e-6),1-1e-6)
        vent.append(wq(u,vent_moyen*s))
    nuage=[]
    c=0.8
    for j in range(366):
        c=0.6*c+0.4*random.uniform(0.55,1.0)
        nuage.append(c)
    def forme(h):
        j=h//24
        d=23.45*math.sin(2*math.pi*(284+j)/365)
        ha=15*(h%24+0.5-12)
        sa=math.sin(math.radians(LAT))*math.sin(math.radians(d))+math.cos(math.radians(LAT))*math.cos(math.radians(d))*math.cos(math.radians(ha))
        return max(sa,0.0)
    brut=[forme(h)*nuage[h//24] for h in range(H)]
    return vent,brut
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 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
def centrale(k,n_eol=0,vent_moyen=6.2,pv=500.0,h_tr=2100.0,prix_pv=0.57,liaison=0.0):
    vent,brut=profils(vent_moyen)
    surf=k**EXPO["surface"]
    n=n_eol*k
    pvt=pv*k
    pvf=24.0*surf
    racc=320.0*k
    pom=64.0*k
    tur=48.0*k
    stmax=588.0*k
    somme=sum(brut)
    eol=[n*p_eol(v*1.09)*(1-0.08)*1.018 for v in vent]
    prod_pv=[min(b/somme*pvt*h_tr,0.85*pvt)+b/somme*pvf*h_tr*0.75 for b in brut]
    pluie=1.9*1000/H*surf
    stock=stmax*0.5
    vendu=turb_t=pompe_t=ecr=pv_t=eol_t=rec=0.0
    for h in range(H):
        heure=h%24
        prod=eol[h]+prod_pv[h]
        pv_t+=prod_pv[h];eol_t+=eol[h]
        stock=min(stmax,stock+pluie)
        exp=min(prod,racc)
        sur=prod-exp
        pp=0.0
        if sur>0 and stock<stmax:
            pp=min(sur,pom,(stmax-stock)/0.87)
            sur-=pp
        if (heure>=23 or heure<6 or (5<=heure<17)) and stock<stmax*0.95:
            d=min(exp,pom-pp,(stmax-stock)/0.87)
            if d>0:
                pp+=d
                exp-=d
        stock+=pp*0.87
        pompe_t+=pp
        ecr+=sur
        t=0.0
        if (heure>=17 or heure<5) and stock>0:
            t=min(tur,stock*0.86,racc-exp)
            stock-=t/0.86
        turb_t+=t
        vendu+=exp+t
        rec+=(exp+t)*PRIX[heure]
    p={}
    p["eolien"]=n*P_EOL*229/240
    p["trackers"]=pvt*prix_pv
    p["flottant"]=pvf*(prix_pv+0.18)
    p["gouttieres"]=pvt*8.3/500
    p["couverture"]=2.0*surf
    p["machines"]=67.8*k**EXPO["machines"]
    p["conduites"]=13.3*k**EXPO["conduites"]
    p["bassins"]=19.1*k**EXPO["bassins"]
    p["acces"]=2.0*k**EXPO["acces"]
    p["raccord"]=29.3*k**EXPO["raccord"]
    p["instrum"]=0.5*k**EXPO["instrum"]
    p["liaison"]=liaison*k**EXPO["liaison"]
    p["ligne"]=20.0*k**EXPO["ligne"]
    machines=p["eolien"]+p["trackers"]+p["flottant"]+p["machines"]+p["raccord"]+p["instrum"]
    civil=p["bassins"]+p["conduites"]+p["acces"]+p["gouttieres"]+p["couverture"]+p["liaison"]
    invest=machines+civil
    opex=1.5*k**0.5+2.8*n*P_EOL/240+0.006*(pvt+pvf)+0.15
    lp=(machines*0.0578+civil*0.0442+opex)*1e6/vendu
    lg=(machines*0.0465+civil*0.0361+opex)*1e6/vendu
    return dict(k=k,pvt=pvt,pvf=pvf,n=n,racc=racc,stock=stmax,vendu=vendu/1000,pv=pv_t/1000,eol=eol_t/1000,pompe=pompe_t/1000,turbine=turb_t/1000,ecrete=ecr/1000,invest=invest,ligne=p["ligne"],postes=p,opex=opex,lcoe_prive=lp,lcoe_public=lg,valeur=rec/vendu,retour=invest/(rec/1e6-opex))
if __name__=="__main__":
    cas=[("V14 Soleil",dict()),("V14 Soleil mixte",dict(n_eol=30,pv=200.0)),("V14 Soleil de poste 2 100 h",dict(liaison=10.0)),("V14 Soleil de poste 1 500 h",dict(liaison=10.0,h_tr=1500.0))]
    for nom,a in cas:
        for k,lab in [(1,"base"),(2,"MAX"),(4,"MAX 2"),(6,"MAX 3")]:
            r=centrale(k,**a)
            tot=r["invest"]+r["ligne"]
            print(f"{nom} {lab}: vendu {r['vendu']:.0f} GWh, pilotable {100*r['turbine']/r['vendu']:.1f} %, ecrete {r['ecrete']:.0f} GWh, invest {r['invest']:.1f} M, avec ligne {tot:.1f} M, {tot/r['vendu']:.3f} M/GWh, sans ligne {r['invest']/r['vendu']:.3f}, LCOE {r['lcoe_prive']:.1f}/{r['lcoe_public']:.1f}, valeur {r['valeur']:.1f}, retour {r['retour']:.1f} ans")