
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% Dynamic coupling between shallow-water sloshing and horizontal vehicle
% motion: Nonlinear coupled problem. ----MATLAB version----
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% Written by H. Alemi Ardakani & T.J. Bridges,
% Department of Mathematics, University of Surrey, UK.


format long

clear
clc

tic

g=9.81;
rho=1000.0;
h0=0.08;
L=0.525;

mf=rho*h0*L;
mv=0.10*mf;

NN=5000;
n=50E+3;
dt=1E-4;
da=L/NN;

% % % % % % % % % % % Preallocating for speed:

q(1:n+1)=0.0; 
t(1:n+1)=0.0;
wn(1:NN+1)=0.0;
a(1:NN+1)=0.0;
H_N(1:n+1)=0.0;
Ev(1:n+1)=0.0;
whalf(1:NN+1)=0.0;
xsec(1:NN+1)=0.0;
wsec(1:NN+1)=0.0;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

for j=1:NN+1
     a(j)=(j-1)*da;
end

t(1)=0.0;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% % % % % % % Initial conditions [near resonance] with q(0) \neq 0 and/or
% qdot(0) \neq 0.

   % Linear problem frequencies:
% % nu=80.0;
% % G=(nu*(L^2))/(4.0*g*h0*mf);
% % R=mv/mf;
% % 
% % s=0.10;
% % ite=1.0;
% % while (ite>1E-14)
% %        sp=s-(((G/s)-(R*s)-tan(s))/((-G/(s^2))-R-(1.0+(tan(s)^2))));
% %        ite=abs(sp-s);
% %        s=sp;
% % end
% % 
% % alpha=2*s/L;
% % omega=alpha*sqrt(g*h0);
% % omega_f=(pi/L)*sqrt(g*h0);
% % omega_v=sqrt(nu/mv);
% % omega_fv=sqrt(nu/(mv+mf));
   %

   % R, G and nu near resonance:
% % R=mv/mf;
% % G=(pi^2)*R;
% % nu=(4*g*h0*mf*G)/(L^2);
   %

nu=80.0;
q(1)=0.0;
qn=q(1);
qdot=0.10;

xn(1:NN+1)=a(1:NN+1);
h(1:NN+1)=h0;
wn(1:NN+1)=qdot;
chi(1:NN+1)=h(1:NN+1);

sigma=rho*da*sum(wn(2:NN).*chi(2:NN));
sigmaN=sigma+(0.5*rho*da*((wn(1)*chi(1))+(wn(NN+1)*chi(NN+1))));

pn=(mv*qdot)+sigmaN;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% % % % % % % % % %  % energy at t=0.0 sec:

sum_1=0.50*rho*da*sum(chi(2:NN).*(wn(2:NN).^2));
sum_2=0.50*rho*g*(da^2)*sum((chi(1:NN).^2)./(xn(2:NN+1)-xn(1:NN)));

H_N(1)=sum_1+(0.25*rho*da*((chi(1)*(wn(1)^2))+(chi(NN+1)*(wn(NN+1)^2))))...
       +sum_2+(0.50*(1.0/mv)*((pn-sigmaN)^2))+(0.50*nu*(qn^2));

Ev(1)=(0.50*(1.0/mv)*((pn-sigmaN)^2))+(0.50*nu*(qn^2));

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

gamma=(rho*da/mv)*(0.50*(chi(1)+chi(NN+1)));

% % % % % % % % % % % %  time loop:

for i=1:n

     t(i+1)=i*dt;
     
     % St\"ormer-Verlet (SV) for the coupled nonlinear equations:
     phalf=pn-((dt/2.0)*nu*qn);
    
     whalf(2:NN)=wn(2:NN)+((0.25*g*dt*da*(1.0./chi(2:NN)))...
                 .*(((chi(1:NN-1).^2)./((xn(2:NN)-xn(1:NN-1)).^2))...
                 -((chi(2:NN).^2)./((xn(3:NN+1)-xn(2:NN)).^2))));
     
     sigmaNhalf=((1.0/(1.0+gamma))*rho*da*sum(chi(2:NN).*whalf(2:NN)))...
                +((gamma/(1.0+gamma))*phalf);
     
     whalf(1)=(1.0/mv)*(phalf-sigmaNhalf);
     whalf(NN+1)=whalf(1);
     
     xsec(1)=xn(1);
     xsec(2:NN)=xn(2:NN)+(dt*(whalf(2:NN)-(phalf/mv)+(sigmaNhalf/mv)));
     xsec(NN+1)=xn(NN+1);
     
     qsec=qn+(dt*((phalf/mv)-(sigmaNhalf/mv)));
     q(i+1)=qsec;
     
     psec=phalf-((dt/2.0)*nu*qsec);
     
     wsec(2:NN)=whalf(2:NN)+((0.25*g*dt*da*(1.0./chi(2:NN)))...
                 .*(((chi(1:NN-1).^2)./((xsec(2:NN)-xsec(1:NN-1)).^2))...
                 -((chi(2:NN).^2)./((xsec(3:NN+1)-xsec(2:NN)).^2))));
     
     sigmaNsec=((1.0/(1.0+gamma))*rho*da*sum(chi(2:NN).*wsec(2:NN)))...
               +((gamma/(1.0+gamma))*psec);

     wsec(1)=(1.0/mv)*(psec-sigmaNsec);
     wsec(NN+1)=wsec(1);
     %

     h(1)=chi(1)/((xsec(2)-xsec(1))/da);
     h(2:NN)=chi(2:NN)./((xsec(3:NN+1)-xsec(1:NN-1))/(2.0*da));
     h(NN+1)=chi(NN+1)/((xsec(NN+1)-xsec(NN))/da);
     
     xn=xsec;
     wn=wsec;
     qn=qsec;
     pn=psec;
     sigmaN=sigmaNsec;
   
% % % % % % % % % %  % energy of t(i+1)=i*dt:

sum_1=0.50*rho*da*sum(chi(2:NN).*(wn(2:NN).^2));
sum_2=0.50*rho*g*(da^2)*sum((chi(1:NN).^2)./(xn(2:NN+1)-xn(1:NN)));

H_N(i+1)=sum_1+(0.25*rho*da*((chi(1)*(wn(1)^2))+(chi(NN+1)*(wn(NN+1)^2))))...
       +sum_2+(0.50*(1.0/mv)*((pn-sigmaN)^2))+(0.50*nu*(qn^2));

Ev(i+1)=(0.50*(1.0/mv)*((pn-sigmaN)^2))+(0.50*nu*(qn^2));

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 


if (i==2000)
       ex_1=xn;
       eh_1=h;
end

if (i==7000)
       exx_1=xn;
       ehh_1=h;
end

if (i==18000)
       exxx_1=xn;
       ehhh_1=h;
end

if (i==26000)
       exxxx_1=xn;
       ehhhh_1=h;
end

if (i==33000)
       ix_1=xn;
       ih_1=h;
end

if (i==39000)
       ixx_1=xn;
       ihh_1=h;
end


if (i==44000)
       ixxx_1=xn;
       ihhh_1=h;
end

if (i==50000)
       ixxxx_1=xn;
       ihhhh_1=h;
end

i

end

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% figure;plot(t,q_num)

figure;plot(t,q,'k','LineWidth',2)

figure;subplot(2,1,1);plot(t(1:end),H_N(1:end)-H_N(1));
subplot(2,1,2);plot(t(1:end),Ev(1:end))

figure;subplot(4,2,1);plot(ex_1,eh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15]);text(0.1,0.12,'t=0.2 s')

subplot(4,2,2);plot(exx_1,ehh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15])

subplot(4,2,3);plot(exxx_1,ehhh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15])

subplot(4,2,4);plot(exxxx_1,ehhhh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15])

subplot(4,2,5);plot(ix_1,ih_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15])

subplot(4,2,6);hold on;plot(ixx_1,ihh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15])

subplot(4,2,7);plot(ixxx_1,ihhh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15])

subplot(4,2,8);plot(ixxxx_1,ihhhh_1,'k','LineWidth',2);
xlim([0 L]);ylim([0 0.15]);text(0.1,0.12,'t=5.0 s')

t_cpu=cputime;

toc
