load hourlyM1.mat

% Convert to mks...
u_20m = u_20m/100;
v_20m = v_20m/100;

% Compute wind/current speed and a number proportional to power
c_speed = (u_20m.^2 + v_20m.^2).^(1/2);
wind_speed = (wind_u.^2 + wind_v.^2).^(1/2);
wind_speed3 = wind_speed .^ 3;
c_speed3 = c_speed .^ 3;

% Running average of days worth of wind & current speed...
windowSize = 24;
c_speed_ave = filter(ones(1,windowSize)/windowSize,1,c_speed);
wind_speed_ave = filter(ones(1,windowSize)/windowSize,1,wind_speed);
c_speed3_ave = filter(ones(1,windowSize)/windowSize,1,c_speed3);
wind_speed3_ave = filter(ones(1,windowSize)/windowSize,1,wind_speed3);

% ok = find(wind_speed_ave < 1e5);
ok = find(wind_speed_ave < 1e5 & c_speed_ave < 1e5);
wind_u = wind_u(ok);
wind_v = wind_v(ok);
wind_speed = wind_speed(ok);
wind_speed3 = wind_speed3(ok);
wind_speed_ave = wind_speed_ave(ok);
wind_speed3_ave = wind_speed3_ave(ok);
wind_time = time(ok);

% ok = find(c_speed_ave < 1e5);
u_20m = u_20m(ok);
v_20m = v_20m(ok);
c_speed = c_speed(ok);
c_speed3 = c_speed3(ok);
c_speed_ave = c_speed_ave(ok);
c_speed3_ave = c_speed3_ave(ok);
c_time = time(ok);

% Compute power in wind and ocean motion...
wind_power = 1.168 * wind_speed3/2;
c_power = 1027 * c_speed3/2;
%...and 24hr average....
wind_power_ave = 1.168 * wind_speed3_ave/2;
c_power_ave = 1027 * c_speed3_ave/2;

% Plot wind versus ocean power...
plot(c_power_ave, wind_power_ave, '.')

% Plot power
figure; plot(wind_time, wind_power, '.'); hold on; plot(c_time, c_power, 'r+'); hold off;

% [wsa, ix_wsa] = sort(wind_speed_ave);
% [csa, ix_csa] = sort(c_speed_ave);

% Plot the results...
figure; plot(wind_time, wind_speed_ave, '.'); hold on; plot(wind_time, wind_speed3_ave.^(1/3), 'r+');
datetick('x',2)
figure; plot(wind_speed_ave, wind_speed3_ave.^(1/3)./wind_speed_ave, '.');

figure; plot(c_time, c_speed_ave, '.'); hold on; plot(c_time, c_speed3_ave.^(1/3), 'r+');
datetick('x',2)
figure; plot(c_speed_ave, c_speed3_ave.^(1/3)./c_speed_ave, '.');