close all
clear all

config = 3;

g = 9.81;

WaveAmplitudes = 1:5;


NperLb = 4.45;

d = 24;

Cd = 0.8;
Ac = 0.25*pi*(d*.0254)^2;  %m^2.  Cross section of 24" diameter assembly
rho = 1025;


if(config == 1)
    m = 590;  %ESP effective mass.
    Ku = 200;  %N/m, large WHOI stretch hose prestretched to epsilon = 0.5 (i.e. 10m hose under 225lbs 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.
    Kl = 68;   %N/m  3/4" chain wet weight per meter for chain-catenary case
end

if(config == 2)  
    m = 590;  %ESP effective mass.
    Ku = 100000;  %Stiff connection between buoy and ESP. "Potted Chain"
    Kl = 68;   %N/m  3/4" chain wet weight per meter for chain-catenary case
end

if(config == 3)
    m = 1030;  %ESP effective mass.
    Ku = 200;  %N/m, large WHOI stretch hose prestretched to epsilon = 0.5 (i.e. 10m hose under 225lbs 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.
    Kl = 50000;  %N/m  ex. 100m of 1/4" Nilspin
end


%Compute and plot at a range of wave amplitudes and wave frequencies
T = 3:.1:20;  %Wave Period
w = 2*pi./T;  %Wave Ang. Freq
f = 1./T;  %Wave Freq;
lambda = 2*pi*g./(w.^2);

wRes = sqrt((Ku+Kl)/m);
fRes = wRes/(2*pi)
TRes = 2*pi/wRes
figure('Position',[1000 246 1028 1252],'PaperType','usletter','PaperOrientation','portrait','PaperPositionMode','auto','PaperPosition',[.5 .5 7.5 10.5]);
for j = 1:length(WaveAmplitudes)
A = WaveAmplitudes(j);  %Wave amplitude

Steepness = 2*A./lambda;

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

wResDamped = wRes*(1-2*zeta.^2).^2;
fResDamped = wResDamped/(2*pi)
TResDamped = 2*pi./wResDamped

k = 20;
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');

idx = find(Steepness < 1/7);
jdx = find(Steepness >= 1/7);

figure('Position',[1000 246 1028 1252]);
subplot(2,1,1)
plot(T,zeta,T,ones(length(T),1)/sqrt(2),'--')
xlabel('Excitiation Period (s)');
ylabel('Damping Ratio');
subplot(2,1,2)
plot(T,TResDamped,T,TRes*ones(length(T),1),'--')
xlabel('Excitiation Period (s)');
ylabel('Damped Resonance Period (s)');
legend('Damped Resonance Period (s)',['Undamped Resonance Period = ' num2str(TRes) ' s']);
ylim([min(T) max(T)]);


%Y_RAO = -Ku./(m*w.^2-(Ku+Kl)+i*w.*b);
%Fu = Ku*A*(1-Y_RAO);
%Fl = -Kl*A*Y_RAO;
[Z Fu Fl] = ComplexAmplitude(A,w,Ku,Kl,b,m);

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(4,1,1)
plot(T,abs(Z))
hold all
subplot(4,1,2)
plot(T,(1/g)*abs(Z.*w.^2))
hold all
subplot(4,1,3)
plot(T,abs(Fu)/NperLb)
hold all
subplot(4,1,4)
plot(T,abs(Fl)/NperLb)
hold all
get(gca,'ColorOrder')
end

figure(1)
subplot(4,1,1)
xlabel('Wave Period (s)');
ylabel('Motion Amplitude (m)');
legend(leg_str,'Location','NorthWest');
title(['Configuration ' num2str(config)  '    Effective mass = ' num2str(m) ' kg']);
subplot(4,1,2)
xlabel('Wave Period (s)');
ylabel('Peak Heave Acceleration (g''s)');
legend(leg_str,'Location','NorthEast');
if(config == 2)
    ylim([0 .75])
end
subplot(4,1,3)
xlabel('Wave Period (s)');
ylabel('Upper Spring Force Amplitude (lbs)');
legend(leg_str,'Location','NorthEast');
title(['Upper Spring Constant = ' num2str(Ku/NperLb,3) ' lbs/m']);
if(config == 2)
    ylim([0 2000])
end
subplot(4,1,4)
xlabel('Wave Period (s)');
ylabel('Lower Spring Force Amplitude (lbs)');
title(['Lower Spring Constant = ' num2str(Kl/NperLb,3) ' lbs/m']);
legend(leg_str,'Location','NorthWest');

print('-dpdf',['LinearlyDampedResults_Configuration' num2str(config)])



QuadResults = load(['QuadDampedResults_Config' num2str(config) '.mat']);

figure(1)
subplot(4,1,1)
plot(QuadResults.Periods,QuadResults.MotionAmp(1,:),'bx');
plot(QuadResults.Periods,QuadResults.MotionAmp(2,:),'gx');
plot(QuadResults.Periods,QuadResults.MotionAmp(3,:),'rx');
plot(QuadResults.Periods,QuadResults.MotionAmp(4,:),'cx');
plot(QuadResults.Periods,QuadResults.MotionAmp(5,:),'mx');

subplot(4,1,2)
plot(QuadResults.Periods,(1/g)*QuadResults.MotionAmp(1,:).*((2*pi./QuadResults.Periods).^2),'bx');
plot(QuadResults.Periods,(1/g)*QuadResults.MotionAmp(2,:).*((2*pi./QuadResults.Periods).^2),'gx');
plot(QuadResults.Periods,(1/g)*QuadResults.MotionAmp(3,:).*((2*pi./QuadResults.Periods).^2),'rx');
plot(QuadResults.Periods,(1/g)*QuadResults.MotionAmp(4,:).*((2*pi./QuadResults.Periods).^2),'cx');
plot(QuadResults.Periods,(1/g)*QuadResults.MotionAmp(5,:).*((2*pi./QuadResults.Periods).^2),'mx');


subplot(4,1,3)
plot(QuadResults.Periods,QuadResults.UpperForceAmp(1,:),'bx');
plot(QuadResults.Periods,QuadResults.UpperForceAmp(2,:),'gx');
plot(QuadResults.Periods,QuadResults.UpperForceAmp(3,:),'rx');
plot(QuadResults.Periods,QuadResults.UpperForceAmp(4,:),'cx');
plot(QuadResults.Periods,QuadResults.UpperForceAmp(5,:),'mx');
%ylim([0 500]);

subplot(4,1,4)
plot(QuadResults.Periods,QuadResults.LowerForceAmp(1,:),'bx');
plot(QuadResults.Periods,QuadResults.LowerForceAmp(2,:),'gx');
plot(QuadResults.Periods,QuadResults.LowerForceAmp(3,:),'rx');
plot(QuadResults.Periods,QuadResults.LowerForceAmp(4,:),'cx');
plot(QuadResults.Periods,QuadResults.LowerForceAmp(5,:),'mx');


print('-dpdf',['QuadAndLinearlyDampedResults_Configuration' num2str(config)])
%print('-dpng',['QuadAndLinearlyDampedResults_Configuration' num2str(config)])


