

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% Shallow-water sloshing on the London Eye with rotation about the
% suspension point.
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

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


clear
clc

tic

format long


g=9.81;
L=0;
R=1.0;
m=100;
Dx=R/m;
B=0;
h0=0.16;
T=4.0*h0;

radius=1.5;

           deltaT=0.01;
           omegafre=pi/R*sqrt(g*h0);

           epsilon=0.20 ;
           omega_n=omegafre;        
                   
           t(1)=0;
           
           L3=0;
           d1=-R/2;
           d2=-T-L3;
           
           
for i=1:m+1
    X(i)=(i-1)*((R-L)/m);
end

EH(1:m+1)=h0;
Ufirst(1:m+1)=0.0;           

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

aa=[L+d1;B+d2];
bbb=[R+d1;B+d2];
cc=[L+d1;T+d2];
dd=[R+d1;T+d2];

ee=[-2*radius-R;-radius+1.3*d2];
ff=[R;-radius+1.3*d2];
gg=[-2*radius-R;radius+h0];
hh=[R;radius+h0];

thetac=linspace(0,2*pi,2000);
xc(1:length(thetac))=-radius.*(1-cos(thetac(1:end)));
yc(1:length(thetac))=radius.*sin(thetac(1:end));   


theta0=0.0;

lin1x=linspace(aa(1),bbb(1),101);
lin1y=linspace(aa(2),bbb(2),101);

lin2x=linspace(aa(1),cc(1),101);
lin2y=linspace(aa(2),cc(2),101);

lin4x=linspace(bbb(1),dd(1),101);
lin4y=linspace(bbb(2),dd(2),101);



lin5x=linspace(ee(1),ff(1),101);
lin5y=linspace(ee(2),ff(2),101);

lin6x=linspace(ee(1),gg(1),101);
lin6y=linspace(ee(2),gg(2),101);

lin7x=linspace(gg(1),hh(1),101);
lin7y=linspace(gg(2),hh(2),101);

lin8x=linspace(ff(1),hh(1),101);
lin8y=linspace(ff(2),hh(2),101);

xfixed=[-radius radius-radius*cos(theta0)];
yfixed=[0 radius*sin(theta0)];

lin9x=linspace(0,cc(1),101);
lin9y=linspace(0,cc(2),101);

lin10x=linspace(0,dd(1),101);
lin10y=linspace(0,dd(2),101);

lin11x=linspace(cc(1),dd(1),101);
lin11y=linspace(cc(2),dd(2),101);

lin1xbase=lin1x;
set(gcf,'Color',[1,1,1])
hold on
plot(xfixed,yfixed,'r','LineWidth',2.5)
plot(xc,yc,'k','LineWidth',2.5)
plot(lin1x,lin1y,'k','LineWidth',2.5)
plot(lin2x,lin2y,'k','LineWidth',2.5)
plot(lin4x,lin4y,'k','LineWidth',2.5)
plot(lin5x,lin5y,'w','LineWidth',2.5)
plot(lin6x,lin6y,'w','LineWidth',2.5)
plot(lin7x,lin7y,'w','LineWidth',2.5)
plot(lin8x,lin8y,'w','LineWidth',2.5)
plot(lin9x,lin9y,'k','LineWidth',2.5)
plot(lin10x,lin10y,'k','LineWidth',2.5)
plot(lin11x,lin11y,'k','LineWidth',2.5)
plot(lin1x,EH+d2,'b','LineWidth',2.5)
plot(0,0,'o','MarkerFaceColor','w','MarkerEdgeColor','k','MarkerSize',7)
plot(-radius,0,'o','MarkerFaceColor','b','MarkerEdgeColor','b','MarkerSize',8)
axis off
hold off
set(gca,'DataAspectRatio',[1 1 1])
xlim([-2*radius-R R]);ylim([-radius+1.2*d2 radius+h0])
Image(1)=getframe;

constant=0.50;
for jj=1:1000
    
    t(jj+1)=jj*deltaT;
    tT=t(jj+1);
    
    thetaCircle=constant*(tT)+theta0;
    
    theta=epsilon*sin(omega_n*tT);
    bigomega=epsilon*omega_n*cos(omega_n*tT);
    bigomegadot=-epsilon*(omega_n^2)*sin(omega_n*tT);
    q1=-radius*(cos(theta0)-cos(thetaCircle));
    qddot1=-radius*(constant^2)*cos(thetaCircle);
    q2=radius*(sin(thetaCircle)-sin(theta0));
    qddot2=-radius*(constant^2)*sin(thetaCircle);
    
    
    alpha=(g*cos(theta))+(bigomegadot.*(X+d1))-((bigomega^2)*d2)-(qddot1*sin(theta))+(qddot2*cos(theta));
    beta=(-g*sin(theta))+(bigomegadot*d2)+((bigomega^2).*(X+d1))-(qddot1*cos(theta))-(qddot2*sin(theta));
    
A=cell(m+1);
b=cell(m+1,1);
Uzero=Ufirst;
Hzero=EH;


    id=1;
    qpp=1;
    while qpp>1E-8
    
    alphatilde=alpha-((bigomega^2).*Hzero)+(2*bigomega.*Uzero);
    
        
A{1,1}=[1 0;0 1];
A{1,2}=(Hzero(1)*deltaT/Dx)*[0 1;0 0];
b{1,1}=[1 0;0 0]*[EH(1);Ufirst(1)];

for i=3:m+1
    A{1,i}=zeros(2);
end


for ii=2:m
    
    if ii>2
       for i=1:ii-2
           A{ii,i}=zeros(2);
       end
    end
    
    Atilde=(deltaT/(2*Dx))*[alphatilde(ii) (Uzero(ii)+(2*bigomega*Hzero(ii)));Uzero(ii) Hzero(ii)];
    
    A{ii,ii-1}=-Atilde;
    A{ii,ii}=[-bigomegadot*deltaT 1;1 0];
    A{ii,ii+1}=Atilde;
    
    for i=ii+2:m+1
        A{ii,i}=zeros(2);
    end
    b{ii,1}=([0 1;1 0]*[EH(ii);Ufirst(ii)])+(deltaT*beta(ii)*[1;0]);
end

for i=1:m-1
    A{m+1,i}=zeros(2);
end

A{m+1,m}=(-Hzero(m+1)*deltaT/Dx)*[0 1;0 0];
A{m+1,m+1}=[1 0;0 1];
b{m+1,1}=[1 0;0 0]*[EH(m+1);Ufirst(m+1)];


AA=cell2mat(A);
bb=cell2mat(b);

x=AA\bb;
hsec=x(1:2:end);
usec=x(2:2:end);
Usec=usec';
Hsec=hsec';

qpp=max(abs(Hsec-Hzero)+abs(Usec-Uzero));

id=id+1;

Uzero=Usec;
Hzero=Hsec;

idd(jj)=id;

    end % while
   
Ufirst=Usec;
EH=Hsec;


    Q=[cos(theta) -sin(theta);sin(theta) cos(theta)];
    
    ar=(Q*aa)+[q1;q2];
    br=(Q*bbb)+[q1;q2];
    cr=(Q*cc)+[q1;q2];
    dr=(Q*dd)+[q1;q2];
    
    for j=1:m+1
        EHp{1,j}=(Q*[lin1xbase(j);EH(j)+d2])+[q1;q2];
    end
    

lin1x=linspace(ar(1),br(1),101);
lin1y=linspace(ar(2),br(2),101);

lin2x=linspace(ar(1),cr(1),101);
lin2y=linspace(ar(2),cr(2),101);

lin4x=linspace(br(1),dr(1),101);
lin4y=linspace(br(2),dr(2),101);


lin5x=linspace(ee(1),ff(1),101);
lin5y=linspace(ee(2),ff(2),101);

lin6x=linspace(ee(1),gg(1),101);
lin6y=linspace(ee(2),gg(2),101);

lin7x=linspace(gg(1),hh(1),101);
lin7y=linspace(gg(2),hh(2),101);

lin8x=linspace(ff(1),hh(1),101);
lin8y=linspace(ff(2),hh(2),101);

xfixed=[-radius -radius*(cos(theta0)-cos(thetaCircle))];
yfixed=[0 radius*(sin(thetaCircle)-sin(theta0))];

lin9x=linspace(-radius*(cos(theta0)-cos(thetaCircle)),cr(1),101);
lin9y=linspace(radius*(sin(thetaCircle)-sin(theta0)),cr(2),101);

lin10x=linspace(-radius*(cos(theta0)-cos(thetaCircle)),dr(1),101);
lin10y=linspace(radius*(sin(thetaCircle)-sin(theta0)),dr(2),101);

lin11x=linspace(cr(1),dr(1),101);
lin11y=linspace(cr(2),dr(2),101);


EHPX=cell2mat(EHp);
set(gcf,'Color',[1,1,1])
plot(xfixed,yfixed,'r','LineWidth',2.5)
hold on
plot(xc,yc,'k','LineWidth',2.5)
plot(lin1x,lin1y,'k','LineWidth',2.5)
plot(lin2x,lin2y,'k','LineWidth',2.5)
plot(lin4x,lin4y,'k','LineWidth',2.5)
plot(lin5x,lin5y,'w','LineWidth',2.5)
plot(lin6x,lin6y,'w','LineWidth',2.5)
plot(lin7x,lin7y,'w','LineWidth',2.5)
plot(lin8x,lin8y,'w','LineWidth',2.5)
plot(lin9x,lin9y,'k','LineWidth',2.5)
plot(lin10x,lin10y,'k','LineWidth',2.5)
plot(lin11x,lin11y,'k','LineWidth',2.5)
plot(EHPX(1,:),EHPX(2,:),'b','LineWidth',2.5)
plot(xfixed(2),yfixed(2),'o','MarkerFaceColor','w','MarkerEdgeColor','k','MarkerSize',7)
plot(-radius,0,'o','MarkerFaceColor','b','MarkerEdgeColor','b','MarkerSize',8)
hold off
axis off
set(gca,'DataAspectRatio',[1 1 1])
xlim([-2*radius-R R]);ylim([-radius+1.2*d2 radius+h0])
Image(jj+1)=getframe;


jj
end % jj

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% movie(Image)
% map=colormap;
% mpgwrite(Image,map,'Sloshing.mpg');
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

toc
ElapsedTime=toc;
