clear
close all
clc

%% IMPORTANT MODEL SETTINGS
Ts = 0.1; %Time step. For better accuracy = 0.01

decimation = 60/Ts; %Inside display and scope blocks. It saves a value per minute.
% The higher the decimation the lower the number of data to be stored
SOC_initial = 100; %Battery initial state of charge

pause_weight = 1;
first_day = '01-Jul';
last_day = '01-Jul';

plot_weather = 0;
plot_propeller = 0;
plot_model = 0;
plot_vel_drag_pow = 0;

model_name = 'Musarra_Tubi_SIM.slx';
save_name = ['RUN_',first_day,'_to_',last_day,'_pw_',num2str(pause_weight)];

%% GRAPHICAL SETTINGS

setDefaults
%colors vectors

blue=[0 0.44706 0.74118];
orange=[0.85098 0.32549 0.098039];
ocra=[0.92941 0.69412 0.12549];
beige = [0.92157 0.74118 0.53333];
grey = [0.50196 0.50196 0.50196];
l_green = [0.74902 0.83137 0.63137];
green = [ 0.46667 0.67451 0.18824];
purple = [0.71765 0.27451 1];
red = [0.63529 0.078431 0.18431];
brown = [0.4902 0.3098 0.3098];
l_brown = [0.67843 0.50196 0.50196];
l_blue = [0.71765 0.86667 0.92941];

%% BATTERY PARAMETERS

%Charging data
load('SOC_R0_V0.mat')
SOC_v = SOC_R0_V0(:,1); %State of charge vector [0,1]
R_int_v = SOC_R0_V0(:,2); %Resistance vector as a function of SOC (from Simulink)
V0_v = SOC_R0_V0(:,3); %Voltage vector as a function of SOC (from Simulink)

% Datasheet parameters
N_ser_batt = 2; %Series
N_par_batt = 4; %Parallel

E_single_battery = 1.84; %[kWh]
E_total_battery = E_single_battery*N_par_batt*N_ser_batt; %[kWh]
V_full = 27.3688*N_ser_batt; %[V]
V_nom = 25.6*N_ser_batt; 
V_exp = 26.4259*N_ser_batt;       
V_cut_off = 20*N_ser_batt;   
V_charge_max = 29.2*N_ser_batt; 

Lambda = 1.038999999968; % From Matlab generic battery model
Q_rated = 71.84*N_par_batt; %[Ah]
Q_max = 1.05*Q_rated; 
Q_exp = 0.0357*Q_rated;
Q_nom = 0.83*Q_rated;   
Q_MAX = Q_max*Lambda;
Q = Q_MAX;
Q_v = SOC_v*Q_rated;

I_nom = 50; %[A]
I_max_charge = 30; %To preserve battery life
I_charge = -Q_rated/10; %10% of rated capacity

R_int_single = mean(R_int_v); %[Ohm]
R_int = R_int_single*N_ser_batt/N_par_batt;

T_response = 0.001; %[s] to obtain a faster simulation

% Extraction of empirical parameters
A = (V_full - V_exp)*1.005;%[V]
V0 = V_full + R_int*I_nom - A;
B_batt  = 3/Q_exp; %[Ah^-1]
K = (-V_exp + V0 -R_int*I_nom + A.*exp(-3)) * (Q-Q_exp)/(Q*(Q_exp+I_nom)); %[Ohm] or [V/Ah]
it_initial = (1-SOC_initial/100)*Q/Lambda*3600; %[Ah]

%% MOTOR PARAMETERS

%Datasheet
P_rated = 12000; %Electrical power [W] 
P_propulsive = 6720; %Propulsive power [W] 
eta_total = 0.56; %From electrical input to propulsive power
Va_rated = 48; %[V] 
eta_motor = 0.96; %Assumed from website
w_rated = 1400*pi/30; %[rad/s]

%Derived
Ia_rated = P_rated/Va_rated; %[A] 
Te_rated = P_rated*eta_motor/w_rated; %[Nm] 
Ra = Va_rated/Ia_rated*(1-eta_motor); %[Ohm] 
La = 1.7e-3; % Armature inductance  [H]  
J = 0.025; % System inertia[kg*m^2] 
B = 0.01; % Viscous friction coefficient [Nm/(rad/s)]

Kt = Te_rated/Ia_rated; %[Nm/A]
Ke = Kt; %[V/(rad/s)]
Ks = 1/Kt; %[A/Nm]

%% BOAT RESISTANCE AND PROPELLER PARAMETERS

%Boat resistance data
load('vel_drag_pow.mat')
L_boat = 9; %[m] overall boat length

c1_v = vel_drag_pow(:,2)./(vel_drag_pow(:,1).^2); %Coefficient

%Data propeller
D = 0.32; %[m] Diameter
Pitch = 0.330; %[m]
Z = 5;  %Number of blades (2 - 7)
w_max_prop = 1400; %[rpm]
np_max = w_max_prop/60; %[rev/s]
PD = Pitch/D;  %Pitch/diameter ratio (0.6 - 1.4)
kp = 1; %Number of propellers

%Assumed
t = 0.15; % [-] Thrust deduction factor
w = 0.2; %[-] Wake factor
AeA0 = 0.75; %Blade area ratio (0.30 - 1.05)
eta_r = 0.99; %relative rotative efficiency
eta_hull = (1-t)/(1-w); %Hull efficiency
v_s_v = vel_drag_pow(:,1); % m/s Boat velocity vector
v_a_v = v_s_v*(1-w); % m/s Velocity of advance vector 

rho = 1025; %[kg/m^3] water density
mu = 1e-3; %[Pa*s] water dynamic viscosity
Re_v = rho.*(v_s_v).*L_boat/mu; %[-] Max Reynolds number vector (2e6 - 2e9)

% N.B. There is basically no difference apported by the Re number at these speeds, so it could not be taken into account

J_min_v = v_a_v./(np_max*D); %Advance ratio based on max propeller speed

R_v = vel_drag_pow(:,2); % [N] Boat resistance vector = f(v_s)

T_v = R_v./(kp*(1-t)); % [N] Required thrust vector

c8_v = c1_v./(kp*(1-t)*(1-w)^2); %Coefficient
c7_v = c8_v./(rho*D^2); %Coefficient

J_v = (0:0.001:1.16)'; %Advance ratio vector

%Calculations 

%Obtainining open water diagram curves thanks to propcalc
[Kt_v, Kq_v, eta0_v] = propcalc(Z, AeA0, PD, J_v, mean(Re_v(2:end))); %Obtain 

Kt_ship_mat = zeros(length(J_v),length(v_s_v));
np_int_v = zeros(length(v_s_v),1);

J_int_v = np_int_v;
Kt_int_v = np_int_v;
Kq_int_v = np_int_v;
eta0_int_v = np_int_v;
T_prop_v = np_int_v;
Q_prop_v = np_int_v;
P_prop_v = np_int_v;
np_rpm_v = np_int_v;

for ii=2:length(v_s_v)

    Kt_ship_mat(:,ii) = c7_v(ii).*J_v.^2;

    int=find((Kt_v-Kt_ship_mat(:,ii))<=0 ,1,'first');

    J_int_v(ii) = J_v(int);
    np_int_v(ii)=v_a_v(ii)/J_int_v(ii)/D;
    
    Kt_int_v(ii) = Kt_v(int);
    Kq_int_v(ii) = Kq_v(int);
    eta0_int_v(ii) = eta0_v(int);
    T_prop_v(ii) = (Kt_v(int)*rho*np_int_v(ii)^2*D^4); %Propeller Thrust [N]
    Q_prop_v(ii) = (Kq_v(int)*rho*(np_int_v(ii))^2*D^5)/eta_r; %Propeller Torque [Nm]
    P_prop_v(ii) = 2*pi*Q_prop_v(ii)*np_int_v(ii); %Propeller load [W]
    np_rpm_v(ii) = np_int_v(ii)*60; %Propeller speed [W]

end

%To be fed to Simulink lookup table
v_s_and_Q_prop = [v_s_v Q_prop_v];
v_s_and_np_rpm = [v_s_v np_rpm_v];
v_s_and_P_prop = [v_s_v P_prop_v];
v_s_and_eta0 = [v_s_v eta0_int_v];
v_s_and_T_prop = [v_s_v T_prop_v];

%% MOTOR CONTROL PARAMETERS
w_ref_3 = 75.088;  %rad/s at 3 m/s  (717.039 rpm)
w_ref_2 = 49.867; %rad/s at 2 m/s   (476.1905 rpm)

Q_ref_3 = 19.855; %Torque load at 3 m/s
Q_ref_2 = 8.6966; %Torque load at 2 m/s

% s = speed  c = current

fs_chopper = 1e3; % 1 [kHz]
w_c = fs_chopper/20; %Frequency inner control loop
w_s = w_c/5; %Frequency outer control loop

Kpc = La*w_c; %Proportional gain coefficient of inner loop
Kic = Ra*w_c; %Integral gain coefficient of inner loop
Kac = 1/Kpc; %Anti wind-up coefficient

Kps = J*5; %Proportional gain coefficient of outer loop
Kis = J*w_s; %Integral gain coefficient of outer loop
Kas = 1/Kps; %Anti wind-up coefficient

%% PV PARAMETERS

Ns_module = 3; %Number of series connected PV modules on a single column
Np_strings = 3; %Number of parallel connected PV modules on a single row

Irr_min = 80; %W/m^2 needed to activate the panels
% 1000 W/m^2 Reference Irradiance
T_0 = 298.15; %K -> 25degC

%Fixed  
Kb = 1.3805e-23; %J/K Boltzmann constant
q = 1.6022e-19; %C electron charge
Eg0 = 1.1; %[eV] Band  gap  energy  for  silicon
dEgdT = -0.0002678; %Temperature coefficient of band gap energy

% Datasheet
Vmp_0 = 30.88; % [V] Module maximum power point voltage at 25degC and 1000 W/m^2
Imp_0 = 8.29; % [A] Module maximum power point current at 25degC and 1000 W/m^2
Voc_0 = 38; % [V] Module open circuit voltage at 25degC and 1000 W/m^2
Isc_0 = 8.85;  % [A] Module short circuit current at 25degC and 1000 W/m^2
beta_Voc = -0.36/100; %[K^-1] Temperature coefficient of open circuit voltage
alpha_Isc = 0.05/100; %[K^-1] Current coefficient of open circuit voltage
Ns = 60; %Number of cells
n_d = 1.26; %Ideality factor diode

% Parameter derivation according to reference paper
delta_0 = (1-beta_Voc*T_0)/(50.1-alpha_Isc*T_0);
w_0 = lambertw(exp(1/delta_0+1));
a_0 = delta_0*Voc_0;
Rs_0 = (a_0*(w_0-1)-Vmp_0)/Imp_0;
Rsh_0 = a_0*(w_0-1)/(Isc_0*(1-1/w_0)-Imp_0);
Iph_0 = (1+Rs_0/Rsh_0)*Isc_0;
Is_0 = Iph_0*exp(-1/delta_0);
alpha_Imp = alpha_Isc + (beta_Voc - 1/T_0)/(w_0 - 1);
beta_Vmp = Voc_0/Vmp_0 * ( beta_Voc/(1+delta_0) + (delta_0*(w_0-1)-1/(1+delta_0))/T_0 );
epsilon_0 = delta_0/(1+delta_0)*Voc_0/Vmp_0;
epsilon_1 = delta_0*(w_0-1)*Voc_0/Vmp_0 - 1;

%% MAIN BUCK CONVERTER

fs = 25000; % Switching frequency [Hz]
eta_buck = 0.9; % Efficiency

% L and C derivation:
P_mps = Vmp_0*Imp_0*9;
V_mps = Vmp_0*3;
I_mps = Imp_0*3;

P_mpw = 0.05*P_mps;
V_mpw = 0.9*V_mps;
I_mpw = P_mpw/V_mpw;

R_mps = V_mps/I_mps;
R_mpw = V_mpw/I_mpw;

Ro = R_mps/2.5;

D_mps = sqrt(Ro/R_mps);
D_mpw = sqrt(Ro/R_mpw);

V_os = D_mps*V_mps;
I_os = V_os/Ro;

V_ow = D_mpw*V_mpw;
I_ow = V_ow/Ro;

d_V_is = 0.002*V_mps;
d_V_os = 0.002*V_os;
d_I_is = 0.2*I_mps;
d_I_os = 0.4*I_os;

d_V_iw = 0.002*V_mpw;
d_V_ow = 0.002*V_ow;
d_I_iw = 0.2*I_mpw;
d_I_ow = 0.4*I_ow;

L_buck = round((2*V_mps*D_mps*(1-D_mpw))/(d_I_ow*fs),3,"significant");
Ci_buck = round(d_I_is/(d_I_iw*fs),3,"significant");
Co_buck = round(d_I_os/(32*d_V_ow*fs),3,"significant");

%% 12V BUCK CONVERTER

fs_buck_small = fs; %Switching frequency
eta_buck_small = 0.87; %Efficiency
I_max_small = 40; %[A] %40 -> According to datasheet of selected converter

%Input
V_in_buck_small = 53; % [V] Usual
V_out_ref_buck_small = 12; % [V] Desired

Dref = V_out_ref_buck_small/V_in_buck_small;

V_ripple_buck_small = 0.002;
I_ripple_buck_small = 0.2;
%____________________________

Delta_i_L_small = I_ripple_buck_small*I_max_small; 

D_buck_small = V_out_ref_buck_small/(56*eta_buck_small);
L_buck_small = round(V_out_ref_buck_small*(V_in_buck_small-V_out_ref_buck_small)/(Delta_i_L_small*fs_buck_small*52),1,"significant");
C_buck_small = round(Delta_i_L_small/(8*fs_buck_small*V_ripple_buck_small*V_out_ref_buck_small),1,"significant");

% Buck converter controller parameters
Ki_buck_small = 0.25; %After tuning
Kp_buck_small = 0.15; %After tuning

%% AUXILIARY LOADS

%Current intended as the sum of all the current request. 
% It may also be a time series object if data on load usage are provided

I_load_12V = 20; %[A]

%% BUCK-BOOST CONVERTER

eta_bb = 0.9; %Efficiency

Vin_bb = 53; % Usual value [V] 
fsw = fs; % Switching frequency [Hz]

Vout_ref_bb = 56; % Desired [V] 
D_bb = Vout_ref_bb/(Vout_ref_bb+Vin_bb);

Io_bb = 160; %[A] max
i_L_bb = Io_bb/(1-D_bb);

Delta_i_L_bb = 0.2*i_L_bb;
Delta_v_bb = 0.002*Vout_ref_bb;

L_bb = round(Vin_bb*D_bb/(Delta_i_L_bb*fsw),1,"significant");
C_bb = round(Io_bb*D_bb/(Delta_v_bb*fsw),1,"significant");

%% WEATHER AND LOAD, n_days, End_Time

% Obtain data on ambient conditions for the selected time interval
% and create a random load request through load_and_weather function
[t_load,load,time_weather,irr_v,temp_v,n_days] = load_and_weather(first_day,last_day,pause_weight);

End_Time = 24*3600*n_days; %Fed to the simulink model call function

%Useful plot calculations
t_max = max(time_weather*3600);
t_plot = strcat(int2str((repmat((0:23),1,n_days))'),{':00'})';
t_plot(1,n_days*24+1)={'00:00'};

%% RUN

h = waitbar(0,'Please wait...'); %To visualize simulation time advance
myoptions = simset('SrcWorkspace','current','DstWorkspace','current','Solver','ode23t'); %Selection of solver

tic
out=sim(model_name,End_Time,myoptions); %Simulation run
time_elapsed = toc; %[s]

close(h)
save(save_name) %The current workspace values are stored to a .mat file

%% WEATHER AND LOAD PLOTS

if plot_weather == 1
    figure(45)
    plot(t_load,load)
    xlabel('Hour')
    ylabel('[m/s]')
    title('Velocity load')
    ylim([0 5.5])
    xlim([0 t_max/3600])
    xticks((0:1:t_max/3600));
    xticklabels(t_plot)

    figure(46)
    yyaxis left
    plot(time_weather,temp_v,'Color',ocra)
    ylabel('Temperature [degC]')
    set(gca,'ycolor',ocra)

    yyaxis right
    plot(time_weather,irr_v,'Color',orange)
    ylabel('Irradiance [W/m^2]')
    ylim([0 1100])
    set(gca,'ycolor',orange)
    legend('Temperature','Irradiance')
    xlabel('Hour')
    xlim([0 t_max/3600])
    xticks((0:1:t_max/3600));
    xticklabels(t_plot)
end

%% RESISTANCE AND PROPULSIVE POWER PLOTS

if plot_vel_drag_pow == 1

figure(51)

    yyaxis left
    plot(vel_drag_pow(:,1),vel_drag_pow(:,2),'Color',blue,'LineStyle','-')
    hold on
    xlabel('Boat velocity [m/s]')
    ylabel('[N]')
    set(gca,'ycolor',blue)

    yyaxis right
    plot(vel_drag_pow(:,1),vel_drag_pow(:,3),'Color',red,'LineStyle','-')
    ylabel('[W]')
    set(gca,'ycolor',red)

    legend('Resistance','Propulsive power')
end

%% MOTOR - PROPELLER PLOTS

if plot_propeller == 1
  
    figure(24)
    subplot(2,2,1)
    plot(v_s_and_Q_prop(:,1),v_s_and_Q_prop(:,2))
    xlabel('Boat velocity [m/s]')
    ylabel('[Nm]')
    title('Propeller torque')

    figure(24)
    subplot(2,2,2)
    plot(v_s_and_np_rpm(:,1),v_s_and_np_rpm(:,2))
    xlabel('Boat velocity [m/s]')
    ylabel('[rpm]')
    title('Propeller angular speed')

    figure(24)
    subplot(2,2,3)
    plot(v_s_and_P_prop(:,1),v_s_and_P_prop(:,2)/1000)
    xlabel('Boat velocity [m/s]')
    ylabel('[kW]')
    title('Propeller load')

    figure(24)
    subplot(2,2,4)
    plot(v_s_and_eta0(:,1),v_s_and_eta0(:,2))
    xlabel('Boat velocity [m/s]')
    ylabel('[-]')
    title('Propeller open water efficiency')

    v_s_plot = 5; %Design speed

    ii = round(v_s_plot*2)+1;

    J_int_str = num2str(J_int_v(ii),3);
    np_rpm_str = num2str(np_rpm_v(ii),4);
    Q_prop_str = num2str(Q_prop_v(ii),3);
    P_prop_str = num2str(P_prop_v(ii)/1000,3);
    v_s_str = num2str(v_s_plot,2);
    title_1 = ['Open water propeller diagram',',',' v_{ship} = ',v_s_str,' [m/s]'];
    subtitle=['J_{match} = ',J_int_str,', ','n_{prop} = ',np_rpm_str, ' [rpm]',', ','Q_{prop} = ',Q_prop_str, ' [Nm]', ',',' P_{prop} = ',P_prop_str, ' [kW]'];

    figure(28)
    H(4)=plot(J_v,Kt_ship_mat(:,ii),'k');
    hold on
    H(1)=plot(J_v,Kt_v,'Color',blue);
    H(2)=plot(J_v,10*Kq_v,'Color',orange);
    H(3)=plot(J_v,eta0_v,'Color',ocra);
    H(5)=plot([J_min_v(ii) J_min_v(ii)],[0 0.9],'Color',grey,'LineStyle','-.','LineWidth',1);
    H(6)=plot([J_int_v(ii) J_int_v(ii)],[0 0.9],'Color',green,'LineStyle','-.','LineWidth',1);
    plot(J_int_v(ii),Kt_int_v(ii),'Color',green,'MarkerSize',30,'Marker','.')

    ylim([0 0.8])
    title(title_1,subtitle)
    xlabel('Advance ratio J')
    legend(H,'K_t','10K_q','\eta_0','Kt_{ship}','J_{min}','J_{match}','location','westoutside')

    v_s_plot = 2; %Any other value in between [0:0.5:5] can be chosen

    ii = round(v_s_plot*2);

    J_int_str = num2str(J_int_v(ii),3);
    np_rpm_str = num2str(np_rpm_v(ii),3);
    Q_prop_str = num2str(Q_prop_v(ii),3);
    P_prop_str = num2str(P_prop_v(ii)/1000,3);
    v_s_str = num2str(v_s_plot,2);
    title_1 = ['Open water propeller diagram',',',' v_{ship} = ',v_s_str,' [m/s]'];
    subtitle=['J_{match} = ',J_int_str,', ','n_{prop} = ',np_rpm_str, ' [rpm]',', ','Q_{prop} = ',Q_prop_str, ' [Nm]', ',',' P_{prop} = ',P_prop_str, ' [kW]'];

    figure(29)
    H(4)=plot(J_v,Kt_ship_mat(:,ii),'k');
    hold on
    H(1)=plot(J_v,Kt_v,'Color',blue);
    H(2)=plot(J_v,10*Kq_v,'Color',orange);
    H(3)=plot(J_v,eta0_v,'Color',ocra);
    H(5)=plot([J_min_v(ii) J_min_v(ii)],[0 0.9],'Color',grey,'LineStyle','-.','LineWidth',1);
    H(6)=plot([J_int_v(ii) J_int_v(ii)],[0 0.9],'Color',green,'LineStyle','-.','LineWidth',1);
    plot(J_int_v(ii),Kt_int_v(ii),'Color',green,'MarkerSize',30,'Marker','.')
    ylim([0 0.8])
    title(title_1,subtitle)
    xlabel('Advance ratio J')
    legend(H,'K_t','10K_q','\eta_0','Kt_{ship}','J_{min}','J_{match}','location','westoutside')

end

%% MODEL PLOTS

if plot_model == 1

    % Comparison between ideal and actual pv power
    figure(1)

    yyaxis left
    plot(out.P_PV_ideal.Time,out.P_PV_ideal.Data,'LineWidth',6,'Color',blue)
    hold on
    plot(out.P_PV.Time,out.P_PV.Data,'Color',ocra,'LineStyle','-')
    ylabel('[W]')
    set(gca,'ycolor','k')

    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)

    %xlabel('Time [Hour]')
    title('PV Power')
    legend('MPP','Actual','SOC battery')
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)

    % Overview of main power output and input involved (PV,battery,loads)
    figure(2)

    yyaxis left
    plot(out.P_loads.Time,medfilt1(out.P_loads.Data,6),'Color',orange,'LineStyle','-','LineWidth',6)
    hold on
    plot(out.P_charge.Time,-medfilt1(out.P_charge.Data,6),'Color',l_brown,'LineStyle','-','LineWidth',3)
    plot(out.P_buck.Time,medfilt1(out.P_buck.Data,6),'Color',blue,'LineStyle','-','LineWidth',3)

    plot(out.P_battery.Time,medfilt1(out.P_battery.Data,6),'Color',green,'LineStyle','-','LineWidth',3)

    %xlabel('Time [Hour]')
    ylabel('[W]')
    ylim([-1700 14000])
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    set(gca,'ycolor','k')
    ylim([-1700 9000])


    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)

    title('Power distribution')
    legend('Motor + auxiliaries','Shore charger','Buck converter','Battery','SOC battery')

    % Effective and requested speed
    figure(3)
    
    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    hold on
    plot([0 86400*n_days],[5 5],'LineStyle','-','Color',[grey,0.7],'LineWidth',1)
    plot([0 86400*n_days],[10 10],'LineStyle','-','Color',[grey,0.7],'LineWidth',1)
    plot([0 86400*n_days],[15 15],'LineStyle','-','Color',[grey,0.7],'LineWidth',1)
    plot([0 86400*n_days],[20 20],'LineStyle','-','Color',[grey,0.7],'LineWidth',1)
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)

    yyaxis left
    plot(out.Vel_requested.Time,out.Vel_requested.Data,'Color',blue,'LineStyle','-','LineWidth',6)
    hold on
    plot(out.Vel_effective.Time,out.Vel_effective.Data,'Color',ocra,'LineStyle','-')
    ylim([0 5.5])
    %xlabel('Time [Hour]')
    ylabel('[m/s]')
    set(gca,'ycolor','k')

    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    title('Boat velocity')
    legend('Requested','Actual','SOC battery','Low SOC [%]')

    % Motor and load speed
    figure(4)

    yyaxis left
    plot(out.w_load.Time,medfilt1(out.w_load.Data,6)*30/pi,'Color',blue,'LineStyle','-','LineWidth',6)
    hold on
    plot(out.w_motor.Time,medfilt1(out.w_motor.Data,6)*30/pi,'Color',ocra,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[rpm]')
    set(gca,'ycolor','k')

    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    title('Load and motor rotational speed')
    legend('Load','Motor','SOC battery')

    % Motor and load torque
    figure(5)

    yyaxis left
    plot(out.Q_load.Time,medfilt1(out.Q_load.Data,6),'Color',blue,'LineStyle','-','LineWidth',6)
    hold on
    plot(out.Q_motor.Time,medfilt1(out.Q_motor.Data,6),'Color',ocra,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[Nm]')
    set(gca,'ycolor','k')

    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)

    title('Load and motor torque')
    legend('Load','Motor','SOC battery')

    % Relevant voltages
    figure(6)

    yyaxis left
    plot(out.V_battery.Time,medfilt1(out.V_battery.Data,6),'Color',green,'LineStyle','-')
    hold on
    plot(out.V_oc_battery.Time,medfilt1(out.V_oc_battery.Data,6),'Color',orange,'LineStyle','-')
    plot(out.V_12V.Time,medfilt1(out.V_12V.Data,6),'Color',purple,'LineStyle','-')
    plot(out.V_motor_drive.Time,medfilt1(out.V_motor_drive.Data,6),'Color',red,'LineStyle','-')
    plot(out.V_PV.Time,medfilt1(out.V_PV.Data,6),'Color',ocra,'LineStyle','-','Marker','none')
    set(gca,'ycolor','k')
    ylabel('[V]')

    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    %xlabel('Time [Hour]')
    title('Relevant Voltages')
    legend('Battery voltage','Open circuit battery voltage','12V loads voltage','Motor drive voltage','PV array voltage','SOC battery')

    % Relevant currents
    figure(7)

    yyaxis left
    plot(out.I_motor.Time,medfilt1(out.I_motor.Data,6),'Color',[orange,0.8],'LineStyle','-','LineWidth',5)
    hold on
    plot(out.I_in_12V.Time,medfilt1(out.I_in_12V.Data,6),'Color',purple,'LineStyle','-','LineWidth',5)
    plot(out.I_PV.Time,medfilt1(out.I_PV.Data,6),'Color',ocra,'LineStyle','-')
    plot(out.I_out_buck.Time,medfilt1(out.I_out_buck.Data,6),'Color',blue,'LineStyle','-','LineWidth',2,'Marker','none')
    plot(out.I_battery.Time,medfilt1(out.I_battery.Data,6),'Color',green,'LineStyle','-','Marker','none','LineWidth',2)
    plot([0 max(out.I_battery.Time)],[-30 -30],'Color',red,'LineStyle','-.','LineWidth',1,'Marker','none')
    set(gca,'ycolor','k')
    ylabel('[A]')
    ylim([-35 160])

    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)

    %xlabel('Time [Hour]')
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    title('Relevant Currents')
    legend('Motor system input current','12V loads input current','PV array current','Buck converter output current','Battery current','-30 [A]','SOC battery')

    % Motor power distribution
    figure(8)
    yyaxis left
    plot(out.P_el_motor.Time,medfilt1(out.P_el_motor.Data,6),'Color',orange,'LineStyle','-')
    hold on
    plot(out.P_mechanical.Time,medfilt1(out.P_mechanical.Data,6),'Color',blue,'LineStyle','-')
    plot(out.P_propulsive.Time,medfilt1(out.P_propulsive.Data,6),'Color',green,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[W]')
    ylim([0 9000])
    set(gca,'ycolor','k')

    yyaxis right
    fill_curves(out.eta_electric.Time,zeros(length(out.SOC.Time),1),medfilt1(out.eta_electric.Data,6)*100,red,0.1);
    hold on
    fill_curves(out.eta_propulsive.Time,zeros(length(out.SOC.Time),1),medfilt1(out.eta_propulsive.Data,6)*100,blue,0.15);
    fill_curves(out.eta_tot.Time,zeros(length(out.SOC.Time),1),medfilt1(out.eta_tot.Data,6)*100,green,0.2);
    ylabel('[%]')
    %xlabel('Time [Hour]')
    ylim([0 100])
    set(gca,'ycolor','k')
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    title('Motor power distribution')
    legend('Electrical input','Mechanical output','Propulsive output','\eta_{electrical}','\eta_{propulsive}','\eta_{tot}')

    % Relevan duty cycles
    figure(9)

    yyaxis left
    plot(out.D_MPPT.Time,medfilt1(out.D_MPPT.Data,6),'Color',[green,0.6],'LineStyle','-','LineWidth',5)
    hold on
    plot(out.D_max_current.Time,medfilt1(out.D_max_current.Data,6),'Color',[orange,0.8],'LineStyle','-','LineWidth',5)
    plot(out.D_to_follow_load.Time,medfilt1(out.D_to_follow_load.Data,6),'Color',[ocra,0.6],'LineStyle','-','LineWidth',5)
    plot(out.D.Time,medfilt1(out.D.Data,6),'Color',blue,'LineStyle','-.','LineWidth',2)
    ylim([0.38 0.72])
    ylabel('[-]')
    set(gca,'ycolor','k')

    yyaxis right
    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.1);
    ylabel('[%]')
    ylim([0 100])
    set(gca,'ycolor',l_blue)

    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    %xlabel('Time [Hour]')
    title('Relevant Duty Cycles')
    legend('D_{MPPT}','D_{max current}','D_{load following}','D','SOC battery')

    % Energy shortcomings: DC MOTOR

    eta_total_str = num2str(eta_total);
    figure(10)
    yyaxis left
    plot(out.P_not_av.Time,medfilt1(out.P_not_av.Data/1000/eta_total,6),'Color',orange,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[kW]')
    set(gca,'ycolor',orange)
    ylim([ 0 max(out.P_not_av.Data/1000/eta_total)+1])

    yyaxis right
    plot(out.E_not_av.Time,medfilt1(out.E_not_av.Data/1000/3600/eta_total,6),'Color',ocra,'LineStyle','-')

    ylabel('[kWh]')
    set(gca,'ycolor',ocra)
    ylim([0 max(out.E_not_av.Data/1000/3600/eta_total)+1])
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)

    subtitle = [num2str(max(out.E_not_av.Data/1000/3600/eta_total),'%.0f'),' [kWh]',' @ \eta_{tot} = ',eta_total_str,' [-]'];
    title('DC motor: Energy shortcomings',subtitle)
    legend('Motor power not available','Energy not available')

    % Battery energy and cycles
    figure(11)

    n_cycles = max(out.E_battery_abs.Data/1000/3600)/E_total_battery;
    n_cycles_str = num2str(round(n_cycles));
    str_title = ['Number of cycles = ',n_cycles_str];

    yyaxis left
    plot(out.E_battery.Time,medfilt1(out.E_battery.Data/1000/3600,6),'Color',orange,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[kWh]')
    set(gca,'ycolor',orange)

    yyaxis right
    plot(out.E_battery_abs.Time,medfilt1(out.E_battery_abs.Data/1000/3600,6),'Color',ocra,'LineStyle','-')

    ylabel('[kWh]')
    set(gca,'ycolor',ocra)
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)

    ylim([ 0 max(out.E_battery_abs.Data/1000/3600)+1])

    title('Battery energy and cycles',str_title)
    legend('Energy exchange','Cumulated energy')

    % Potential PV energy lost
    figure(12)

    yyaxis left
    plot(out.P_PV_loss.Time,medfilt1(out.P_PV_loss.Data/1000,6),'Color',orange,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[kW]')
    set(gca,'ycolor',orange)

    yyaxis right
    plot(out.E_PV_loss.Time,medfilt1(out.E_PV_loss.Data/1000/3600,6),'Color',ocra,'LineStyle','-')

    ylabel('[kWh]')
    set(gca,'ycolor',ocra)
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    subtitle = [num2str(max(out.E_PV_loss.Data/1000/3600),'%.0f'),' [kWh]'];
    title('Potential PV energy lost',subtitle)
    legend('Power','Energy')

    % Energy shortcomings: 12V grid
    figure(13)
    yyaxis left
    plot(out.P_12V_not_av.Time,medfilt1(out.P_12V_not_av.Data/1000,6),'Color',purple,'LineStyle','-')
    %xlabel('Time [Hour]')
    ylabel('[kW]')
    set(gca,'ycolor',purple)
    ylim([ 0 max(out.P_not_av.Data/1000/eta_total)+1])
    yyaxis right
    plot(out.E_12V_not_av.Time,medfilt1(out.E_12V_not_av.Data/1000/3600,6),'Color',ocra,'LineStyle','-')

    ylabel('[kWh]')
    set(gca,'ycolor',ocra)
    ylim([ 0 max(out.E_12V_not_av.Data/1000/3600)+1])
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    subtitle = [num2str(max(out.E_12V_not_av.Data/1000/3600),'%.0f'),' [kWh]'];

    title('Energy shortcomings - 12V grid',subtitle)
    legend('12V power not available','Energy not available')

    % Energy balance

    E_PV_ideal = max(out.E_PV_ideal.Data)/1000/3600;
    E_PV = max(out.E_PV.Data)/1000/3600;
    E_buck = max(out.E_buck.Data)/1000/3600;
    E_PV_loss = max(out.E_PV_loss.Data)/1000/3600;
    E_12V = max(out.E_12V.Data)/1000/3600;
    E_el_motor = max(out.E_el_motor.Data)/1000/3600;
    E_buck_directly = max(out.E_buck_directly.Data)/1000/3600;
    E_buck_stored = max(out.E_buck_stored.Data)/1000/3600;
    E_battery_discharge = max(out.E_battery_discharge.Data)/1000/3600;
    E_battery_shore_charge = max(out.E_battery_shore_charge.Data)/1000/3600;
    E_not_av_motor = max(out.E_not_av.Data/1000/3600/eta_total);
    E_not_av_12V = max(out.E_12V_not_av.Data/1000/3600);

    figure(14)

    X1 = categorical({'PV - Buck energy'});
    X2 = categorical({'Energy source'});
    X3 = categorical({'Energy requested'});

    Y1_1 = E_PV_ideal;
    Y1_2 = [E_PV E_PV_loss];
    Y1_3 = [E_buck E_PV-E_buck];
    Y2 = [E_buck_directly E_buck_stored E_battery_shore_charge];
    Y3 = [E_el_motor E_12V E_not_av_motor+E_not_av_12V];

    H1 = bar(X1,Y1_1,0.8, 'FaceColor','flat');
    set(H1, 'FaceColor', 'Flat')
    H1(1).CData = l_green;
    hold on

    H2 = bar(X1,Y1_2,0.7,'stacked','FaceColor','flat');
    H2(1).CData = beige;
    H2(2).CData = red;

    H3 = bar(X1,Y1_3,0.6,'stacked','FaceColor','flat');
    H3(1).CData = ocra;
    H3(2).CData = brown;

    H4=bar(X2,Y2,'stacked', 'FaceColor','flat');
    H4(1).CData = green;
    H4(2).CData = blue;
    H4(3).CData = l_brown;

    H5=bar(X3,Y3,'stacked', 'FaceColor','flat');
    H5(1).CData = orange;
    H5(2).CData = purple;
    H5(3).CData = [0 0 0];

    title('Energy Balance')
    ylim([0 1.3*max([sum(Y1_1),sum(Y2),sum(Y3)])])
    legend([H1,H2,H3],'PV ideal energy','PV energy','Potential PV energy lost','Buck converter energy','PV energy lost in conversion','Location','NorthWest')
    ylabel('[kWh]')
    a=axes('position',get(gca,'position'),'visible','off');
    legend(a,H4,'Buck energy directly used','Buck energy stored','Shore charger','Location','North');
    b=axes('position',get(gca,'position'),'visible','off');
    legend(b,H5,'DC motor','12V grid','Energy shortcoming');

    % SOC average
    SOC_time = out.SOC_av_h.Time(2:end)/3600;
    SOC_av = out.SOC_av_h.Data(2:end);

    count_on =0;
    count_SOC_on =0;
    count_off=0;
    count_SOC_off =0;

    for zz=1:n_days
        for ii=1:24

            if (SOC_time(ii+(zz-1)*24)>7+(zz-1)*24 && SOC_time(ii+(zz-1)*24)<=22+(zz-1)*24)
                count_on = count_on+1;
                count_SOC_on = count_SOC_on+SOC_av(ii+(zz-1)*24);
            else
                count_off = count_off+1;
                count_SOC_off = count_SOC_off+SOC_av(ii+(zz-1)*24);
            end
        end
    end


    SOC_on_av = count_SOC_on/count_on;
    SOC_off_av = count_SOC_off/count_off;

    SOC_on_av_str = num2str(SOC_on_av,3);
    SOC_off_av_str = num2str(SOC_off_av,3);

    subtitle=['SOC_{mean} during the day = ',SOC_on_av_str,' [%], ',...
        'SOC_{mean} during the night charge = ',SOC_off_av_str,' [%], '];

    figure(15)

    fill_curves(out.SOC.Time,zeros(length(out.SOC.Time),1),out.SOC.Data,blue,0.3);
    xlim([0 t_max])
    xticks((0:3600:t_max));
    xticklabels(t_plot)
    hold on
    H1 = bar(out.SOC_av_h.Time-3600/2,out.SOC_av_h.Data,0.5, 'FaceColor',ocra);
    plot([0 t_max],[SOC_on_av SOC_on_av],'Color',orange,'LineStyle','--','LineWidth',2)
    plot([0 t_max],[SOC_off_av SOC_off_av],'Color',l_brown,'LineStyle','-.','LineWidth',2)
    ylabel('[%]')
    ylim([0 100])
    %xlabel('Time [Hour]')
    title('Battery State of Charge [SOC]',subtitle)
    legend('SOC battery','Hourly average','SOC_{mean} day','SOC_{mean} charge')

end


%% PLOTS DISTRIBUTION

% 1 Comparison between ideal and actual PV power
% 2 Overview of main power output and input involved (PV,battery,loads,shore charger)
% 3 Effective and requested boat velocity
% 4 Load and motor rotational speed
% 5 Load and motor torque
% 6 Relevant voltages
% 7 Relevant currents
% 8 Motor power distribution
% 9 Relevant duty cycles
% 10 Energy shortcomings DC motor
% 11 Battery energy and cycles
% 12 Potential PV energy lost
% 13 Energy shortcomings 12V grid
% 14 Energy balance
% 15 SOC and SOC average

% 24 Propeller torque load
% 24 Propeller angular speed
% 24 Propeller load
% 24 Propeller open water efficiency
% 28 Open water diagram v = 2 [m/s]
% 29 Open water diagram v = 5 [m/s]
% 45 Weather
% 46 Load
% 51 Velocity Drag and Power

distFig('Only',24)
distFig('Position','E','Only',[28 29])
distFig('Position','E','Only',[25 27])
distFig('Position','SW','Only',3)
distFig('Position','W','Only',[1 2])
distFig('Position','E','Only',[6 7])
distFig('Position','SE','Only',9)
distFig('Position','E','Only',[3 8])
distFig('Position','E','Only',[7 9])
distFig('Position','W','Only',[10,11])
distFig('Position','E','Only',[12,13])
distFig('Position','E','Only',[45 46])
distFig('Position','E','Only',[4 5])
distFig('Only',14)
distFig('Only',15)
