###### ###### ###### ###### ###### ###### ###### ###### ###### ###### ######
# This script allows to reproduce fig. 3 from the paper
# "A robust active power control algorithm to maximize wind farm power tracking margins in waked conditions"
# by Tamaro, Campagnolo and Bottasso, 2025
###### ###### ###### ###### ###### ###### ###### ###### ###### ###### ######

import numpy as np
import matplotlib.pyplot as plt
plt.close("all")
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "serif",
    "font.serif": ["Helvetica"],
})
plt.rcParams.update({'font.size': 18})
plt.ion()
plt.close("all")
import colormaps as cmaps
from scipy.interpolate import RegularGridInterpolator

def get_pitch(d,u,gamma):
    zero = np.load('data/luts_iea_modelTamaro_v2.npz')
    
    gamma_arr = zero['arr_0']
    u_array   = zero['arr_1']
    curtailer = zero['arr_2']
    pitch_angle  = zero['arr_6']
    
    if u < u_array[0]:
        u = u_array[0]
    elif u > u_array[-1]:
        u = u_array[-1]
    if d < curtailer[0]:
        d = curtailer[0]
    elif d > curtailer[-1]:
        d = curtailer[-1]
    if gamma < gamma_arr[0]:
        gamma = gamma_arr[0]
    elif gamma > gamma_arr[-1]:
        gamma = gamma_arr[-1]
    
    interp = RegularGridInterpolator((u_array, curtailer, gamma_arr), pitch_angle,
                                 bounds_error=False, fill_value=None,method="linear")
    return interp((u,d,gamma))

def get_tsr(d,u,gamma):
    zero = np.load('data/luts_iea_modelTamaro_v2.npz')
    
    gamma_arr = zero['arr_0']
    u_array   = zero['arr_1']
    curtailer = zero['arr_2']
    tippo  = zero['arr_7']
    
    if u < u_array[0]:
        u = u_array[0]
    elif u > u_array[-1]:
        u = u_array[-1]
    if d < curtailer[0]:
        d = curtailer[0]
    elif d > curtailer[-1]:
        d = curtailer[-1]
    if gamma < gamma_arr[0]:
        gamma = gamma_arr[0]
    elif gamma > gamma_arr[-1]:
        gamma = gamma_arr[-1]
        
    interp = RegularGridInterpolator((u_array, curtailer, gamma_arr), tippo,
                                 bounds_error=False, fill_value=None,method="linear")
    return interp((u,d,gamma))

def get_cp_mis(d,u,gamma):
    zero = np.load('data/luts_iea_modelTamaro_v2.npz')
    
    gamma_arr = zero['arr_0']
    u_array   = zero['arr_1']
    curtailer = zero['arr_2']
    cp_g1  = zero['arr_4']

    if u < u_array[0]:
        u = u_array[0]
    elif u > u_array[-1]:
        u = u_array[-1]
    if d < curtailer[0]:
        d = curtailer[0]
    elif d > curtailer[-1]:
        d = curtailer[-1]
    if gamma < gamma_arr[0]:
        gamma = gamma_arr[0]
    elif gamma > gamma_arr[-1]:
        gamma = gamma_arr[-1]
    
    interp = RegularGridInterpolator((u_array, curtailer, gamma_arr), cp_g1,
                                 bounds_error=False, fill_value=None,method="linear")
    return interp((u,d,gamma))


fig,ax=plt.subplots(figsize=(9,3),nrows=1,ncols=2,sharey=False, layout='constrained')

looper = np.arange(0.3,1.01,0.1)
colori = cmaps.gothic_r(looper)
#
GAMMA = np.linspace(-30,30,31)
u = 8

rho = 1.17
P_A = 0.5*rho*np.pi*65**2*0.45*u**3

P_A_diviso_P_R = P_A/3.37e6
# gamma = 10
conto = 0
for d in looper:
    # get_pitch(d,U,gamma)    
    passo = np.zeros(len(GAMMA))
    tipspi = np.zeros(len(GAMMA))
    c=0
    for gamma in GAMMA:
        passo[c] = get_pitch(d,u,gamma)
        tipspi[c] = get_tsr(d,u,gamma)
        c+=1
        
    ax[0].plot(GAMMA,tipspi,color=colori[conto])
    ax[1].plot(GAMMA,passo,color=colori[conto])
    conto += 1

ax[0].set_xlabel('$\gamma$ [$^\circ$]')
ax[1].set_xlabel('$\gamma$ [$^\circ$]')

ax[0].set_ylabel('$\lambda$ [-]')
ax[1].set_ylabel(r'$\theta$ [$^\circ$]')

ax[0].set_xlim([-30,30])
ax[1].set_xlim([-30,30])

ax[0].set_ylim([5,9])
ax[1].set_ylim([0,16])

ax[0].set_yticks([5,6,7,8,9])
ax[1].set_yticks([0,4,8,12,16])

ax[0].set_xticks([-30,-15,0,15,30])
ax[1].set_xticks([-30,-15,0,15,30])

ax[0].text(22.5,8.5,'(a)')
ax[1].text(22.5,14,'(b)')

import matplotlib as mpl
norm = mpl.colors.Normalize(vmin=0, vmax=np.max(looper))
cbar=fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmaps.gothic_r),
             ax=ax[1], orientation='vertical', label='$\epsilon$ [-]',location='right',ticks=[0,0.25,0.5,0.75,1])

cbar.ax.set_yticklabels(np.round(np.array([0,0.25,0.5,0.75,1])*P_A_diviso_P_R,1))  # vertically oriented colorbar