
name = 'BuoyA2';

results_filenm = create_filenm(['FlatPlate_' name],Hs,Ls,Lc,dt,tf,b)

rho = 1025;
g = 9.81;

N = tf/dt;      %Number of timesteps to compute

%Set up "Buoy"
Buoy = RigidBody;
Buoy = set(Buoy,'name','Surface Buoy');
Buoy = set(Buoy,'dt',dt,'N',N);  %Set timesteps
Buoy = set(Buoy,'dof',6);  %With 6 DOF.
Buoy = set(Buoy,'Fixed',0);  %Set free to move
Buoy = set(Buoy,'mass', 738.8,'S',2.01,'disp',1.06,'g',g);
Buoy = set(Buoy,'x0',[0 0 0]);  %Equilibrium Position in global coordinates.
Buoy = set(Buoy,'xG',[0 0 -0.67]);
Buoy = set(Buoy,'xB',[0 0 -0.30]);
Buoy = set(Buoy,'S11',0.32,'S22',.32);  %S11 = S22 = (-1)(-1)(3.14159) m^4
I = zeros(3);
I(1,1) = 546.6; I(2,2) = 546.6; I(3,3) = 102.1;      
Buoy = set(Buoy,'I',I);
clear I;
Buoy = set(Buoy,'Area',[1.1 1.1 2.01]);
Buoy = set(Buoy,'Cd',[21.5 1.5 1]);  %Large Horizontal Drag component to keep it from drifting due to start-up transient.
Buoy = set(Buoy,'AreaW',[.64 .64 0]);
Buoy = set(Buoy,'CdW',[1.5 1.5 0]);  %Cd on pi*r^2 from Hoerner Figure 14, pp5-8.  Needs refinement

%Read Frequency Domain Hydrodynamic Data from file.

filename = ['WAMIT_Results\' name];
disp('Reading frequency domain hydrodynamic coefficients');
FD_coefs = ReadWamitFD(filename);
for i = 1:6 
   for j = 1:6 
      mu_inf(i,j) = rho*FD_coefs.X{i,j}(2); 
   end
end
Buoy = set(Buoy,'mu_inf',mu_inf);
clear mu_inf;

%Compute Impulse Response Functions from Frequency Domain Results

disp('Computing impulse response functions'); 
L = ComputeWamitIRF_coefs(filename);


%Set up Plate
S = Lc*Ls;
thickness = 0.0164;  %(For flat plate of thickness equal to mean thickness of NACA0015 foil with 40cm chord)
displ = Lc*Ls*thickness;  %m^3 
wet_weight = 347.7*g;  %Newtons
mass = wet_weight/g+displ*rho;  %kg 
Plate = RigidBody;
Plate = set(Plate,'name','Plate');
Plate = set(Plate,'x0',[0  0 -20]);  %Equilibrium Position in global coordinates.
Plate = set(Plate,'dt',dt,'N',N);  %Set timesteps
Plate = set(Plate,'dof',6);  %With 6 DOF.
Plate = set(Plate,'Fixed',0);  %Set free to move
Plate = set(Plate,'mass',mass,'S',0,'disp',displ,'g',9.81); 
Plate = set(Plate,'xG',[0.0 0 0]);
Plate = set(Plate,'xB',[0.0 0 0]);
Plate = set(Plate,'S11',0.0,'S22',0.0);  %m^4


I = zeros(3);
I(1,1) = 35.5; I(2,2) = 1.88; I(3,3) = 37.3;
Plate = set(Plate,'I',I);
clear I;
Plate = set(Plate,'Area',[0  0  S]);
Plate = set(Plate,'Cd',[0  0  1.2]);

mu_inf = zeros(6,6);     
mu_inf(1,1) = Ls*pi*rho*(thickness/2)^2; %m11 
mu_inf(3,3) = Ls*pi*rho*(Lc/2)^2;   %added mass for 2-d ellipse plate.  Newman pp 145
mu_inf(5,5) = (1/8)*Ls*pi*rho*((Lc/2)^2-(thickness/2)^2)^2;                     %m55

Plate = set(Plate,'mu_inf',mu_inf);


%Compute incident wave profile

t = get(Buoy,'t');

Type = 0;  %Type refers monochromatic(Type = 1) or PM spectrum (Type = 0) wave
A = Hs/2;
wind_spd = sqrt(Hs*g/0.21); %Wind appropriate to wave-height in PM spectrum
wave = IncWave;
wave = set(wave,'A',A);
wave = set(wave,'Type',Type);
wave = set(wave,'beta',-1);

disp('Computing incident wave profile');
t1= dt:dt:L{7,1}.breaks(end);
eta0 = eta(wave,0,[t t(end)+t1]);

%Add forces to Buoy
Buoy = add_Force_fcn(Buoy,@F_HydroStatic);
Buoy = add_Force_fcn(Buoy,@F_Weight);
Buoy = add_Force_fcn(Buoy,@(body,RB,n) Mem_Radiation(body,RB,n,L));
Buoy = add_Force_fcn(Buoy,@(body,RB,n) F_Exciting(body,RB,n,L,eta0));
Buoy = add_Force_fcn(Buoy,@(body,RB,n) F_Drag_Inc(body,RB,n,'Surface Buoy',wave));
%Buoy = add_Force_fcn(Buoy,@(body,RB,n) F_Wind(body,RB,n,[-wind_spd 0 0]));

%Add forces to Plate
Plate = add_Force_fcn(Plate,@F_HydroStatic);
Plate = add_Force_fcn(Plate,@F_Weight);
Plate = add_Force_fcn(Plate,@(body,RB,n) F_Drag_Inc(body,RB,n,'Surface Buoy',0));
%Add forces that act on two bodies
p_aB = [0.0  0.0  -1.25];  %Tether attachment point on Buoy (local CoordSys).
p_aF = [0.0  0.0  0.0];  %Tether attachment point on Plate (local CoordSys).

k = wet_weight/1;  %Spring Elongates 1m to support static wet weight
[Spring_Force_Coeffs] = LinearSpringForce(k);
%[Damping_Force_Coeffs] = LinearDampingForce(b);
[Damping_Force_Coeffs] = QuadraticDampingForce(b);

Buoy = add_Force_fcn(Buoy,@(body,RB,n) Hydraulic_tether(body,RB,n,'Plate',p_aB,p_aF,Spring_Force_Coeffs,Damping_Force_Coeffs));
Plate = add_Force_fcn(Plate,@(body,RB,n) Hydraulic_tether(body,RB,n,'Surface Buoy',p_aF,p_aB,Spring_Force_Coeffs,Damping_Force_Coeffs));



%Set Initial Conditions  %Offset position from equilibrium position.
Buoy = set_pos(Buoy,[0 0 eta0(1) 0 0 0],1);
Plate = set_pos(Plate,[0 0 0 0 0 0],1);

RB{1} = Buoy;
RB{2} = Plate;

%Solve for motion
RB = Solve_RK4(RB);


Buoy = RB{1};
Plate = RB{2};


save([results_filenm '.mat'],'Buoy','Plate','rho','g','dt','tf','N','eta0','wave','Lc','Ls','p_aF','p_aB','Spring_Force_Coeffs','Damping_Force_Coeffs');










