%One-way Wave equation Approximation Methods
%18.086 - Gilbert Strang
%written by: Joseph Sikora III
%February 16, 2006

clear all
close all
clc
set(0,'defaulttextfontsize',14)
set(0,'defaultaxesfontsize',14)

%reduce dx to improve result appearance (dx=1 or 0.5)
dx=4;
%dx=10;

%movie shutter speed (sec)
dT=0.005;
%movie [on=1]/[off=0]
movie_on=1;

c=1500; %sound speed

%Courant number - chooses dt such that r<=1;     % ************************
dt = dx / c * (1-2/100);                        % TRY ALSO WITH  *********
%dt = dx / c * (1+2/100);                        % (1+2/100) and  *********
%dt = dx / c * (1);                              % (1) !!!  ***************
dt_exact = dx / c * (1);                         % ************************

%axes --------
x  = -2000 : dx : 0;
t  =    0  : dt : 1;
t2 =    0  : dt_exact : 1;

%initialize u(x,t) step function
u=zeros(length(x),length(t));
u(:,1)=[zeros(round(length(x)*3/4),1) ;ones(round(length(x)/4),1)];
for count=1:length(t2)
    u(1:end-1,count+1)=(1-c*dt_exact/dx)*u(1:end-1,count)+c*dt_exact/dx*u(2:end,count);
end

u_up=u;
u_LF=u;
u_LW=u;

%remove sequence duration (500 sec)
x=x+500;

figure
for count=1:length(t)
    
    %upwind ----------------------------------
    u_up(1:end-1,count+1)=(1-c*dt/dx)*u_up(1:end-1,count)+c*dt/dx*u_up(2:end,count);
    
    %Lax-Wendroff -----------------------------
    u_LW(2:end-1,count+1)=u_LW(2:end-1,count)+c*dt/dx/2*(u_LW(3:end,count)-u_LW(1:end-2,count))+...
        +c^2*dt^2/dx^2/2*(u_LW(3:end,count)-2*u_LW(2:end-1,count)+u_LW(1:end-2,count));
    
    %Lax-Friedrichs -----------------------------   
    u_LF(2:end-1,count+1)=1/2*(u_LF(3:end,count)+u_LF(1:end-2,count))+...
        c*dt/2/dx*(u_LF(3:end,count)-u_LF(1:end-2,count));
    if movie_on==1
        subplot(311),area(x,u_up(:,count)),ylabel('upwind'),set(gca,'xticklabel',[],...
            'ylim',[min([-.2 min(u_up)]) max([1.2 max(u_up)])])
        title(sprintf('t=%.2f sec, c=%d m/s',t(count),c))
        subplot(312),area(x,u_LW(:,count)),ylabel('LW'),set(gca,'xticklabel',[],...
            'ylim',[min([-.2 min(u_LW)]) max([1.2 max(u_LW)])])
        subplot(313),area(x,u_LF(:,count)),ylabel('LF'),set(gca,...
            'ylim',[min([-.2 min(u_LF)]) max([1.2 max(u_LF)])])
        xlabel('x [meters]')
%        pause(dT)
    end
end

%picks a reasonable value of t to look at
t_ind=round(median(find(t>=.75 & t<.85)));
t_ind2=round(median(find(t2>=.75 & t2<.85)));

%picks a reasonable region of x to look at
x_ind=find(u_up(:,t_ind)>0);
L=find(x>min(x(x_ind))-7.5*dx & x<min(x(x_ind))+75*dx^.25);

%plot results 
figure
plot(x(L),u(L,t_ind2),'g',x(L),u_up(L,t_ind),'k',...
    x(L),u_LW(L,t_ind),'--b',x(L),u_LF(L,t_ind),'-.r','linewidth',2)
%set(gca,'ylim',[min([-.1 min(u_up) min(u_LW) min(u_LF)]) max([1.2 max(u_up) max(u_LW) max(u_LF)])])
set(gca,'ylim',[-.1 1.2])
xlabel('x [meters]')
title(sprintf('t=%.2f sec, c=%d m/s',t(t_ind),c))
legend('u(x,t)','upwind','LW','LF',4)