
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% Dynamic coupling between shallow-water sloshing and horizontal vehicle
% motion: Linear 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_exa(1:n+1)=0.0;
q_num(1:n+1)=0.0; 
t(1:n+1)=0.0;
x_num(1:NN+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 for comparison with the exact solutions:

nu=110.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));

C=0.03;

qhat=-(C/(alpha*h0))*cot(s);
Hhat(1:NN+1)=C*(cos(alpha.*a(1:NN+1))-(cot(s)*sin(alpha.*a(1:NN+1))));
Xhat(1:NN+1)=(C/(alpha*h0)).*((cot(s)*(1.0-cos(alpha.*a(1:NN+1))))...
             -sin(alpha.*a(1:NN+1)));

wn(1:NN+1)=(omega.*Xhat(1:NN+1))+(omega*qhat);

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

pn=(mv*qhat*omega)+((mf/L)*sigmaN);

q_exa(1)=0.0;
x_exa(1:NN+1)=0.0;
h_exa(1:NN+1)=h0;

Xn=x_exa;
qn=q_exa(1);
q_num(1)=qn;
h_num=h_exa;

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

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

% R=mv/mf;
% G=(pi^2)*R;
% nu=(4*g*h0*mf*G)/(L^2);
% q_num(1)=0.15;
% qn=0.15;
% qdot=0.0;
% 
% Xn(1:NN+1)=0.0;
% h_num(1:NN+1)=h0;
% wn(1:NN+1)=qdot;
% 
% sigma=da*sum(wn(2:NN));
% sigmaN=sigma+(0.5*da*(wn(1)+wn(NN+1)));
% 
% pn=(mv+mf)*qdot;

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


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

sum_1=0.50*da*(mf/L)*sum(wn(2:NN).^2);
sum_2=0.5*(mf/L)*g*h0*(1.0/da)*sum((Xn(2:NN+1)-Xn(1:NN)).^2);

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

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

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

gamma=1.0+(da/(L*R));

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

for i=1:n

     t(i+1)=i*dt;
     
     % exact solutions:
     q_exa(i+1)=qhat*sin(omega*t(i+1));
     h_exa(1:NN+1)=h0+(Hhat(1:NN+1)*sin(omega*t(i+1)));
     x_exa(1:NN+1)=a(1:NN+1)+(Xhat(1:NN+1)*sin(omega*t(i+1)));
     %
     
     % St\"ormer-Verlet (SV) for the coupled linear equations:
     phalf=pn-((dt/2.0)*nu*qn);
    
     whalf(2:NN)=wn(2:NN)+((dt/2.0)*g*h0*((Xn(3:NN+1)-2.0*Xn(2:NN)...
                 +Xn(1:NN-1))/(da^2)));
     
     sigma=da*sum(whalf(2:NN));
     sigmaNhalf=(1.0/gamma)*(sigma+((da/mv)*phalf));
     
     whalf(1)=(phalf/mv)-(sigmaNhalf/(L*R));
     whalf(NN+1)=whalf(1);
    
     Xsec(1)=0.0;
     Xsec(2:NN)=Xn(2:NN)+(dt*(whalf(2:NN)-(phalf/mv)+(sigmaNhalf/(L*R))));
     Xsec(NN+1)=0.0;
     
     qsec=qn+(dt*((phalf/mv)-(sigmaNhalf/(L*R))));
     q_num(i+1)=qsec;
     
     psec=phalf-((dt/2.0)*nu*qsec);
     
     wsec(2:NN)=whalf(2:NN)+((dt/2.0)*g*h0*((Xsec(3:NN+1)...
                -2.0*Xsec(2:NN)+Xsec(1:NN-1))/(da^2)));
     
     sigma=da*sum(wsec(2:NN));
     sigmaNsec=(1.0/gamma)*(sigma+((da/mv)*psec));

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

     x_num(1:NN+1)=a(1:NN+1)+Xsec(1:NN+1);
     h_num(1)=h0*(1.0-((Xsec(2)-Xsec(1))/da));
     h_num(2:NN)=h0*(1.0-((Xsec(3:NN+1)-Xsec(1:NN-1))/(2.0*da)));
     h_num(NN+1)=h0*(1.0-((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*da*(mf/L)*sum(wn(2:NN).^2);
sum_2=0.5*(mf/L)*g*h0*(1.0/da)*sum((Xn(2:NN+1)-Xn(1:NN)).^2);

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

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


if (i==2000)
       ex_1=x_num;
       eh_1=h_num;
end

if (i==7000)
       exx_1=x_num;
       ehh_1=h_num;
end

if (i==18000)
       exxx_1=x_num;
       ehhh_1=h_num;
end

if (i==26000)
       exxxx_1=x_num;
       ehhhh_1=h_num;
end

if (i==33000)
       ix_1=x_num;
       ih_1=h_num;
end

if (i==39000)
       ixx_1=x_num;
       ihh_1=h_num;
end


if (i==44000)
       ixxx_1=x_num;
       ihhh_1=h_num;
end

if (i==50000)
       ixxxx_1=x_num;
       ihhhh_1=h_num;
end

i

end

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

% figure;plot(t,q_num)

figure;plot(t,q_num,t,q_exa,'--')

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

