close all

A = 1;
g = 9.81;
d = 5.0; %.64;
rho = 1025;
S = pi*d^2/4;
vol = pi*d^3/12; %.35;  %m^3  pi*d^3/12;
kd = k*d;
i = sqrt(-1);

F_b = rho*vol

mR = .5*vol*rho;  %Reaction Mass (kg)
mB = rho*vol-mR;  %Buoy Mass (kg)

Fex = F{7,3};
a22 = real(F{3,3});
b22 = imag(F{3,3});

kk_range = mR*g*.0258/d;  % Spring Coefficient N/m  
bb_range = [30000 40000 50000];  %Damping Coefficient N/m/s

PwResults = cell(1,length(bb_range));
PdResults = cell(1,length(bb_range));
effResults = cell(1,length(bb_range));

for j = 1:length(bb_range)
j
bb = bb_range(j)


xi = (A*Fex)./(-w.^2.*(a22+mB)+i*w.*(b22+bb)+rho*g*S+kk-(kk+i*w*bb).^2./(-mR*w.^2+kk+i*w*bb));
zeta = xi.*(kk+i*w*bb)./(-mR*w.^2+kk+i*w*bb);

%Power Dissipation
Pw = d*0.25*rho*g^2*A^2./w;
Pd = A^2*0.5*bb*w.^2.*abs(xi-zeta).^2;
eff = Pd./Pw;

PwResults{j} = Pw;
PdResults{j} = Pd;
effResults{j} = eff;

end