% Assumes file mtm_ndbc20042005.mat loaded....

% -------------------------------------------------------------------------
% Plot solar, wind, wave power running averages
% (ndbc data every 1 hour)

iii = find(~isnan(ndbc_rad_swpower1));
solar = ndbc_rad_swpower1(iii);
time = ndbc_rad_time(iii);

ll = length(time);  % length of time series
dd = 24*10/2;  % ten day interval
rave = smooth(solar, 24*10);  % compute average over 10 days
plot(time(dd:ll-dd), rave(dd:ll-dd), '.r')
datetick('x', 2);
xlabel('Time');
ylabel('Solar Power Delivered')
title('Solar(r), wind(g), wave(b) Power 10 Day Running Average');
hold on

% Plot wind power running averages

iii = ~isnan(ndbc_wind_power);
wind = ndbc_wind_power(iii);
time = ndbc_wind_time(iii);
rave = smooth(wind, 24*10);  % compute average over 10 days
ll = length(time);  % length of time series
plot(time(dd:ll-dd), rave(dd:ll-dd), '.g')

% Do the same with waves from NDBC...
% Plot wave power running averages 

iii = (~isnan(ndbc_wave_power));
wave = ndbc_wave_power(iii);
time = ndbc_wave_time(iii);
rave = smooth(wave, 24*10);  % compute average over 10 days
ll = length(time);  % length of time series
plot(time(dd:ll-dd), rave(dd:ll-dd), '.b')


% -------------------------------------------------------------------------
% 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')