clear all
close all


FD_coefs = ReadData('box1.out');

F = FD_coefs.F;
w = FD_coefs.w;

figure('Position', [100 100 900 600])
set(gcf,'color','white');


subplot(2,2,1),plot(w,real(F{1,1}),w,real(F{3,3}))
xlabel('w');
ylabel('Added Mass');
legend('Surge','Heave');
set(gca,'FontName','times');

subplot(2,2,2),plot(w,imag(F{1,1}),'x-',w,imag(F{3,3}))
xlabel('w');
ylabel('Damping');
legend('Surge','Heave');
set(gca,'FontName','times');

subplot(2,2,3),plot(w,real(F{7,1}),w,real(F{7,3}));
xlabel('w');
ylabel('Re[Exciting Force]');
legend('Surge','Heave');
set(gca,'FontName','times');

subplot(2,2,4),plot(w,imag(F{7,1}),w,imag(F{7,3}));
xlabel('w');
ylabel('Im[Exciting Force]');
legend('Surge','Heave');
set(gca,'FontName','times');


%Calculation of Response Amplitude Operator as a function of frequency

amp = 1.0;      %wave amplitude (m)
S = 14.4;       %Waterplane area of floating body (m)
rho = 1025;     %Density of seawater (kg/m^3)
g = 9.81;       %acceleration due to gravity (m/s^2)
Vbar = 14.4;    %displacement of floating body (m^3)

stiff = 1000;   %spring stiffness (N/m)
damp  = 200;     %damping coefficient (N/(m/s))

m_fl = 12760;   %mass of floating body
m_sub = 2000;   %mass of submerged square plate.....rho*Vbar = m_fl + m_sub

len = 2.4;      %side of square
add_sub = rho*pi*len.^2/4;% m33 = (pi.rho.len^2)/4

add_fl = real(F{3,3});
dc_fl = imag(F{3,3});

fex_fl = F{7,3};
fex_sub = 0;

% response amplitude operator of Floating Body



for j = 1:80;
 
    
T(j) =2*pi/w(j);


%Steps to calculate linearized drag using damping coefficient
kin_visc = 1.05*(10.^-6); %in m^2/ sec at 20 deg. C
Re(j) = len*abs(amp*w(j))/kin_visc;
DCoef(j) = (2*0.427)/((log10(Re(j))-0.407)^2.64); %from page 2-5 in FD Drag-Hoerner...Cd = 2*Cf

area = len^2;
dc_sub(j)= (0.5)*DCoef(j)*rho*area*abs(amp*w(j));

%Calculate Floating body RAO 
num_1(j) =  fex_sub*(stiff +i*w(j)*damp); 

num_2(j) =  (fex_fl(j))*(m_sub*w(j)*w(j) + w(j)*w(j)*add_sub - i*w(j)*dc_sub(j) + (stiff + i*w(j)*damp));

den_1(j) = ((stiff + i*w(j)*damp).^2);

den_2(j) = (-m_sub*w(j)*w(j) -w(j)*w(j)*add_sub + i*w(j)*dc_sub(j) - (stiff + i*w(j)*damp));
    
den_3(j) = (-w(j)*w(j)*m_fl + rho*g*S -w(j)*w(j)*add_fl(j)+ i*w(j)*dc_fl(j) + (stiff + i*w(j)*damp));

rao_fl(j) = (num_1(j) - num_2(j))/(den_1(j) + den_2(j)*den_3(j));
          
         
          
%Calculate Submerged body RAO


rao_sub(j) = (-w(j)*w(j)*m_fl + rho*g*S -w(j)*w(j)*add_fl(j)+ i*w(j)*dc_fl(j) + stiff + i*w(j)*damp)*rao_fl(j) - fex_fl(j);

          
  rao_sub(j) = rao_sub(j)/(stiff + i*w(j)*damp); 
  
  delv(j) = amp*(rao_fl(j) - rao_sub(j));
  
  avg_power(j) = 0.5*damp*w(j)*w(j)*abs((delv(j)).^2);
  
end

arr = zeros(80,4);
arr(:,1) = w;
arr(:,2) = T;
arr(:,3) = abs(rao_fl);
arr(:,4) = abs(rao_sub);


avg_power = avg_power';

figure('Position', [100 100 900 600])
set(gcf,'color','white');

subplot(2,1,1),plot(w, abs(rao_fl),'Color','b','linewidth',2.5);

hold on
plot(w, abs(rao_sub),'Color','r','linewidth',2.5);legend('Floating body','Submerged plate')
%xlabel('w (rad/s)','FontName','times', 'FontSize',18);
ylabel('Heave RAO (m/m)','FontName','times', 'FontSize',22,'FontWeight','bold');
set(gca,'FontName','times', 'FontSize',20);
xlim([0 4])

subplot(2,1,2),plot(w,avg_power,'linewidth',2.5);
xlabel('w (rad/s)','FontName','times', 'FontSize',22,'FontWeight','bold');
ylabel('Avg Power over 1 cycle(W)','FontName','times', 'FontSize',22,'FontWeight','bold');
set(gca,'FontName','times', 'FontSize',20);
xlim([0 4])












