close all
clear all

rho = 1025;
g = 9.81;

c0 = rho*pi*g/4;

cinf = 245;


H = 1:3;
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 = (5000/2.2)/1025;  
%V = 2.5;
f = 1./T;
Pb = c0*V*H*(T.^-1);
Pa = cinf*(H.^2)*(T.^3);


figure('Position',[ 321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(T,min(Pb,Pa),'LineWidth',2.5);
set(gca,'fontsize',18);

xlabel('Wave Period (s)');
ylabel('Power (W)','fontsize',18);
%title(['Upper Bound for ' num2str(V) ' m^3 point absorber in heave - regular waves']);
title(['Upper bound for 5000 lb displacement point absorber in heave - regular waves']);
legend('H = 1m','H = 2m','H = 3m');
grid on

print('-dpdf','BudalUpperBound');
print('-djpeg','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
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper bound for 5000 lb displacement point absorber in heave - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),4) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');
 
print('-dpdf','PM_BudalUpperBound');
print('-djpeg','PM_BudalUpperBound');



%Compute upper bound for each spectrum for WaveGlider
V = (200/2.2)/1025;
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper power bound for WaveGlider - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),3) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');
 
print('-dpdf','PM_WaveGliderBudalUpperBound');
print('-djpeg','PM_WaveGliderBudalUpperBound');






%Compute upper bound for each spectrum for BuoyA2
V = 1.3*1.06;
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper power bound for BuoyA2 - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),4) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');
 


%Compute upper bound for each spectrum for BuoyA2
V = 1.3*1.06;
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper power bound for BuoyA2 - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),4) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');


%Compute upper bound for each spectrum for BuoyA5
V = 1.3*2.2;
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper power bound for BuoyA5 - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),4) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');


%Compute upper bound for each spectrum for BuoyA10
V = 1.3*5.26;
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper power bound for BuoyA10 - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),4) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');


%Compute upper bound for each spectrum for BuoyA10
V = 1.3*10.37;
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
    i
     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',[321         109        1196         779],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,S_Hs,'LineWidth',2.5);
set(gca,'fontsize',18);
title(['Upper power bound for BuoyA20 - PM spectrum']);
for i = 1:length(Hs)
    leg{i} = ['H_s = ' num2str(Hs(i),1) ' m     P_b = ' num2str(Pb(i),4) ' W'];
end
legend(leg);
xlabel('Wave Frequency (Hz)');
ylabel('m^2/Hz');




