% Assumes file mtm_ndbc20042005.mat loaded....
 
% -------------------------------------------------------------------------
% Plot solar and wind power running averages
 
rave = smooth(mtm2_solarpower, 24*6*30); % compute average over 30 days
ll = length(mtm2_time); % length of time series
dd = 24*6*30/2; % interval of a month
plot(mtm2_time(dd:ll-dd)-mtm2_time(1), rave(dd:ll-dd))
datetick('x', 2);
xlabel('Time');
ylabel('Solar Power Delivered')
title('Solar Power 3,10 and 30 Day Running Average');
hold on
 
rave = smooth(mtm2_solarpower, 24*6*10); % compute average over 10 days
plot(mtm2_time(dd:ll-dd)-mtm2_time(1), rave(dd:ll-dd), 'g')
 
rave = smooth(mtm2_solarpower, 24*6*3); % compute average over 3 days
plot(mtm2_time(dd:ll-dd)-mtm2_time(1), rave(dd:ll-dd), 'r')
 
% Plot wind power running averages
 
rave = smooth(mtm2_windpower, 24*6*30); % compute average over 30 days
ll = length(mtm2_time); % length of time series
dd = 24*6*30/2; % interval of a month
plot(mtm2_time(dd:ll-dd)-mtm2_time(1), rave(dd:ll-dd))
datetick('x', 2);
xlabel('Time');
ylabel('Wind Power Delivered')
title('Wind Power 3,10 and 30 Day Running Average');
hold on
 
rave = smooth(mtm2_windpower, 24*6*10); % compute average over 10 days
plot(mtm2_time(dd:ll-dd)-mtm2_time(1), rave(dd:ll-dd), 'g')
 
rave = smooth(mtm2_windpower, 24*6*3); % compute average over 3 days
plot(mtm2_time(dd:ll-dd)-mtm2_time(1), rave(dd:ll-dd), 'r')
 
% -------------------------------------------------------------------------
% Compute battery-load-generation performance envelope
 
power = mtm2_windpower;
ave_power = mean(power);
battery = battery_predict(mtm2_time, power, ave_power/2, ave_power*12);
frac_battery_flat = length(find(battery == 0))/length(battery);
 
power = mtm2_solarpower;
ave_power = mean(power);
battery = battery_predict(mtm2_time, power, ave_power/2, ave_power*12);
frac_battery_flat = length(find(battery == 0))/length(battery);
 
power = mtm2_solarpower/mean(mtm2_solarpower) + mtm2_windpower/mean(mtm2_windpower);
ave_power = mean(power);
battery = battery_predict(mtm2_time, power, ave_power/2, ave_power*12);
frac_battery_flat = length(find(battery == 0))/length(battery);
 
% -------------------------------------------------------------------------
% Do a range of load and battery capacity conditions
power = (mtm2_solarpower/mean(mtm2_solarpower) + mtm2_windpower/mean(mtm2_windpower))/2;
ave_power = mean(power);
 
size = 4;
clear load capacity frac_empty battery
for i = 1:size*2
    for j = 1:size
        generation(i,j) = 1.3^(i-1)/2; % power generated normalized against load
        capacity(i,j) = 12 * 2^j;
 
        ii = i+(j-1)*4;
        battery(ii,:) = battery_predict(mtm2_time, power*generation(i,j), 1, capacity(i,j));
 
        frac_empty(i,j) = length(find(battery(ii,:) == 0))/length(battery(ii,:));
    end;
end
 
pcolor(capacity, generation, frac_empty)
shading interp
colorbar