%  This Matlab program is used in 
%  MATH 365 - Computational Fluid Dynamics at 
%  James Madison University. This program was
%  written by James Sochacki of the 
%  Department of Mathematics at James Madison 
%  University as part of the NSF Grant - 
%  A Collaborative Computational Science Program 
%  between JMU and North Carolina Central University. 
  
%  Jim Sochacki
%  Department of Mathematics
%  James Madison University
%  Harrisonburg, VA   22807
%  jim@math.jmu.edu

       disp(' This code looks at the leap frog finite difference solution ') 
       disp(' to the advection diffusion equation ')
       disp(' w_t + c w_x + d w_(xx) = 0 ; -L < x < L, 0 < t < T '  )
       disp(' w(x,0) = q arctan(s*x)+r ')

       disp(' Give c and d ')
       c = input(' c ');
       d = input(' d ');
       disp(' Give L ')
       L = input(' ');

       disp(' Give a_left and a_right ')
       a_left = input(' a_left ');
       a_right = input(' a_right ');
       disp(' Give s ')
       s = input(' ');
 
    dx = input(' Give the grid size. ');
    dt = input(' Give the time step. ');
    T = input(' Give the time T to run the model. ');

    ntstps = T/dt;

      disp(' Give an integer indicating how often a snapshot ')
      disp(' of the profile should be taken. ')
      interval=input(' ');
      
%
% Set up the parameters for the leap frog scheme.
%
       c2 = -2.*d*dt/(dx*dx);
       c1 = -c*dt/dx;
       jmax = (2*L)/dx+1;

%  set up the initial condition

       for j = 1:jmax

         w0(j) = 0;
         w1(j) = 0;
         w2(j) = 0;
         x(j) = -L+j*dx;

       end

       q = (a_right-a_left)/pi;  % q and r gaurantee w->a_left as x->-infinity
       r = (a_right+a_left)/2;   % and w->a_right as x->+infinity

       for j=1:jmax

         w0(j) = q*atan(s*x(j))+r;

       end

       for j=1:jmax

         w1(j) = w0(j);
         w2(j) = w0(j);

       end
            
        plot(x,w0);
        pause(0.25);

       
       if (c<0) then    

         for j = 2:jmax-1
 
           w1(j) = w0(j)+(c*dt/dx)*(w0(j+1)-w0(j))+c2*(w0(j+1)-2.*w0(j)+w0(j-1));

         end

       else    

         for j = 2:jmax-1
 
           w1(j) = w0(j)+(c*dt/dx)*(w0(j)-w0(j-1))+c2*(w0(j+1)-2.*w0(j)+w0(j-1));

         end

       end

       count = 0;

%  start the time iterations

       for n = 1:ntstps 

         for j=2:jmax-1

           w2(j) = w0(j)+c1*(w1(j+1)-w1(j-1))+c2*(w1(j+1)-2.*w1(j)+w1(j-1));

         end

         count = count+1;

%
% Create the snapshots.

          if (count == interval)

            plot(x,w2);
            pause(0.25);
            count = 0;
 
          end
%
%  Update the calculations in time
%
          for j=1:jmax
            w0(j)=w1(j);
            w1(j)=w2(j);
          end

      end
