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 = .1*vol*rho;  %Reaction Mass (kg)
mB = rho*vol-mR;  %Buoy Mass (kg)

kk = mR*g*.0258/d;  % Spring Coefficient N/m  
bb = 1000;  %Damping Coefficient N/m/s

kd0 = kk*d/(mR*g);

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


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));
%xi = (A*Fex)./(-w.^2.*(a22+mB)+i*w.*(b22)+rho*g*S);
zeta = xi.*(kk+i*w*bb)./(-mR*w.^2+kk+i*w*bb);

%Time domain results
for j = 1:length(k);  %Frequency selection
t(j,:) = 0:max(T)/200:2*max(T);
eta(j,:) = cos(w(j)*t(j,:));  %Incident wave
y1(j,:) = real(xi(j)*exp(i*w(j)*t(j,:)))/A;
y2(j,:) = real(zeta(j)*exp(i*w(j)*t(j,:)))/A;
end





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




figure
plot(kd,abs(xi)/A,kd,abs(zeta)/A,kd,a22/(rho*vol),kd,b22./(rho*vol*w));
xlabel('w^2*d/g');
legend('Buoy RAO','Mass RAO','a22/(rho*V)','b22/(rho*V*w)');

figure
plot(kd,180*angle(xi)/pi,kd,180*angle(zeta)/pi);
xlabel('w^2*d/g');
ylabel('phase angle relative to incident wave');
legend('Buoy Motion','Mass Motion');

figure
plot(kd,Pw,kd,Pd,kd,Pr,kd,Pe)
xlabel('w^2*d/g');
ylabel('Power (W)');
legend('Power In Wave impacting buoy','Power Dissipated in Damper','Power Dissipated in Radiated Waves','Power in Wave Exciting Force');

figure
plot(kd,Pw,kd,Pd,kd,Pr,kd,Pe)
xlabel('w^2*d/g');
ylabel('Power (W)');



figure
plot(kd,100*eff);
xlabel('w^2*d/g');
ylabel('efficiency (%)');

init_time = round(length(T)/2);

fh=figure;
plot(t(init_time,:),eta(init_time,:),t(init_time,:),y1(init_time,:),t(init_time/2,:),y2(init_time,:))
xlabel('time (s)');
ylabel('y/A');
legend('Incident Wave','Buoy Motion','Mass Motion');
minx = 0;
maxx = max(T);
miny = -max([max(eta) max(y1) max(y2)]);
maxy = -miny;
axis([minx maxx miny maxy]);
text(0.02*maxx,.90*maxy,['Period = ' num2str(T(init_time)) ' sec // foo']);
text(0.02*maxx,.80*maxy,['efficiency = ' num2str(100*eff(init_time)) ' %']);

sh = uicontrol(fh,'Style','slider',...
                'Max',length(T),'Min',1,'Value',init_time,...
                'SliderStep',[1/length(T) 5/length(T)],...
                'Position',[350 2 200 20],...
                'Callback', 'j = round(get(sh,''Value''));plot(t(j,:),eta(j,:),t(j,:),y1(j,:),t(j,:),y2(j,:));axis([minx maxx miny maxy]);xlabel(''time (s)'');ylabel(''y/A'');legend(''Incident Wave'',''Buoy Motion'',''Mass Motion'');text(0.02*maxx,.90*maxy,[''Period = '' num2str(T(j)) '' sec'']);text(0.02*maxx,.80*maxy,[''efficiency = '' num2str(100*eff(j)) '' %'']);');