# -*- coding: utf-8 -*-
import numpy as np, pickle, matplotlib
matplotlib.use('Agg'); import matplotlib.pyplot as plt
from datetime import date, timedelta
from srad_fix import srad_mask
plt.rcParams.update({'figure.facecolor':'white','axes.facecolor':'white','savefig.facecolor':'white',
 'font.family':'DejaVu Sans','font.size':9,'axes.labelsize':9,'axes.titlesize':9,'legend.fontsize':8,
 'xtick.labelsize':8,'ytick.labelsize':8,'axes.grid':True,'grid.alpha':0.3,'grid.linewidth':0.5,
 'axes.spines.top':False,'axes.spines.right':False,'lines.linewidth':1.3})
THR=20.0; PMIN=1.0; C=8.6
days=[date(2025,1,1)+timedelta(k) for k in range(365)]
prices=pickle.load(open('prices2025.pkl','rb')); PM=np.array([prices[d] for d in days])
A=srad_mask(days)
DA=np.load('disp_arb.npz')['disp']; DC=np.load('disp_const.npz')['disp']
DS=np.load('fix_disp_coupled.npz')['disp']; DB=np.load('disp_band_real.npz')['disp']
def save(f,n):
    f.savefig(n+'.png',dpi=300,bbox_inches='tight'); f.savefig(n+'.pdf',bbox_inches='tight'); plt.close(f); print(' ->',n)

# F1 despacho: ahora el mensaje es el regimen dual dentro/fuera de ventana
fig,axes=plt.subplots(1,2,figsize=(7.2,3.0),sharey=True)
for ax,dd in zip(axes,[date(2025,1,15),date(2025,5,15)]):
    i=days.index(dd); h=np.arange(25); ext=lambda v:np.append(v,v[-1])
    for x0 in range(24):
        if A[i,x0]: ax.axvspan(x0,x0+1,color='#c0392b',alpha=0.07,lw=0)
    ax.step(h,ext(DC[i]),where='post',color='0.45',ls='-',label='const')
    ax.step(h,ext(DA[i]),where='post',color='#1f77b4',ls='--',label='arb')
    ax.step(h,ext(DS[i]),where='post',color='#c0392b',ls='-',lw=1.6,label='band+srad')
    ax.axhline(PMIN+C,color='#c0392b',lw=0.7,alpha=0.7)
    ax.text(0.3,PMIN+C+0.25,r'$P^{min}+C^{srad}$',color='#c0392b',fontsize=7)
    ax2=ax.twinx(); ax2.step(h,ext(PM[i]),where='post',color='0.25',lw=0.8,ls='-.')
    ax2.set_ylim(-25,140); ax2.grid(False); ax2.spines['top'].set_visible(False)
    if ax is axes[1]: ax2.set_ylabel('Day-ahead price (EUR/MWh)',fontsize=8)
    else: ax2.set_yticklabels([])
    ax.set_xlim(0,24); ax.set_xticks([0,6,12,18,24]); ax.set_xlabel('Hour of day')
    ax.set_title(dd.strftime('%d %b %Y'),fontsize=8.5)
axes[0].set_ylabel('Electrolyser power (MW)'); axes[0].set_ylim(-0.5,11)
axes[0].legend(loc='lower left',ncol=3,framealpha=0.9,handlelength=2.2,columnspacing=0.9,borderpad=0.3)
save(fig,'fig1_dispatch')

# F2 nuevo: LCOH frente a capacidad comprometida, con el secuencial en C=0
sw={}
import glob
for f in glob.glob('srad_fix_*.pkl'): sw.update(pickle.load(open(f,'rb')))
FIX=2006184.0; kg=2800*365; H=int(A.sum())
Cs=sorted(sw); L=[(FIX+sw[c]['ec']-sw[c]['br']-42.62*c*H)/kg for c in Cs]
fig,ax=plt.subplots(figsize=(5.0,3.2))
ax.plot(Cs,L,'-o',color='#c0392b',ms=4,label='joint co-optimisation')
ax.plot([0],[L[0]],'s',color='#1f77b4',ms=7,zorder=5)
ax.annotate('sequential:\nzero committable capacity',xy=(0.05,L[0]),xytext=(1.6,4.42),
            fontsize=7.5,color='#1f77b4',arrowprops=dict(arrowstyle='->',color='#1f77b4',lw=0.8))
imin=int(np.argmin(L))
ax.annotate(f'optimum {Cs[imin]:.1f} MW\n{L[imin]:.2f} EUR/kg',xy=(Cs[imin],L[imin]),xytext=(5.0,3.72),
            fontsize=7.5,color='#c0392b',arrowprops=dict(arrowstyle='->',color='#c0392b',lw=0.8))
ax.set_xlabel('Period-wide committed SRAD capacity (MW)'); ax.set_ylabel('LCOH (EUR/kg)')
ax.set_xlim(-0.4,9.2)
save(fig,'fig2_joint_vs_sequential')

# F4 estacionalidad
months=np.array([d.month for d in days])
fig,ax=plt.subplots(figsize=(5.4,3.0))
for lab,D,c,ls,mk in [('const',DC,'0.45','-','o'),('arb',DA,'#1f77b4','--','s'),
                      ('band',DB,'#2ca02c','-.','^'),('band+srad',DS,'#c0392b',':','D')]:
    ov=[]
    for m in range(1,13):
        sel=months==m; Pm=PM[sel]; Dm=D[sel]; pr=Dm>1e-6; t=Pm<=THR
        ov.append(100*np.sum(Dm[pr&t])/max(np.sum(Dm[pr]),1e-9))
    ax.plot(range(1,13),ov,color=c,ls=ls,marker=mk,ms=3.2,label=lab)
hb=[100*np.mean(PM[months==m]<=THR) for m in range(1,13)]
ax.plot(range(1,13),hb,color='0.15',ls=(0,(1,1)),lw=1.0,label='hours below threshold')
ax.set_xlabel('Month of 2025'); ax.set_ylabel('RFNBO overlap, energy-weighted (%)')
ax.set_xticks(range(1,13)); ax.set_xlim(0.7,12.3); ax.set_ylim(0,95)
ax.legend(loc='upper right',ncol=2,framealpha=0.9)
save(fig,'fig4_seasonality')

# F5 sensibilidad precio
ps=np.array([10,20,30,42.62,50,65,80]); best=[]
for p in ps:
    v=[(FIX+sw[c]['ec']-sw[c]['br']-p*c*H)/kg for c in Cs]
    j=int(np.argmin(v)); best.append((Cs[j],v[j]))
fig,ax=plt.subplots(figsize=(5.0,3.2))
ax.plot(ps,[b[1] for b in best],'-o',color='#c0392b',ms=4,label='optimal commitment')
ax.axhline(4.5927,color='#1f77b4',ls=':',lw=1.0); ax.text(11,4.62,'arb',color='#1f77b4',fontsize=7.5)
ax.axhline(4.5307,color='#2ca02c',ls='-.',lw=1.0); ax.text(11,4.36,'band',color='#2ca02c',fontsize=7.5)
for xv,lab in [(42.62,'H2-2026'),(65,'H1-2026')]:
    ax.axvline(xv,color='0.55',lw=0.8,ls=(0,(2,2))); ax.text(xv+1,2.35,lab,rotation=90,fontsize=7,color='0.35')
ax2=ax.twinx(); ax2.plot(ps,[b[0] for b in best],'--s',color='0.45',ms=3,lw=1.0)
ax2.set_ylabel('Optimal $C^{srad}$ (MW)',fontsize=8,color='0.45'); ax2.set_ylim(-0.5,10); ax2.grid(False)
ax.set_xlabel('SRAD availability price (EUR/MW.h)'); ax.set_ylabel('LCOH (EUR/kg)')
ax.set_xlim(8,84); ax.set_ylim(2.1,4.85)
save(fig,'fig5_srad_sensitivity')
print('OK')