import csv,sys,collections
import pycountry
PLAT=sys.argv[1] if len(sys.argv)>1 else "plateaux_ceinture_soleil.csv"
REB=sys.argv[2] if len(sys.argv)>2 else "plateaux_rebord_vent.csv"
CRETES=sys.argv[3] if len(sys.argv)>3 else "cretes_vent_v14.csv"
V21=sys.argv[4] if len(sys.argv)>4 else "v21_retenues_1mw.csv"
PRUDENCE_VENT=float(sys.argv[5]) if len(sys.argv)>5 else 0.9
AVEC_V21=(sys.argv[6]!="0") if len(sys.argv)>6 else True
SOLEIL=[(1800,871.3),(1900,917.4),(2000,959.8),(2100,998.5),(2200,1033.7),(2300,1065.7),(2400,1094.8),(2500,1121.5),(2600,1145.6),(2700,1167.1)]
MIXTE=[(1800,940.1),(2100,1002.3),(2400,1054.8),(2700,1096.3)]
VENT=[(5.5,662.8),(5.75,715.1),(6.0,766.0),(6.25,815.2),(6.5,862.7),(6.75,908.1),(7.0,951.0),(7.5,1030.3),(8.0,1099.1),(8.5,1155.5),(9.0,1201.0)]
CALAGE_V14=885.0/899.4
def interp(t,h):
    h=min(max(h,t[0][0]),t[-1][0])
    for (a,x),(b,y) in zip(t,t[1:]):
        if a<=h<=b: return x+(y-x)*(h-a)/(b-a)
    return t[-1][1]
COUT_SOLEIL=465.0;COUT_MIXTE=461.0;COUT_CRETE=364.0;EAU_APPOINT=15.0;PACHA=(155.6,31.1);GOLD=(56.0,6.4);FRAC_GOLD=0.47;PART_V21=2.0/3.0
LIGNES={"V14A":(5000,4425.0,1820.0),"V11":(4600,1047.0,501.0),"V56":(4500,597.0,334.0),"V3":(20000,496.0,254.0),"V14F":(1500,1328.0,546.0),"V14S":(4110,4102.0,1911.0),"V14M":(4110,4118.0,1895.0),"P2N":(312,48.5,9.7),"GN":(147,8.2,0.9),"PAN":(2700000,244.0,144.0),"JOU":(0,14.0,8.0),"V15":(1000,405.0,747.0),"TC":(1500,43.0,49.0)}
REMPLACE=["V3","V56","V11","PAN","JOU"]
reb={(round(float(r["lat"]),3),round(float(r["lon"]),3)):float(r["rebord_km_6ms"]) for r in csv.DictReader(open(REB,encoding="utf-8"))}
def pacha(p): return (PACHA[0]+FRAC_GOLD*GOLD[0],PACHA[1]+FRAC_GOLD*GOLD[1]) if p>=1500 else (0.0,0.0)
cand=[]
for r in csv.DictReader(open(PLAT,encoding="utf-8")):
    if r["classe"]=="C" or r["occupation"]=="occupe" or int(r["dist_grande_ville_km"])>150: continue
    n=int(r["sites"]) if r["occupation"]=="libre" else int(r["sites"])//2
    h=int(r["heures_tracker"]);p=int(r["pluie_mm"]);km=reb.get((round(float(r["lat"]),3),round(float(r["lon"]),3)),0.0)
    nmix=int(km//12) if p>=200 else 0
    gp,kp=pacha(p)
    for i in range(n):
        gs=interp(SOLEIL,h);ks=COUT_SOLEIL+(EAU_APPOINT if p<75 else 0.0)
        typ,g,k="soleil",gs,ks
        if i<nmix:
            gm=interp(MIXTE,h)
            if COUT_MIXTE/gm<ks/gs: typ,g,k="mixte",gm,COUT_MIXTE
        cand.append(dict(type=typ,g=g+gp,k=k+kp,poids=1.0,pacha=p>=1500,eau=(typ=="soleil" and p<75),lat=r["lat"],lon=r["lon"],pays=r["pays"],ville=r["ville"],grande_ville=r["grande_ville"],dist=r["dist_grande_ville_km"],h=h,p=p,occ=r["occupation"],vent="",g0=g,k0=k))
for r in csv.DictReader(open(CRETES,encoding="utf-8")):
    if int(r["empl_AB"])==0 or r["occupation"]=="occupe" or int(r["dist_grande_ville_km"])>150 or r["vent_150m_moyen"]=="": continue
    n=int(r["empl_AB"]) if r["occupation"]=="libre" else int(r["empl_AB"])//2
    v=float(r["vent_150m_moyen"])*PRUDENCE_VENT;p=int(r["pluie_mm"]);g=interp(VENT,v)*CALAGE_V14;gp,kp=pacha(p)
    for i in range(n):
        cand.append(dict(type="crete",g=g+gp,k=COUT_CRETE+kp,poids=1.0,pacha=p>=1500,eau=False,lat=r["lat"],lon=r["lon"],pays=r["pays"],ville=r["ville"],grande_ville=r["grande_ville"],dist=r["dist_grande_ville_km"],h="",p=p,occ=r["occupation"],vent=round(v,2),g0=g,k0=COUT_CRETE))
if AVEC_V21:
    for r in csv.DictReader(open(V21,encoding="utf-8")):
        if r["type"]!="V21 Soleil" or float(r["gwh_v21"])<=0: continue
        g=float(r["gwh_v21"]);k=float(r["meur_v21"])
        pc=pycountry.countries.get(alpha_3=r["pays"])
        cand.append(dict(type="v21",pv=float(r["pv_mwc"])/1000.0,g=g,k=k,poids=PART_V21,pacha=False,eau=False,lat=r["lat"],lon=r["lon"],pays=pc.alpha_2 if pc else r["pays"],ville=r["nom"],grande_ville="",dist="",h=r["heures_tracker"],p="",occ="",vent="",g0=g,k0=k))
cand.sort(key=lambda c:c["k"]/c["g"])
cible=LIGNES["V14S"][1]+LIGNES["V14M"][1]+LIGNES["P2N"][1]+LIGNES["GN"][1]+LIGNES["V14A"][1]+sum(LIGNES[k][1] for k in REMPLACE)
cout_avant=LIGNES["V14S"][2]+LIGNES["V14M"][2]+LIGNES["P2N"][2]+LIGNES["GN"][2]+LIGNES["V14A"][2]+sum(LIGNES[k][2] for k in REMPLACE)
A3=dict(twh=10675.0+4425.0,md=4586.0+1820.0+11.9,sites=9761+5000)
twh=md=0.0;pris=[];T=collections.defaultdict(lambda:[0.0,0.0,0.0,0.0,0.0,0.0]);i=0
while twh<cible and i<len(cand):
    c=cand[i];i+=1;w=c["poids"]
    twh+=c["g"]*w/1000.0;md+=c["k"]*w/1000.0;pris.append(c)
    t=T[c["type"]];t[0]+=w;t[1]+=c["g0"]*w/1000.0;t[2]+=c["k0"]*w/1000.0
    if c["pacha"]: t[3]+=w;t[4]+=(c["g"]-c["g0"])*w/1000.0;t[5]+=(c["k"]-c["k0"])*w/1000.0
dernier=pris[-1]["k"]/pris[-1]["g"]
n_sol=round(T["soleil"][0]);n_mix=round(T["mixte"][0]);n_cr=round(T["crete"][0]);n_v21=sum(1 for c in pris if c["type"]=="v21");n_v21p=round(T["v21"][0])
n_eau=sum(1 for c in pris if c["eau"]);n_pacha=round(sum(T[k][3] for k in T));n_gold=round(FRAC_GOLD*n_pacha)
twh_pg=sum(T[k][4] for k in T);md_pg=sum(T[k][5] for k in T)
n_sky=round(n_gold*225.0/1126);sky_md=n_sky*87.2/1000.0
print(f"Candidats : {len(cand)} (plateaux {sum(1 for c in cand if c['type'] in ('soleil','mixte'))}, dont mixtes possibles {sum(1 for c in cand if c['type']=='mixte')} ; cretes {sum(1 for c in cand if c['type']=='crete')} ; barrages V21 au soleil {sum(1 for c in cand if c['type']=='v21')})")
print(f"Cible : {cible:.0f} TWh, cout avant {cout_avant:.0f} Md€ ; Artesia 3 sur ces lignes : {A3['twh']:.0f} TWh, {A3['md']:.0f} Md€ avec ses Sky")
for k in ("crete","soleil","mixte","v21"):
    t=T[k];print(f"{k} : {t[0]:.0f} centrales, {t[1]:.0f} TWh, {t[2]:.0f} Md€ ({t[2]/max(t[1],1e-9):.3f} M€/GWh) ; Pacha {t[3]:.0f}, {t[4]:.0f} TWh, {t[5]:.1f} Md€")
print(f"V21 : {n_v21} barrages, comptes pour deux tiers ({n_v21p}) ; eau d'appoint {n_eau} ; Pacha 2 {n_pacha}, Gold {n_gold}, {twh_pg:.0f} TWh, {md_pg:.1f} Md€ ; Sky 2 nouvelles {n_sky}, {sky_md:.1f} Md€")
print(f"Flotte : {twh:.0f} TWh pour {md:.0f} Md€ ({md/twh:.3f}) ; dernier pris {dernier:.3f} M€/GWh ; economie sur Artesia 3 pour ces lignes {A3['md']-md-sky_md:.0f} Md€")
hs=[c["h"] for c in pris if c["type"] in ("soleil","mixte")];vs=[c["vent"] for c in pris if c["type"]=="crete"]
print(f"Heures moyennes des plateaux pris {sum(hs)/max(len(hs),1):.0f}, min {min(hs)} ; vent prudent moyen des cretes prises {sum(vs)/max(len(vs),1):.2f} m/s, min {min(vs) if vs else 0}")
garde={k:v for k,v in LIGNES.items() if k in ("V14F","V15","TC")}
pacha_twh=420.0+twh_pg*PACHA[0]/(PACHA[0]+FRAC_GOLD*GOLD[0]);pacha_md=84.0+md_pg*PACHA[1]/(PACHA[1]+FRAC_GOLD*GOLD[1])
gold_twh=63.0+twh_pg-(pacha_twh-420.0);gold_md=7.2+md_pg-(pacha_md-84.0)
sky_tot=45.1+sky_md
ts={k:T[k][1] for k in T};ms={k:T[k][2] for k in T}
total_twh=sum(v[1] for v in garde.values())+sum(ts.values())+pacha_twh+gold_twh
B=sum(v[2] for v in garde.values())+sum(ms.values())+pacha_md+gold_md+sky_tot+1000.0
print(f"Tableau B : V14 de crete {n_cr} {ts['crete']:.0f} TWh {ms['crete']:.0f} Md€ ; V14 fossiles 1500 1328 TWh 546 Md€ ; V14 Soleil {n_sol} {ts['soleil']:.0f} TWh {ms['soleil']:.0f} Md€ ; mixtes {n_mix} {ts['mixte']:.0f} TWh {ms['mixte']:.0f} Md€ ; V21 {n_v21} barrages {ts['v21']:.0f} TWh {ms['v21']:.0f} Md€ ; Pacha {2700+n_pacha} {pacha_twh:.0f} TWh {pacha_md:.0f} Md€ ; Gold {1126+n_gold} {gold_twh:.0f} TWh {gold_md:.1f} Md€ ; Sky {356+n_sky} {sky_tot:.1f} Md€ ; V15 405 TWh 747 Md€ ; Titan 43 TWh 49 Md€ ; reseaux 1000 Md€")
print(f"Total : {total_twh:.0f} TWh, scenario B {B:.0f} Md€ ({B/17360:.3f} M€ par GWh annuel) ; economie sur Artesia 3 {8897-B:.0f} Md€ ; marge sous 10 000 : {10000-B:.0f}")
sites=n_sol+n_mix+n_cr+n_v21+1500+1000+1500+(2700+n_pacha)+(1126+n_gold)+(356+n_sky)
hyd=45220-37320+n_sol+n_mix-5000+n_cr
francis=116.0*4425.0/(4425.0+1047.0)*n_cr/5000.0+(1500+n_cr+n_sol+n_mix)*13.5/1000.0
print(f"Sites Hydroloop {hyd} (Artesia 3 : 17 661) plus {n_v21} barrages existants ; Francis {francis:.0f} Md€ ; scenario A {B-francis:.0f} Md€")
eol=2779-5000*0.24+n_cr*0.24-4880*0.180+n_mix*0.180
pv=4230-5000*0.016+n_cr*0.016-4881*0.5-4880*0.2+n_sol*0.5+n_mix*0.2
pv_v21=sum(c.get("pv",0.0)*c["poids"] for c in pris)
pil=1096-(9761+5000)*0.048+(n_sol+n_mix+n_cr)*0.048
sto=14.0-(9761+5000)*0.00065+(n_sol+n_mix+n_cr)*0.00065
eau=629+(n_pacha-1444)*410.0/2700
reserve=7.6+(n_sky-136)*0.3/29
pluie_turb=17.0*(2700+n_pacha)/4144
volet_eau=179.0+(n_sky-136)*130.6/1000.0
print(f"Eoliennes {eol:.0f} GW ; PV plateaux et cretes {pv:.0f} GWc, barrages V21 {pv_v21:.0f} GWc ; pilotable le soir {pil:.0f} GW (hors V21) ; stockage quotidien {sto:.1f} TWh ; reserve {reserve:.1f} TWh ; pluie turbinable {pluie_turb:.0f} TWh ; eau regulee {eau:.0f} km3 ; volet eau {volet_eau:.0f} Md€ ; Sky {356+n_sky} ; Station Castle 2 {254+n_sky}")
print(f"Par an : B {B/25:.0f} Md€, A {(B-francis)/25:.0f} ; part du budget militaire mondial : {B/25/2500*100:.0f} % ; chantiers par an en phase 2 : environ {round((hyd*0.7+1000+1500)/10/100)*100}")
with open("plan_mondial_artesia_4_sites.csv","w",newline="",encoding="utf-8") as f:
    w=csv.writer(f,lineterminator="\n");w.writerow(["rang","centrale","lat","lon","pays","ville_ou_barrage","grande_ville","dist_grande_ville_km","heures_tracker","vent_150m_prudent","pluie_mm","occupation","pacha_2","eau_appoint","poids","gwh_vendus","cout_M_euros","M_euros_par_GWh"])
    for j,c in enumerate(pris,1): w.writerow([j,c["type"],c["lat"],c["lon"],c["pays"],c["ville"],c["grande_ville"],c["dist"],c["h"],c["vent"],c["p"],c["occ"],"oui" if c["pacha"] else "non","oui" if c["eau"] else "non",round(c["poids"],3),round(c["g"]),round(c["k"],1),round(c["k"]/c["g"],3)])
P=collections.defaultdict(lambda:[0,0,0,0,0.0])
for c in pris:
    t=P[c["pays"]];t[{"crete":0,"soleil":1,"mixte":2,"v21":3}[c["type"]]]+=1;t[4]+=c["g"]*c["poids"]/1000.0
with open("plan_mondial_artesia_4_pays.csv","w",newline="",encoding="utf-8") as f:
    w=csv.writer(f,lineterminator="\n");w.writerow(["pays","v14_crete","v14_soleil","mixtes","barrages_v21","twh"])
    for k,t in sorted(P.items(),key=lambda x:-x[1][4]): w.writerow([k]+t[:4]+[round(t[4])])