clear all
close all
clc

%% INITIALISATION

% Add paths to needed packages (wafo, toporoll, etc.)
addpath(genpath(fullfile(['C:\Users\luigir\Desktop\backup\scripts\toolbox\stochastic generation toolbox\wafo25'])));  
addpath(genpath(fullfile(['C:\Users\luigir\Desktop\matlab\data analysis\Curve analysis\Toolboxes\toporoll']))); 
addpath(genpath(fullfile(['C:\Users\luigir\Desktop\wafo-wafo-v2017.1\wafo-wafo-v2017.1\wafo'])));
addpath(genpath(fullfile(['C:\Users\luigir\Desktop\backup\scripts\toolbox\stochastic generation toolbox\auxiliary functions'])));

% Save to designation directory
designationDir = 'C:\Users\luigir\Desktop\sOC for Anandh\generated_dOCs';

% Load an sOC
load('C:\Users\luigir\Desktop\sOC for Anandh\mean_sOC\sOCVästraGötaland.mat')
StatisticalOC = statOC;

%% SET INITIAL CONDITIONS

numberOfDays = 5;
seasonIndex = 4;      

time = 0:1:24*numberOfDays-1;
numSteps = length(time);

% Generate day of the year
if seasonIndex == 1
    population =-30:60;
    rndDay = randsample(population,1);
    if rndDay < 0
        rndDay = -rndDay + 330;
    end
elseif seasonIndex == 2
    population = 60:150;
    rndDay = randsample(population,1);
elseif seasonIndex == 3
    population = 150:240;
    rndDay = randsample(population,1);
elseif seasonIndex == 4
    population = 240:330;
    rndDay = randsample(population,1);
end

%% Generate temperature series

% Generate stochastic component
MdlTemperature = StatisticalOC.Weather.Temperature{seasonIndex,1}.ARIMA;
stocasticTemperature = simulate(MdlTemperature, numSteps);

% Generate deterministic component
meanT = StatisticalOC.Weather.Temperature{seasonIndex,1}.meanT;
T1 = StatisticalOC.Weather.Temperature{seasonIndex,1}.T1;
omega1 = StatisticalOC.Weather.Temperature{seasonIndex,1}.omega1;
phi_T1 = StatisticalOC.Weather.Temperature{seasonIndex,1}.phi_T1;
T2 = StatisticalOC.Weather.Temperature{seasonIndex,1}.T2;
omega2 = StatisticalOC.Weather.Temperature{seasonIndex,1}.omega2;
phi_T2 = StatisticalOC.Weather.Temperature{seasonIndex,1}.phi_T2;

n = 0;
starting_hour = rndDay*24;
for i=1:numSteps
    n = n+1;
    total_hour = starting_hour + time(i);
    determinisitcTemperature = meanT + T1*sin(omega1*total_hour+phi_T1) ...
                               + T2*sin(omega2*total_hour + phi_T2);
end
Temperature = stocasticTemperature + determinisitcTemperature;


%% Generate humidity series

% Generate stochastic component
MdlHumidity = StatisticalOC.Weather.Humidity{seasonIndex,1}.ARIMA;
stocasticHumidity = simulate(MdlHumidity, numSteps);

% Generate deterministic component
meanPsi = StatisticalOC.Weather.Humidity{seasonIndex,1}.meanPsi;
psi1 = StatisticalOC.Weather.Humidity{seasonIndex,1}.psi1;
omega1 = StatisticalOC.Weather.Humidity{seasonIndex,1}.omega1;
phi_psi1 = StatisticalOC.Weather.Humidity{seasonIndex,1}.phi_psi1;
psi2 = StatisticalOC.Weather.Humidity{seasonIndex,1}.psi2;
omega2 = StatisticalOC.Weather.Humidity{seasonIndex,1}.omega2;
phi_psi2 = StatisticalOC.Weather.Humidity{seasonIndex,1}.phi_psi2;

n = 0;
starting_hour = rndDay*24;
for i=1:numSteps
    n = n+1;
    total_hour = starting_hour + time(i);
    determinisitcHumidity = meanPsi + psi1*sin(omega1*total_hour+phi_psi1) ...
                               + psi2*sin(omega2*total_hour + phi_psi2);
end
Humidity = stocasticHumidity + determinisitcHumidity;
for i=1:length(Humidity)
    if Humidity(i) > 100
    Humidity(i) = 100;
    end
end


%% Generate rain series

% n = 0;
% mcRainAutumn = dtmc(StatisticalOC.Weather.Rain{seasonIndex,1}.PP);
% rng(1); % For reproducibility
% XAutumn = simulate(mcRainAutumn,numSteps);
% 
% for i=1:numSteps
%     n = n+1;
%     if XAutumn(i) == 2
%         rndRain(n) = gamrnd(alpha_hp,beta_hp);
%     elseif XAutumn(i) == 1
%         rndRain(n) = 0;
%     end
% end
% OperatingCycle.Weather.PrecipitationLevel.Data = [];
% OperatingCycle.Weather.PrecipitationLevel.Data(1,:) = time;
% OperatingCycle.Weather.PrecipitationLevel.Data(2,:) = rndRain;
n = 0;
P = StatisticalOC.Weather.Rain{seasonIndex,1}.PP;
alpha_hp = StatisticalOC.Weather.Rain{seasonIndex,1}.alpha_hp;
beta_hp = StatisticalOC.Weather.Rain{seasonIndex,1}.beta_hp;
zsim = mcsim(P,numSteps,[]);
for i=1:numSteps
    n = n+1;
    if zsim(i) == 2
        rndRain(n) = gamrnd(alpha_hp,beta_hp);
    elseif zsim(i) == 1
        rndRain(n) = 0;
    end
end

%% Generate atmospheric pressure series

MdlPressure = StatisticalOC.Weather.Pressure{seasonIndex,1}.ARIMA;
Pressure = simulate(MdlPressure,numSteps);


%% Generate wind speed and direction series

MdlWind = StatisticalOC.Weather.Wind{seasonIndex,1}.VARM;
mu_w = StatisticalOC.Weather.Wind{seasonIndex,1}.mu_w;
Wind = simulate(MdlWind, numSteps);

WindSpeed = max(0,Wind(:,1));
WindDirection = 2*pi*(normcdf(Wind(:,2))-0.5)+mu_w;
for i =1:length(WindDirection)
    if WindDirection(i) < 0
        WindDirection(i) = 2*pi+WindDirection(i);
    elseif WindDirection(i) > 2*pi
        WindDirection(i) = 2*pi-WindDirection(i);
    end        
end

%%
save('WeatherSeriesAutumn', 'Temperature','Humidity','Pressure','WindSpeed','WindDirection');