#!/usr/bin/octave -q

%-----------------------------------------------------------------------------
% Jonathan Senning <jonathan.senning@gordon.edu>
% Gordon College
% April 26, 1999
%
% Solve the one-dimensional wave equation
%
%	d  du       d  du
%	-- -- = c^2 -- --
%	dt dt       dx dx
%
% along with boundary conditions
%
%	u(xmin,t) = 0
%	u(xmax,t) = 0
%
% and initial conditions
%
%	u(x,tmin)  = f(x)
%	u'(x,tmin) = g(x)
%
%-----------------------------------------------------------------------------

individual_curves = 1;

% Set the value of the wave velocity

c = 1;

% Set size of the domain.  These values are combined with interval sizes to
% compute h (spatial stepsize) and k (temporal stepsize).

n = 256; 			% Number of spatial intervals
m = 600;			% Number of time steps

% We need three rows (timesteps) of data to be available at any given time.
% The finite difference stencil for this problem looks like a cross, where
% u(i,j+1) is computed using u(i-1,j), u(i,j), u(i+1,j) and u(i,j-1).  We
% compute the j+1 timestep information using the j and j-1 timestep
% information.  As we don't need to keep all previously computed data, only
% that for the most recent two timesteps, we can use a three-row matrix to
% hold the data being computed (j+1) and the data being used (j and j-1).  A
% separate index array, called z in this program, will be used to cycle
% through the rows of u so we don't actually need to copy them as the
% iteration progresses.  z(3) will always refer to the j+1 timestep
% information, z(2) refers to the j timestep and z(1) refers to the j-1
% timestep.

u = zeros( n+1, 3 );
z = 1:3;

% X position of left and right endpoints

xmin = 0;
xmax = 1;

% Interval of time: tmin should probably be left at zero

tmin = 0;
tmax = 2;

% Compute step sizes.

h = ( xmax - xmin ) / n;
k = ( tmax - tmin ) / m;
rho = ( c * c ) * ( k * k ) / ( h * h );

if ( rho >= 1 )
   fprintf( stderr, 'Warning: may be unstable rho = %f >= 1\n', rho );
end

% Generate x and t values.  These aren't really needed to solve the PDE,
% but they are useful for computing boundary/initial conditions and graphing.

x = linspace( xmin, xmax, n+1 );

% Initial conditions

function y = f( x , n )
    %y = zeros( 1, n+1 );
    %y(n/4+1:n/2+1) = 100 * ( 1 .- cos( 8 * pi * x(n/4+1:n/2+1) ) );
    x0 = 0.6;
    y = exp( -200 * ( x .- x0 ).^2 );
end

function y = g( x, n )
    y = 0 * x;
end

u(:,1) = f( x, n )';
u(:,2) = 0.5 * rho * ( f( x .- h, n )' + f( x .+ h, n )' ) ...
		+ ( 1 - rho ) * f( x, n )' + k * g( x, n )';

% Boundary conditions

u(1,:)   = 0.0;			% Left
u(n+1,:) = 0.0;			% Right

%-----------------------------------------------------------------------------
%	Should not need to make changes below this point :)
%-----------------------------------------------------------------------------

% Find likely extremes for u

umax = max( max( abs( u ) ) );
umin = -umax;

% Plot initial condition curve.  The "pause()" is used to allow time for the
% plot to appear on the screen before actually starting to solve the problem
% for t>0.

if ( individual_curves ~= 0 )
    fmt = sprintf( '-;t = %5.2f;', 0 );
    plot( x, u(:,1), fmt );
    axis( [xmin, xmax, umin, umax] );
    drawnow();
    pause( 2 );
end

% Main loop.

for j = 2 : m

    if ( individual_curves ~= 0 )
        fmt = sprintf( '-;t = %5.2f;', ( j - 1 ) * k );
        plot( x, u(:,z(2)), fmt );
        axis( [xmin, xmax, umin, umax] );
        drawnow();
    end

    u(2:n,z(3)) = rho * ( u(1:n-1,z(2)) + u(3:n+1,z(2)) ) ...
			+ 2.0 * ( 1 - rho ) * u(2:n,z(2)) - u(2:n,z(1));

    temp = z(1);
    z(1) = z(2);
    z(2) = z(3);
    z(3) = temp;
    
end

if ( individual_curves ~= 0 )
    fmt = sprintf( '-;t = %5.2f;', m * k );
    plot( x, u(:,z(2)), fmt );
    axis( [xmin, xmax, umin, umax] );
    drawnow();
end

input( 'Press Enter to quit... ' );

% End of file