#!/usr/bin/octave -q

%-----------------------------------------------------------------------------
% Jonathan R. Senning <jonathan.senning@gordon.edu>
% Gordon College
% April 22, 1999
%
% Use FTCS (forward time, centered space) scheme to solve the heat equation
% in a thin rod.
%
% The equation solved is
%
%	du     d  du
%	-- = K -- --
%	dt     dx dx
%
% along with boundary conditions
%
%	u(xmin,t) = a(t)
%	u(xmax,t) = b(t)
%
% and initial conditions
%
%	u(x,tmin) = f(x)
%
%-----------------------------------------------------------------------------

% Determine the desired plot types.  A value of 1 means show the plot and
% a value of 0 means don't show the plot.

individual_curves = 1;		% Show curves during computation
color_image = 1;		% Show pseudocolor plot of solution
surface_mesh = 1;		% Show surface plot of solution
contour_plot = 1;		% Show contour plot of solution

% Set the value of the thermal diffusivity.

K = 0.25;

% Set size of image.  Also, these values are combined with interval sizes to
% compute h (spatial stepsize) and k (temporal stepsize).  Define array
% to hold image.

n = 16; %32;			% Number of spatial intervals
m = 512; %4096;			% Number of temporal intervals
%m = 2 * ( tmax - tmin ) / ( h * h );
u = zeros( n+1, m+1 );

% 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;

% 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 );
t = linspace( tmin, tmax, m+1 );

% Initial condition

u(:,1) = 100 * sin( pi * x )';

% Boundary conditions

u(1,:) = zeros( 1, m+1 );				% Left
u(n+1,:) = 60 * ( ( 1 .- cos( pi * t ) ) / 2.0 );	% Right

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

% Compute step sizes and the value of K * k / ( h * h ).
% Note that if k > ( h * h ) / 2 then this procedure may be unstable.

h = ( xmax - xmin ) / n;
k = ( tmax - tmin ) / m;
a = K * k / ( h * h );

if ( k > ( h * h ) / 2.0 )
    fprintf( stderr, 'Warning: may be unstable; k = %f > (h*h)/2 = %f\n', ...
    	     k, h * h / 2.0 );
end

% Find likely extremes for u

umin = min( min( u ) );
umax = max( max( u ) );

% 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( 'o-;t = %5.2f;', 0 );
    plot( x, u(:,1), fmt );
    axis( [xmin, xmax, umin, umax] );
%    xlabel( 'x' );
%    ylabel( 'Temperature' );
    drawnow();
    pause( 2 );
end

% Main loop.  This consists of updating a single row (timestep t) of the
% temperature function using values from the previous row (timestep t-k).
% We also plot the solution at each time step.

for j = 1 : m

    u(2:n,j+1) = a * u(1:n-1,j) + ( 1 - 2 * a ) * u(2:n,j) + a * u(3:n+1,j);

    if ( individual_curves != 0 )
        fmt = sprintf( 'o-;t = %5.2f;', j * k );
        plot( x, u(:,j+1), fmt );
        axis( [xmin, xmax, umin, umax] );
%        xlabel( 'x' );
%        ylabel( 'Temperature' );
        drawnow();
    end

end

if ( individual_curves != 0 )
    input( 'Press Enter to continue... ' );
end
   
% All done computing solution, now show the desired plots...
%
% Alter the colormap to get a nice one and then display the image.  In the
% image, x is the horizontal axis and t is the vertical axis - however, t
% increases in the downward direction.  Thus, the first line of the image
% corresponds to the initial condition.

if ( color_image != 0 || surface_mesh != 0 || contour_plot != 0 )
    cm = colormap( 'default' );
    cm(:,3)=(1.0.-[0:63]/63)';
    cm(:,2)=4*cm(:,1).*(1.-cm(:,1));
    colormap( cm );
end

if ( color_image != 0 )
    image( x, t, u' );
    xlabel( 'x' );
    ylabel( 't' );
    drawnow();
    input( 'Press Enter to continue... ' );
end

% Display a mesh (surface) plot showing the solution.

if ( surface_mesh != 0 )
    meshc( x, t, u' );
    view( 37.5, 45 );
    xlabel( 'x' );
    ylabel( 't' );
    drawnow();
    input( 'Press Enter to continue... ' );
end

% We could also create a contour graph with the following command.  Note that
% correct contour levels depend on the solution itself.

if ( contour_plot != 0 )
    umin = min( min( u ) );
    umax = max( max( u ) );
    levels = linspace( umin, umax, 21 );
    contour( x, t, u', levels );
    xlabel( 'x' );
    ylabel( 't' );
    drawnow();
    input( 'Press Enter to continue... ' );
end

% End of file