close all
clear all

rho = 1025;
g = 9.81;

c0 = rho*pi*g/4;

cinf = 245;


H = 1:4;
H = H';
T = 4:20;

%Peak Free Surface Speeds
figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(T,.5*H*2*pi*(T.^-1))
xlabel('Wave Period (s)');
ylabel('(m/s)');
title(['Peak Free-Surface Speed  for Monochromatic Waves']);
legend('H = 1m','H = 2m','H = 3m','H = 4m');
grid on


%Plot of wave energy flux for monochromatic waves

J = rho*g^2*(H/2).^2*T/(8*pi);
figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
semilogy(T,J);
xlabel('Wave Period (s)');
ylabel('Power (W/m)');
title(['Wave-Energy Transport for Monochromatic Waves']);
legend('H = 1m','H = 2m','H = 3m','H = 4m');
grid on

print('-dpdf','WaveEnergyTransport');


%Budal Limt
V = 1.25*(5000/2.2)/1025;  
V = 2.5;
f = 1./T;
Pb = c0*V*H*(T.^-1);
Pa = cinf*(H.^2)*(T.^3);


figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(T,min(Pb,Pa));
xlabel('Wave Period (s)');
ylabel('Power (W)');
title(['Upper Bound for ' num2str(V) ' m^3 point absorber in heave - regular waves']);
legend('H = 1m','H = 2m','H = 3m','H = 4m');
grid on

print('-dpdf','BudalUpperBound');


f = 0:.002:.5;
f = f';
alpha = 8.1e-3;
beta = .74;
w = 2*pi*f;

df = f(2)-f(1);

 Hs = 1:5;
 clear S_Hs;
 for i = 1:length(Hs)
   U19 = sqrt(Hs(i)*g/0.21);
   w0 = g/U19;
   S_Hs(i,:) = 2*pi*(alpha*g^2)*exp(-beta*(w0./w).^4)./(w.^5);
   S_Hs(i,1) = S_Hs(i,2);
 end
 
 figure;
 plot(f,S_Hs);
 xlabel('freq (Hz)');
 ylabel('m^2/Hz');
 title('Pierson-Moskowitz Spectrums');
 legend('H_s = 1 m','H_s = 2 m','H_s = 3 m','H_s = 4 m','H_s = 5m');
 


%Compute upper bound for each spectrum 
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
     Pb(i) = BudalLimit(S_Hs(i,:),f,V);
 % Pow(i,:) = 2*df*S_Hs(i,:).*(f'.^2);
  %Pb(i) = 2*V*c0*sqrt(sum(Pow(i,:)));
end


figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs);
for i = 1:length(Hs)
    leg{i} = ['Hs = ' num2str(Hs(i),1) ' m     Pb = ' num2str(Pb(i),5) ' W'];
end
legend(leg);

%zs = (2*S_Hs*df).^0.5;

%figure;plot(f,2*V*c0*zs(1,:).*f')