clear all
close all


FD_coefs = ReadData('tcyl.out');

T = 2;  %Period of oscillation.
dt = 0.01;


tau = 0:dt:10;
L = ComputeIRF_coefs(FD_coefs.w,FD_coefs.F,tau);

I = zeros(3);
I(1,1) = 1613.25; I(2,2) = 1615.53; I(3,3) = 1613.25;

Buoy = RigidBody;
Buoy = set(Buoy,'mass',3220.13,'S',3.14159,'disp',3.14159,'g',9.81);
Buoy = set(Buoy,'xG',[0 -.25 0]);
Buoy = set(Buoy,'xB',[0 -.5 0]);
Buoy = set(Buoy,'S11',0.79,'S22',0.79);  %m^4
Buoy = set(Buoy,'I',I);
Buoy = set(Buoy,'dt',dt);


tf = 40;

t = 0:dt:tf;

Buoy = set(Buoy,'N',length(t));

wave = IncWave;
wave = set(wave,'T',T);

eta0 = eta(wave,0,t);
%F_E(1,:) = [0 0 0 0 0 0];
RB{1} = Buoy;
for n = 1:length(eta0)-length(tau)
    n
F_E(n,:) = F_Exciting(Buoy,RB,n+.5,L,eta0);
end


%Compute Frequency Domain Forces
jj = 3;
f = 1/T;
a = spline(FD_coefs.w,real(FD_coefs.F{7,jj}),2*pi*f);
b = spline(FD_coefs.w,imag(FD_coefs.F{7,jj}),2*pi*f);

ForceFreq = a*cos(2*pi*f*t)+b*sin(2*pi*f*t);

figure
plotyy(t(1:length(eta0)-length(tau)),eta0(1:length(eta0)-length(tau)),t(1:length(eta0)-length(tau)),F_E(:,jj))
figure
plot(t(1:length(eta0)-length(tau)),F_E(:,jj),t(1:length(eta0)-length(tau)),ForceFreq(1:length(eta0)-length(tau)))
