clear all
close all

rho = 1025;
g = 9.8;
d = 2;  %meters

Cd = 2.0;
Cm = 1.0;



A_inf = 3.5;  %1 meter wave Amplitude
T = 18.5;  %10 second period
w = 2*pi/T;



k_inf = w^2/g;
lambda_inf = 2*pi/k_inf;
P_inf = .25*rho*g*A_inf^2*w/k_inf;
Vp_inf = w/k_inf;
Vg_inf = .5*w/k_inf;
Fi_inf = 
Fd_inf = 


h = A_inf:1:250;  %depth range in meters

k = zeros(1,length(h));
for i = 1:length(h)
k(i) = wave_num(h(i),w,g);
end

lambda = 2*pi./k;
A = A_inf*sqrt(k./(2*k_inf))./sqrt(.5+k.*h./sinh(2*k.*h));
P = .5*rho*g*A.^2.*(.5+k.*h./sinh(2*k.*h))*w./k;
Vp = w./k;
Vg = (.5+k.*h./sinh(2*k.*h)).*Vp;
Fi = -rho*Cm*pi*d^2*A.*w^2./(16*k);
Fd = 0.5*rho*Cd*d*(A./2).^2*w^2.*(sinh(2.*k.*h)/4+k.*h./2)./(4.*k.*sinh(k.*h).^2);

figure
plot(h,A./A_inf,h,lambda/lambda_inf);
legend('A/A_inf','lambda/lambda_inf');

figure
plot(h,Vp./Vp_inf,h,Vg./Vg_inf);
legend('Vp/Vp_inf','Vg/Vg_inf');

figure
plot(h,Fi./Fi(end


%R = ((Cd*h)/(Cm*d))*(sinh(2*k*H)/4+k*H/2)/sinh(k*H).^2;



%Fi = -rho*Cm*pi*d^2*h*w^2/(8*k);

%Fd = 0.5*rho*Cd*d*h^2*w^2*(sinh(2*k*H)/4+k*H/2)/(4*k*sinh(k*H)^2);