% EMTH171
% Case Study 2
% North Island
% Frances Power

clear
clc
close all


% Constants
SI_area = 350; % km^2
SI_min_h = 402; % m
SI_max_h = 410; % m
SI_Gen_h = 0; % m
SI_ave_f = 593; % m^3s^-1
SI_MaxGen = 3590; % MW

NI_area = 620; % km^2
NI_min_h = 355.85; % m
NI_max_h = 357.25; % m
NI_Gen_h = 80; % m
NI_ave_f = 345; % m^3s^-1
NI_MaxGen = 1870; % MW

K = 1.55; % m
L = 300; % m
density = 998; % kgm^-3
g = 9.81; % ms^-2


% Arrays
tArray = 1:1:8760;
h_NI = [];
h_SI = [];
V_NI_spill = [];
V_SI_spill = [];
%initial variables
h_NI_0 = 356.55; % m, NI middle water height
h_SI_0 = 406; % m, SI middle water height
V_NI_spill_0 = 345; % average spill rate (m^3s^-1)
V_SI_spill_0 = 593; % average spill rate (m^3s^-1)
% assigning first variable
h_NI(1) = h_NI_0;
h_SI(1) = h_SI_0;
V_NI_spill(1) = V_NI_spill_0;
V_SI_spill(1) = V_SI_spill_0;

windCapacity = 4130; % MW
P_NI_wind = windCapacity* 2/3;
P_SI_wind = windCapacity* 1/3;


for ii = 1:1:length(tArray)
    mean = 5000; % hours
    standardDeviation = 1000; % hours
    P_NI_demand = 4065 + 1.4e6*normpdf(tArray(ii), mean, standardDeviation);
    P_SI_demand = 1940;
    F_NI_in = 345 + 73* sin(2*pi*(tArray(ii) - 3624)./8760);
    F_SI_in = 593 - 183* sin(2*pi*(tArray(ii) - 2320)./8760);
    F_NI_gen = F_NI_in;
    P_NI_hydro = 0.9* F_NI_gen.* density.* g*(h_NI_0 - NI_Gen_h)./10.^6;
    CF = 0.41 + 0.12*sin(2*pi*(tArray(ii)-5660)/8760);
    P_geo = 1525;
    NI_total = P_NI_hydro + P_geo + P_NI_wind;
    P_hvdc = P_NI_demand - P_geo - P_NI_wind - P_NI_hydro;
    P_SI_hydro = P_SI_demand + P_hvdc - P_SI_wind;
    F_SI_gen = P_SI_hydro * 10^6./ (.9* density* g* (h_SI_0 - SI_Gen_h));



end