clear all
close all
%------------------------------
% Solution of the 1d sine-Gordon equation
% u_tt-u_xx+sin(u)=0
% Initial condition: u(0,x)=f(x), u_t(0,x)=g(x);
% x in (0,L) with
% no-flux boundary conditions: 
%-----------------------------
% Analytical kink solution:
% u(t,x)=4*arctan(exp((x-v*t)/sqrt(1-v^2)))
% v is a velocity of the kink, v<>1.
% Antikink solution:
% u(t,x)=4*arctan(exp(-(x-v*t)/sqrt(1-v^2)))
%-----------------------------
%-----------------------------
% The scheme: 
% u_ij+1=-u_i^j-1+2(1-a^2)u_i^j+a^2(u_i+1^j+u_i-1^j)
% a=dt/dx-Courant number;
% Stability condition: a<1, i.e., dt<dx.
%-----------------------------
L=20; % length of the string
M=400; % amount of points on the (0,L)
x=linspace(-L, L, M+1);
dx=x(2)-x(1);% space grid size
dt=0.05; % time grid size
a=dt/dx; % Courant-number
T=2000;%2000 % amount of time steps
tt=50;
rad=0.01;
ys=7;
v=0.2; % velocity of the kink
sh=10; % shift
om=v; % oscillation frequenz of the breather
b1=a^2;
b2=2*(1-a^2);
b3=b1;
%-----------------------------
u=zeros(T+1,M+1);
%Kink solution
 f=inline('4*atan(exp((y+sh)./sqrt(1-v.^2)))', 'y','v','sh'); % an initial wave u(x,)
 g=inline('-2*(v/sqrt(1-v.^2))*sech((y+sh)./sqrt(1-v.^2))', 'y','v','sh'); % an initial wave u(x,)

%Antikink solution
 %f=inline('4*atan(exp(-(y+sh)./sqrt(1-v.^2)))', 'y','v','sh'); % an initial wave u(x,)
 %g=inline('2*(v/sqrt(1-v.^2))*sech((y+sh)./sqrt(1-v.^2))', 'y','v','sh'); % an initial wave u(x,)

% Kink-antikink collision
% f=inline('4*atan(exp((y+sh)./sqrt(1-v.^2)))+4*atan(exp(-(y-sh)./sqrt(1-v.^2)))', 'y','v','sh'); % an initial wave u(x,)
% g=inline('-2*(v/sqrt(1-v.^2))*sech((y+sh)./sqrt(1-v.^2))-2*(v/sqrt(1-v.^2))*sech((y-sh)./sqrt(1-v.^2))', 'y','v','sh'); % an initial wave u(x,)

% Kink-kink collision
%f=inline('4*atan(exp((y+sh)./sqrt(1-v.^2)))+4*atan(exp((y-sh)./sqrt(1-v.^2)))', 'y','v','sh'); % an initial wave u(x,)
%g=inline('-2*(v/sqrt(1-v.^2))*sech((y+sh)./sqrt(1-v.^2))+2*(v/sqrt(1-v.^2))*sech((y-sh)./sqrt(1-v.^2))', 'y','v','sh'); % an initial wave u(x,)

%Breather
 %f=0; % an initial wave u(x,)
 %g=inline('4*sqrt(1-v^2)*sech(y*sqrt(1-v^2))', 'y','v','sh'); % an initial wave u(x,)


u(1,:)=f(x,v,sh);
%u(1,:)=0; %for breather

plot(x,u(1,:),'r','LineWidth',2)
%axis([-L-rad L+rad 0-10*rad 2*pi+10*rad])
axis([-L-rad L+rad -ys ys])
title(['t=' num2str(0)],'Color','b')
       xlabel('x')
       ylabel('u(x,t)')
       hold on 
       pause(0.5)
       hold off
%-----------------------------
% Calculate the matrix A
ah=b2*ones(M+1,1);
ao=b3*ones(M,1);
au=b1*ones(M,1);
A=diag(ah)+diag(au,-1)+diag(ao,1);
A(1,2)=2*b1;
A(M+1,M)=2*b3;
%-----------------------------
% Calculate the second time row
A2=A/2;
for i=1:M+1
    b(i)=sin(u(1,i));
    gg(i)=g(x(i),v,sh);
end;
um=dt*gg'+A2*u(1,:)'-(dt^2/2)*b'; % 
u(2,:)=um';
clear um b
%----------------------------- 
 for t=2:T
 for i=1:M+1
    b(i)=sin(u(t,i));
 end;
 um=-u(t-1,:)'+A*u(t,:)'-(dt^2)*b';
 u(t+1,:)=um';
if mod(t,tt)==0
       plot(x,u(t+1,:),'r','LineWidth',2);
       hold on
       %axis([-L-rad L+rad 0-10*rad 2*pi+10*rad])
       axis([-L-rad L+rad -ys ys])
       title(['t=' num2str(t)],'Color','b')
       xlabel('x')
       ylabel('u(x,t)')
       hold on 
       pause(0.5)
       hold off
    end;
end;
%----------------------------- 
tz=1:tt:T+1;
figure(2)
%mesh(x,tz,u(tz,:))
h=surf(x,tz,u(tz,:))
shading interp
%set(h,'EdgeColor','none','FaceColor','b')
light('Position',[-2,2,20])
lighting phong
material([0.4,0.6,0.5,30])
set(h,'EdgeColor','none','FaceColor',[0.2,0.2,0.9],'BackFaceLighting','lit')
alpha(0.8)
%view(-48,28)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
% tz=1:tt:T+1;
% figure(2)
% %mesh(x,tz,u(tz,:))
% surf(x,tz,u(tz,:))
% shading interp
% colormap jet
% view(-48,28)
% xlabel('x')
% ylabel('t')

