clear all
close all
%------------------------------
% Solution of the 1d wave equation
% u_tt-c^2*u_xx=0
% c-velocity;
% Initial condition: u(0,x)=f(x), u_t(0,x)=g(x);
% x in (0,L) with
% boundary conditions: u(t,0)=0, u(t,L)=0.
%-----------------------------
% Analytical solution:
% u(t,x)=(f(x-ct)+f(x+ct))/2+1/2c*int_x-ct^x+ct(g(y)dy)
%-----------------------------
%-----------------------------
% 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=c*dt/dx-Courant number;
%-----------------------------
%-----------------------------
L=1; % length of the string
c=1;
M=200; % amount of points on the (0,L)
x=linspace(0, L, M+1);
dx=x(2)-x(1);% space grid size
a=0.25;
dt=(a*dx)/c; % time grid size
T=1000; % amount of time steps
tt=10;
s=2; % shift
k=L; % steepness of the bump
rad=0.01;
ys=1.1;%23%1.1;
ob=1;
%-----------------------------
u=zeros(T+1,M+1);



f=inline('sin(ob*pi*y/k)', 'y','k','ob'); % an initial wave u(x,)
u(1,:)=f(x,k,ob);

 

u(1,1)=0; u(1,M+1)=0;

plot(x,u(1,:),'r','LineWidth',2)
hold on
plot(0, 0,'o','MarkerEdgeColor','r',...
                'MarkerFaceColor','r',...
                'MarkerSize',10) 
            hold on
            plot(L, 0,'o','MarkerEdgeColor','r',...
                'MarkerFaceColor','r',...
                'MarkerSize',10) 
axis([0-rad L+rad -ys ys])
       title(['t=' num2str(1)],'Color','b')
       xlabel('x')
       ylabel('u(x,t)')
       hold on 
       pause(0.5)
       hold off
   
%-----------------------------
% Calculate the second time row
u(2,1)=0; u(2,M+1)=0;
for i=2:M
    u(2,i)=(1-a^2)*u(1,i)+(a^2*(u(1,i+1)+u(1,i-1)))/2;
end;
%-----------------------------
%Calculate all other time rows
for t=2:T
    u(t+1,1)=0; u(t+1,M+1)=0;
    for i=2:M
        u(t+1,i)=-u(t-1,i)+2*(1-a^2)*u(t,i)+a^2*(u(t,i+1)+u(t,i-1));
    end;
    if mod(t,tt)==0
       plot(x,u(t+1,:),'r','LineWidth',2);
       hold on
       plot(0, 0,'o','MarkerEdgeColor','r',...
                'MarkerFaceColor','r',...
                'MarkerSize',10) 
            hold on
            plot(L, 0,'o','MarkerEdgeColor','r',...
                'MarkerFaceColor','r',...
                'MarkerSize',10) 
axis([0-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:T+1;
% figure(2)
% surfc(x,tz,u)
% shading interp
% colormap bone
%--------------------------
% figure(3)
% plot(x,u(1,:),x,u(T,:),'k','LineWidth',2)
% hold on
%        plot(0, 0,'o','MarkerEdgeColor','k',...
%                 'MarkerFaceColor','k',...
%                 'MarkerSize',10) 
%             hold on
%             plot(L, 0,'o','MarkerEdgeColor','k',...
%                 'MarkerFaceColor','k',...
%                 'MarkerSize',10)
%              hold on
%             plot(L/ob, 0,'o','MarkerEdgeColor','r',...
%                 'MarkerFaceColor','r',...
%                 'MarkerSize',5)
%             hold on
%             plot(2*L/ob, 0,'o','MarkerEdgeColor','r',...
%                 'MarkerFaceColor','r',...
%                 'MarkerSize',5)
%             hold on
%             plot(3*L/ob, 0,'o','MarkerEdgeColor','r',...
%                 'MarkerFaceColor','r',...
%                 'MarkerSize',5)
%             hold on
%             plot(4*L/ob, 0,'o','MarkerEdgeColor','r',...
%                 'MarkerFaceColor','r',...
%                 'MarkerSize',5)
%             hold on
%             plot(5*L/ob, 0,'o','MarkerEdgeColor','r',...
%                 'MarkerFaceColor','r',...
%                 'MarkerSize',5)
% axis([0-rad L+rad -ys ys])
