close all
clear all

T = 9;  %Specified Wave period

g = 9.81;

WaveAmplitudes = 1:5;


NperLb = 4.45;

Cd = 1.0;
Ac = .456;  %m^2.  Cross section of 30" diameter assembly
rho = 1025;

m = 450+150;  %ESP effective mass.

Kl = 5000;  %N/m  ex. 100m of 1/4" Nilspin
Kl = 7.9*9.81*.8;  %N/m  3/4" chain wet weight per meter for chain-catenary case
Kl = 50:20:8880;
Ku = 270;  %N/m, large WHOI stretch hose prestretched to epsilon = 0.5 (i.e. 10m hose under 300lbs stretches to 15 meters and that behavior is pretty linear up to epsilon = 1, or 20 meters for 10 meter hose.  Stiffens dramatically beyond that.


%Compute and plot at a range of wave amplitudes and wave upper spring
%constants, at specifed wave period.

w = 2*pi/T;
f= 1/T;

figure('Position',[1000 246 1028 1252]);
for j = 1:length(WaveAmplitudes)
A = WaveAmplitudes(j);  %Wave amplitude

b = Cd*0.5*rho*Ac*A*w/sqrt(2);  %Linearized damping
zeta = b./(2*sqrt(m*(Ku+Kl)));

k = 1;
vel = -A*w(k):.1:A*w(k);
figure('Position',[1529 702 987 606]);
plot(vel,-Cd*.5*rho*Ac*vel.*abs(vel)/NperLb,vel,-b(k)*vel/NperLb);
xlabel('Velocity (m/s)')
ylabel('Damping Force (lbs)');
title(['Damping Forces at Wave Amp = ' num2str(A) ' m and Wave Period = ' num2str(T(k)) ' sec']); 
grid on
lh = legend({'Quadratic Damping Force',['Linearized Damping based on RMS speed = $A \omega / \sqrt{2} = ' num2str(A*w(k)/sqrt(2),3) '$ m/s']},'interpreter','latex');


for jj = 1:length(Kl)
  [Z(jj) Fu(jj) Fl(jj)] = ComplexAmplitude(A,w,Ku,Kl(jj),b,m);
end

figure
plot(T,(m*A*w.^2)/NperLb)
xlabel('Wave Period (s)');
ylabel('Force to cause mass to follow waves (lbs per wave amp)');

leg_str{j} = ['Wave Amp = ' num2str(A)];

figure(1)
subplot(3,1,1)
plot(Kl/NperLb,abs(Z))
hold all
subplot(3,1,2)
plot(Kl/NperLb,abs(Fu)/NperLb)
hold all
subplot(3,1,3)
plot(Kl/NperLb,abs(Fl)/NperLb)
hold all

end

figure(1)
subplot(3,1,1)
xlabel('Lower Spring Constant (lb/m) ');
ylabel('Motion Amplitude (m)');
legend(leg_str,'Location','NorthEast');
title(['Effective mass = ' num2str(m) ' kg  ' 'Wave Period = ' num2str(T) ' s  ' 'Upper Spring Constant = ' num2str(Ku/NperLb,3) ' lbs/m']);
subplot(3,1,2)
xlabel('Lower Spring Constant (lb/m) ');
ylabel('Upper Spring Force Amplitude (lbs)');
legend(leg_str,'Location','NorthEast');
subplot(3,1,3)
xlabel('Lower Spring Constant (lb/m) ');
ylabel('Lower Spring Force Amplitude (lbs)');
legend(leg_str,'Location','NorthEast');




