% Assumes file mtm_ndbc20042005.mat loaded....

% -------------------------------------------------------------------------
% Plot solar and wind power running averages

solar = mtm2_solarpower/mean(mtm2_solarpower);
rave = smooth(solar, 24*6*3);  % compute average over 3 days
ll = length(mtm2_time);  % length of time series
dd = 24*6*30/2;  % interval of a month
plot(mtm2_time(dd:ll-dd), rave(dd:ll-dd), 'r')
datetick('x', 2);
xlabel('Time');
ylabel('Solar Power Delivered')
title('Solar Power 3,10 and 30 Day Running Average');
hold on

rave = smooth(solar, 24*6*10);  % compute average over 10 days
plot(mtm2_time(dd:ll-dd), rave(dd:ll-dd), '.k')

rave = smooth(solar, 24*6*30);  % compute average over 3 days
plot(mtm2_time(dd:ll-dd), rave(dd:ll-dd), 'x')

% Plot wind power running averages

wind = mtm2_windpower/mean(mtm2_windpower);
rave = smooth(wind, 24*6*3);  % compute average over 3 days
ll = length(mtm2_time);  % length of time series
dd = 24*6*30/2;  % interval of a month
plot(mtm2_time(dd:ll-dd), rave(dd:ll-dd), 'r')
datetick('x', 2);
xlabel('Time');
ylabel('Wind Power Delivered')
title('Wind Power 3,10 and 30 Day Running Average');
hold on

rave = smooth(wind, 24*6*10);  % compute average over 10 days
plot(mtm2_time(dd:ll-dd), rave(dd:ll-dd), '.k')

rave = smooth(wind, 24*6*30);  % compute average over 30 days
plot(mtm2_time(dd:ll-dd), rave(dd:ll-dd), 'x')

% -------------------------------------------------------------------------
% Do the same with waves from NDBC...first get overlapping data
ii = find(ndbc_wave_time > mtm2_time(1)& ndbc_wave_time < mtm2_time(length(mtm2_time))&~isnan(ndbc_wave_power));
wave = ndbc_wave_power(ii)/mean(ndbc_wave_power(ii));
wt = ndbc_wave_time(ii);

% Plot wave power running averages (ndbc data every 1 hour)

rave = smooth(wave, 24*3);  % compute average over 3 days
ll = length(wt);  % length of time series
dd = 24*30/2;  % interval of a month
plot(wt(dd:ll-dd), rave(dd:ll-dd), 'r')
datetick('x', 2);
xlabel('Time');
ylabel('Wave Power Delivered')
title('Wave Power 3,10 and 30 Day Running Average');
hold on

rave = smooth(wave, 24*10);  % compute average over 10 days
plot(wt(dd:ll-dd), rave(dd:ll-dd), '.k')

rave = smooth(wave, 24*30);  % compute average over 30 days
plot(wt(dd:ll-dd), rave(dd:ll-dd), 'x')

% -------------------------------------------------------------------------
% Do the same with wind from NDBC...first get overlapping data
ii = find(ndbc_wind_time > mtm2_time(1)& ndbc_wave_time < mtm2_time(length(mtm2_time))&~isnan(ndbc_wind_power));
wind = ndbc_wind_power(ii)/mean(ndbc_wind_power(ii));
wt = ndbc_wind_time(ii);

% Plot wind power running averages (ndbc data every 1 hour)

rave = smooth(wind, 24*3);  % compute average over 3 days
ll = length(wt);  % length of time series
dd = 24*30/2;  % interval of a month
plot(wt(dd:ll-dd), rave(dd:ll-dd), 'r')
datetick('x', 2);
xlabel('Time');
ylabel('Wave Power Delivered')
title('NDBC Wind Power 3,10 and 30 Day Running Average');
hold on

rave = smooth(wind, 24*10);  % compute average over 10 days
plot(wt(dd:ll-dd), rave(dd:ll-dd), '.k')

rave = smooth(wind, 24*30);  % compute average over 30 days
plot(wt(dd:ll-dd), rave(dd:ll-dd), 'x')

% -------------------------------------------------------------------------
% Do the same with current from M2...first get overlapping data
ii = find(time > mtm2_time(1)& time < mtm2_time(length(mtm2_time))&~isnan(u_20m));
current = (u_20m(ii).^2 + v_20m(ii).^2).^0.5;
power = 1/2 * 1027 * current.^3;
power = power / mean(power);
ct = time(ii);

% Plot wave power running averages (M2 data every 1 hour)

rave = smooth(power, 24*3);  % compute average over 3 days
ll = length(ct);  % length of time series
dd = 24*30/2;  % interval of a month
plot(ct(dd:ll-dd), rave(dd:ll-dd), 'r')
datetick('x', 2);
xlabel('Time');
ylabel('Current Power (W)')
title('Current Power 3,10 and 30 Day Running Average');
hold on

rave = smooth(power, 24*10);  % compute average over 10 days
plot(ct(dd:ll-dd), rave(dd:ll-dd), '.k')

rave = smooth(power, 24*30);  % compute average over 30 days
plot(ct(dd:ll-dd), rave(dd:ll-dd), 'x')

% -------------------------------------------------------------------------
% Compute battery-load-generation performance envelope

load = 0.25;  % In units of average power...
power = mtm2_windpower;

ave_power = mean(power);
power = power/ave_power;
storage_capacity = 12; % hours of average power ...
battery = battery_predict(mtm2_time, power, load, storage_capacity);
frac_battery_flat = length(find(battery == 0))/length(battery);

plot(mtm2_time, power, 'r'); hold on;
plot(mtm2_time, battery/storage_capacity, 'b')
datetick('x', 2);
xlabel('Time');
ylabel('Power & Stored Energy');
title('Performance under Constant Load');


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);

size = 4;
for i = 1:size
    for j = 1:size
        load = ave_power / 2^(i-1);
        capacity = load * 12 * 2^j;

        battery(i+(j-1)*4,:) = battery_predict(mtm2_time, power, load, capacity);
    end;
end

% ------------------------------
% Compute curves showing fraction of time the battery goes empty as a
% function of storage capacity of battery for for solar, wind, and
% solar+wind

load = 1/3;  % In units of average power...

power = mtm2_solarpower; %----------SUN
ave_power = mean(power);
power = power/ave_power;

for i = 1:4
    storage_capacity(i) = 3*2^((i-1)/2); % hours of average power ...
    battery = battery_predict(mtm2_time, power, load, storage_capacity(i));
    frac_battery_flat_solar(i) = length(find(battery == 0))/length(battery);
end

power = mtm2_windpower; %------------WIND
ave_power = mean(power);
power = power/ave_power;

for i = 1:4
    battery = battery_predict(mtm2_time, power, load, storage_capacity(i));
    frac_battery_flat_wind(i) = length(find(battery == 0))/length(battery);
end

                        %-------------EQUAL SUN AND WIND
power = mtm2_windpower/mean(mtm2_windpower)+mtm2_solarpower/mean(mtm2_solarpower);
ave_power = mean(power);
power = power/ave_power;

for i = 1:4
    battery = battery_predict(mtm2_time, power, load, storage_capacity(i));
    frac_battery_flat_mix(i) = length(find(battery == 0))/length(battery);
end

plot(storage_capacity, frac_battery_flat_mix, 'g'); hold on; 
plot(storage_capacity, frac_battery_flat_wind, 'b')
plot(storage_capacity, frac_battery_flat_solar, 'r')
hold off