clear;

data = xlsread('Hourly_electricity_data_2016.xlsx','H2:N8785');

wind_power = data(:,3); % Hourly wind power production, column 3
consumption = data(:,7); % Hourly consumption, column 7
wind_scale_factor = []; % TODO: add correct factors for different cases here

% Initializing matrices for duration curves and ramps
sorted_lmw = zeros(length(consumption),3);
ramps = zeros(length(consumption)-1,3);

% For-loop goes through all three cases

for j = 1:3
    load_minus_wind = ; % TODO: add formula for load minus wind here
  
    % Power duration curve
    sorted_lmw(:,j) = ; % TODO: sort the data to descending order
    
    % Ramp-up/down

    timestep = 1; % 1 hour
    for i = 1:size(ramps,1) % Calculating delta P / delta t
        ramps(i,j) = ; % TODO: insert formula for calculating ramps
    end   

    rbins = 100*floor(min(min(ramps))/100):20:100*ceil(max(max(ramps))/100); % Ramp-up bins
    freq = hist(ramps,rbins); % Frequencies
end

% Plots duration curves

figure(1);
plot(sorted_lmw,'LineWidth',1.5);
colormap jet
xlabel('Hours (h)');
ylabel('Power (MW)');
title('Power duration curve'); 
legend('Load - 1 x wind','Load - 3 x wind','Load - 10 x wind');
axis([0 9000 -4500 15500]);

% Plots histograms

figure(2);
plot(rbins,freq,'LineWidth',1.5);
xlabel('Ramp-up/down (MW/h)');
ylabel('Frequency (number of 1 h timeslots)');
title('Variability analysis');
legend('Load - 1 x wind','Load - 3 x wind','Load - 10 x wind');

% Saves figures as .png files

print(figure(1), '-dpng', '-r400', 'b_duration');
print(figure(2), '-dpng', '-r400', 'b_variability');