close all
clear all


a = 1;
D = 2*a;
rho = 1025;
g = 9.81;


drafts = [4 2 1 .6 .4 .2];  % 0.1];

for i = 1:length(drafts)

  d = drafts(i);
    V = pi*d*a^2;
  
  filename = ['t_cyl_d' num2str(d)]

  FD_coefs{i} = Readwdra3dResults([filename '.out']);
    
  w_nd = FD_coefs{i}.w*sqrt(D/g);  %The same for all drafts so no need to index
  w_nd_smooth = w_nd(1):(w_nd(end)-w_nd(1))/100:w_nd(end);

  am_nd(i,:) = real(FD_coefs{i}.F{3,3})/(rho*V);
  am_nd_smooth(i,:) = interp1(w_nd,am_nd(i,:),w_nd_smooth,'pchip');
  
  dmp_nd(i,:) = imag(FD_coefs{i}.F{3,3})./(rho*V*FD_coefs{i}.w);
  dmp_nd_smooth(i,:) = interp1(w_nd,dmp_nd(i,:),w_nd_smooth,'pchip');
  
  real_exct_nd(i,:) = real(FD_coefs{i}.F{7,3})./(rho*V*g);
  real_exct_nd_smooth(i,:) = interp1(w_nd,real_exct_nd(i,:),w_nd_smooth,'pchip');
    
  imag_exct_nd(i,:) = -imag(FD_coefs{i}.F{7,3})./(rho*V*g);
  imag_exct_nd_smooth(i,:) = interp1(w_nd,imag_exct_nd(i,:),w_nd_smooth,'pchip');
  
  
  
  roll_am_nd(i,:) = real(FD_coefs{i}.F{5,5})/(rho*V*D^2);
  roll_am_nd_smooth(i,:) = interp1(w_nd,roll_am_nd(i,:),w_nd_smooth,'pchip');
  
  roll_dmp_nd(i,:) = imag(FD_coefs{i}.F{5,5})./(rho*V*D^2*FD_coefs{i}.w);
  roll_dmp_nd_smooth(i,:) = interp1(w_nd,roll_dmp_nd(i,:),w_nd_smooth,'pchip');
   
  roll_real_exct_nd(i,:) = real(FD_coefs{i}.F{7,5})./(rho*V*g*D);
  roll_real_exct_nd_smooth(i,:) = interp1(w_nd,roll_real_exct_nd(i,:),w_nd_smooth,'pchip');
  
  roll_imag_exct_nd(i,:) = -imag(FD_coefs{i}.F{7,5})./(rho*V*g*D);
  roll_imag_exct_nd_smooth(i,:) = interp1(w_nd,roll_imag_exct_nd(i,:),w_nd_smooth,'pchip');
  
  leg_str{i} = ['$d/D = ' num2str(d/D) ' $ \hspace{.1in}'];
  
  
  %Determine natural heave frequency and amplitude for each d/D ratio.
  S = 0.25*pi*D^2;
  m = rho*V;
  
  [vm] = Heave_VM(w_nd_smooth*sqrt(g/D),m,rho,g,S,FD_coefs{i}.w,am_nd(i,:)*rho*V);
  
  figure
  plot(w_nd_smooth,vm)
  
  
  fun = @(x) Heave_VM(x,m,rho,g,S,FD_coefs{i}.w,am_nd(i,:)*rho*V);    % function of x alone
  x = fzero(fun,3);
  
  NatFreq_H(i) = x*sqrt(D/g);
  
  %HeaveExt = interp1(FD_coefs{i}.w,FD_coefs{i}.F{7,3},NatFreq_H(i)*sqrt(g/D));
  HeaveExt = interp1(FD_coefs{i}.w,(real_exct_nd(i,:)+sqrt(-1)*imag_exct_nd(i,:))*rho*V*g,NatFreq_H(i)*sqrt(g/D));
  HeaveDmp = interp1(FD_coefs{i}.w,dmp_nd(i,:)*rho*V.*FD_coefs{i}.w,NatFreq_H(i)*sqrt(g/D));
    
  MagRes_H(i) = abs(HeaveExt./(sqrt(-1)*NatFreq_H(i)*sqrt(g/D)*HeaveDmp));
  
end


%Plot Heave Resonances
figure('Position',[1000 402 950 1096],'PaperOrientation','portrait','PaperSize',[7.5 6],'PaperPositionMode','manual','PaperPosition',[-.75 0 9 6]);
for i = 1:length(NatFreq_H)
  plot(NatFreq_H(i),MagRes_H(i),'o');
  hold all
  leg_str{i} = ['d/D = ' num2str(drafts(i)/D)];
end
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Heave Resonance Frequency $\left( \omega_n \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Heave Resonance Magnitude $\left( \frac{\left| Z_{33} \right| }{A} \right)$');
set(alh,'Interpreter','Latex');
grid on

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













%Plot Heave Hydrodynamic Coefficients

figure('Position',[1000 402 950 1096],'PaperOrientation','portrait','PaperSize',[7.5 10],'PaperPositionMode','manual','PaperPosition',[-.5 0 8.5 10]);
subplot(2,1,1)
plot(w_nd_smooth,am_nd_smooth);
%hold on
%plot(w_nd,am_nd,'x');
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Added-Mass $\left( \frac{a_{33}}{\rho \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
%th = title('Heave Added Mass for Cylindrical Buoy');
%set(th,'Interpreter','Latex');
grid on
xlim([0 5]);

subplot(2,1,2)
plot(w_nd_smooth,dmp_nd_smooth);
%hold on
%plot(w_nd,dmp_nd,'x')
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Damping $\left( \frac{b_{33}}{\rho \omega \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
%th = title('Heave Damping for Cylindrical Buoy');
%set(th,'Interpreter','Latex');
grid on
xlim([0 5]);


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

figure('Position',[1000 402 950 1096],'PaperOrientation','portrait','PaperSize',[7.5 10],'PaperPositionMode','manual','PaperPosition',[-.5 0 8.5 10]);
subplot(2,1,1)
plot(w_nd_smooth,real_exct_nd_smooth);
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Real Exciting Force $\left( \frac{| X_{33} |}{\rho g \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
%th = title('Heave Exciting Force Magnitude for Cylindrical Buoy');
%set(th,'Interpreter','Latex');
grid on
xlim([0 5]);

subplot(2,1,2)
plot(w_nd_smooth,imag_exct_nd_smooth);
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Imag Exciting Force $\left( \frac{| X_{33} |}{\rho g \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
%th = title('Heave Exciting Force Magnitude for Cylindrical Buoy');
%set(th,'Interpreter','Latex');
grid on
xlim([0 5]);

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




figure('Position',[1000 402 950 1096],'PaperOrientation','portrait','PaperSize',[7.5 10],'PaperPositionMode','manual','PaperPosition',[-.5 0 8.5 10]);
subplot(2,1,1)
plot(w_nd_smooth,roll_am_nd_smooth);
%hold on
%plot(w_nd,roll_am_nd,'x');
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Added-Mass $\left( \frac{a_{55}}{\rho D^2 \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
%th = title('Roll Added Mass for Cylindrical Buoy');
%set(th,'Interpreter','Latex');
grid on
xlim([0 5]);

subplot(2,1,2)
plot(w_nd_smooth,roll_dmp_nd_smooth);
%hold on
%plot(w_nd,roll_dmp_nd,'x')
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Non-dimensional Damping $\left( \frac{b_{55}}{\rho \omega D^2 \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
%th = title('Roll Damping for Cylindrical Buoy');
%set(th,'Interpreter','Latex');
grid on
xlim([0 5]);


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



figure('Position',[1000 402 950 1096],'PaperOrientation','portrait','PaperSize',[7.5 10],'PaperPositionMode','manual','PaperPosition',[-.5 0 8.5 10]);
subplot(2,1,1)
plot(w_nd_smooth,roll_real_exct_nd_smooth);
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Real Exciting Force $\left( \frac{| X_{33} |}{\rho g \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');

grid on
xlim([0 5]);

subplot(2,1,2)
plot(w_nd_smooth,roll_imag_exct_nd_smooth);
lh = legend(leg_str);
set(lh,'Interpreter','Latex');
alh = xlabel('Non-dimensional Frequency $\left( \omega \sqrt{\frac{D}{g}} \right)$');
set(alh,'Interpreter','Latex');
alh = ylabel('Imag Exciting Force $\left( \frac{| X_{33} |}{\rho g \bigtriangledown} \right)$');
set(alh,'Interpreter','Latex');
grid on
xlim([0 5]);

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

