
name = 'BuoyA20';

results_filenm = create_filenm(['FlappingFoil_' 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', 9038,'S',19.63,'disp',10.37,'g',g);
Buoy = set(Buoy,'x0',[0 0 0]);  %Equilibrium Position in global coordinates.
Buoy = set(Buoy,'xG',[0 0 -0.63]);
Buoy = set(Buoy,'xB',[0 0 -0.30]);
Buoy = set(Buoy,'S11',30.66,'S22',30.66);  %S11 = S22 = (pi*r^4)/4  for circle
I = zeros(3);
I(1,1) = 12541.3; I(2,2) = 12541.3; I(3,3) = 12105.0;      
Buoy = set(Buoy,'I',I);
clear I;
Buoy = set(Buoy,'Area',[3.5 3.5 19.63]);
Buoy = set(Buoy,'Cd',[1.5 1.5 1]);
Buoy = set(Buoy,'AreaW',[2.5 2.5 19.63]);
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 "Foil"  
S = Lc*Ls;
thickness = 0.15*Lc;
displ = Ls*0.0164;  %m^3 (For NACA0015 Foil)
wet_weight = 1589.2*g;  %Newtons
mass = wet_weight/g+displ*rho;  %kg 
Foil = RigidBody;
Foil = set(Foil,'name','Foil');
Foil = set(Foil,'x0',[0  0 -20]);  %Equilibrium Position in global coordinates.
Foil = set(Foil,'dt',dt,'N',N);  %Set timesteps
Foil = set(Foil,'dof',6);  %With 6 DOF.
Foil = set(Foil,'Fixed',0);  %Set free to move
Foil = set(Foil,'mass',mass,'S',0,'disp',displ,'g',9.81); 
Foil = set(Foil,'xG',[0.0 0 0]);
Foil = set(Foil,'xB',[0.0 0 0]);
Foil = set(Foil,'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;
Foil = set(Foil,'I',I);
clear I;
Foil = set(Foil,'Area',[0  0  S]);


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

Foil = set(Foil,'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 Foil
Foil = add_Force_fcn(Foil,@F_HydroStatic);
Foil = add_Force_fcn(Foil,@F_Weight);
kk = 20*Ls;
bb = 7*Ls;
Foil = add_Force_fcn(Foil,@(body,RB,n) F_Wing(body,RB,n,'Surface Buoy',0,kk,bb));

%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 Foil (local CoordSys).

%gamma = 1.4;
%[Spring_Force_Coeffs] = AccumSpringForce(gamma);
%b = 3000;
%[Damping_Force_Coeffs] = DampingForce(b);

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,'Foil',p_aB,p_aF,Spring_Force_Coeffs,Damping_Force_Coeffs));
Foil = add_Force_fcn(Foil,@(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);
Foil = set_pos(Foil,[0 0 0 0 0 0],1);

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

%Solve for motion
RB = Solve_RK4(RB);


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

results_filenm


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










