close all
clear all

config = 3;

global A T Ku Kl m b Cd rho Ac;

NperLb = 4.45;

Cd = 0.8;
rho = 1025;
d = 24*.0254;
Ac = 0.25*pi*d^2;

Amplitudes = 1:5;
Periods = 3:20;

for i = 1:length(Amplitudes)
A = Amplitudes(i);
leg_str{i} = ['A = ' num2str(A) 'm'];
for j = 1:length(Periods)

T = Periods(j);

w = 2*pi/T;

if(config == 1)
    m = 590;  %kg
    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;  %kg
    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;  %kg
    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

b = Cd*0.5*rho*Ac*A*w/sqrt(2);



t=0:0.001:10*T;   % time scale

x = A*cos(2*pi*t/T);

% SOLVE  d2y/dt2+5 dy/dt - 4 y = sin(10 t)
% initial conditions: y(0) = 0, y'(0)=0


initial_y    = 0;
initial_dydt = 0;

[t,y]=ode45( @rhs, t, [initial_y initial_dydt] );


Fu = Ku*(x'-y(:,1));
Fl = Kl*y(:,1);

idx = round(2*length(t)/3):length(t);

MotionAmp(i,j) = 0.5*(max(y(idx,1))-min(y(idx,1)));
UpperForceAmp(i,j) = 0.5*(max(Fu(idx)/NperLb)-min(Fu(idx)/NperLb));
LowerForceAmp(i,j) =0.5*(max(Fl(idx)/NperLb)-min(Fl(idx)/NperLb));

figure('Position',[1000 246 1028 1252]);
subplot(3,1,1)
plot(t,y(:,1));
xlabel('Wave Period (s)');
ylabel('Motion  (m)');
title(['A = ' num2str(A) ' m   T = ' num2str(T) ' s   Effective mass = ' num2str(m) ' kg']);
legend(['Amp = ' num2str(MotionAmp(i,j)) ' m'],'Location','NorthWest');
subplot(3,1,2)
plot(t,Fu/NperLb);
ylim([-1000 1000]);
xlabel('Wave Period (s)');
ylabel('Upper Spring Force  (lbs)');
title(['Upper Spring Constant = ' num2str(Ku/NperLb,3) ' lbs/m']);
legend(['Amp = ' num2str(UpperForceAmp(i,j)) ' lbs'],'Location','NorthWest');
subplot(3,1,3)
plot(t,Fl/NperLb);
xlabel('Wave Period (s)');
ylabel('Lower Spring Force  (lbs)');
title(['Lower Spring Constant = ' num2str(Kl/NperLb,3) ' lbs/m']);
legend(['Amp = ' num2str(LowerForceAmp(i,j)) ' lbs'],'Location','NorthWest');


end
end

figure('Position',[1000 246 1028 1252],'PaperType','usletter','PaperOrientation','portrait','PaperPositionMode','auto','PaperPosition',[.5 .5 7.5 10.5]);
ah1 = subplot(3,1,1)
plot(Periods,MotionAmp)
xlabel('Wave Period (s)');
ylabel('Motion Amplitude (m)');
title(['NonLinear Analysis, Cd = ' num2str(Cd) '   A = ' num2str(A) ' m     Effective mass = ' num2str(m) ' kg']);
legend(leg_str,'Location','NorthWest');
ah2 = subplot(3,1,2)
plot(Periods,UpperForceAmp);
xlabel('Wave Period (s)');
ylabel('Upper Force Amplitude (lbs)');
title(['Upper Spring Constant = ' num2str(Ku/NperLb,3) ' lbs/m']);
ah3 = subplot(3,1,3)
plot(Periods,LowerForceAmp);
xlabel('Wave Period (s)');
ylabel('Lower Force Amplitude (lbs)');
title(['Lower Spring Constant = ' num2str(Kl/NperLb,3) ' lbs/m']);

linkaxes([ah1 ah2 ah3],'x');

save(['QuadDampedResults_Config' num2str(config) '.mat'], 'Amplitudes','Periods', 'MotionAmp','UpperForceAmp','LowerForceAmp');

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



