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:20;



%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
Cp = g./w;  %Phase speed
Cg = Cp/2;  %Group speed
Lp = 2*pi*g./(w.^2);  %Wavelength

%Quantities that depend on wind speed and frequency range  wmin < w < wmax 
E_bar = zeros(1,length(U));
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);
E_bar(i) = rho*g*trapz(w,S);  %Average energy contained in this spectrum  (W/m^2);
end





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