close all
clear all

alpha = 0.0081;
beta = 0.74;

rho = 1025;  %kg/m^3
g = 9.81;  %m^2/s
df = 0.001;
f = df:df:.3;
w = 2*pi*f;
dw = 2*pi*df;

figure
hold on;


U = 1:30;



%Quantities that depend on wind speed
Hsig = 0.21*(U.^2)/g;  %Significant waveheight.
wp = 0.877*g./U;  %freq of spectrum peak
Tp = 1.14*U;  %Period of spectrum peak
E_barTotal = rho*g*0.00274*(U.^4)/(g^2);  %Average energy of waves in entire spectrum at range of wind speeds.

%Quantities that depend on frequency
Vp = g./w;  %Phase speed
Vg = Vp/2;  %Group speed
L = 2*pi*g./(w.^2);  %Wavelength

%Quantities that depend on wind speed and frequency range  fmin < f < fmax
fmin = 0.0  ; wmin = 2*pi*fmin;
fmax = .4 ; wmax = 2*pi*fmax;
idx = find(w<=wmax & w>=wmin);

E_bar = zeros(1,length(U));
eff = zeros(1,length(w));
eff(idx) = 1;
for i = 1:length(U);
S = alpha*g^2*exp(-beta*(g./(U(i)*w)).^4)./(w.^5);  % P-M Wave spectrum (M^2/freq) 
plot(f,2*pi*S);
plot(f(idx),2*pi*S(idx),'r');
E_bar(i) = rho*g*trapz(w,S.*eff);  %Average energy contained in this spectrum  (W/m^2);
dEdt_bar(i) = rho*g*trapz(w,S.*eff.*Vg);
end



figure
plot(U,E_barTotal,U,E_bar);
xlabel('Wind Speed (m/s)');
ylabel('Wave Energy per m^2  (J/m^2)');


figure
plot(U,dEdt_bar);
xlabel('Wind Speed (m/s)');
ylabel('Wave Power  (W/m)');

figure
plot(Hsig,dEdt_bar,Hsig,rho*g^2*Hsig.^2./(16*wp));
xlabel('Significant Wave Height');
ylabel('Wave Power  (W/m)');

