%  
%   HT1D_SS3.M   Heat Transfer in a 1-D Finite Bar Using the STATE SPACE Method
%       (1-D Space-Time Heat Transfer Example #3 in Class Notes - Time Dependent BCs)
%  
%   State Space Numerical Solution to the following problem: 
%     ut(x,t) = alf*uxx(x,t)   with u(0,t) = uL   u(L,t) = g(t)   u(x,0) = f(x) 
%   where
%     u(0,t) = uL                      --- fixed temp on left   
%     u(L,t) = g(t) = uR*exp(-mu*t)    --- decaying temperature on right side
%     u(x,0) = f(x) = (uR-uL)x/L + uL  --- linear profile 
%
%   Uses function file HT1D_SS3A.M within ODE23
%
%   File prepared by J. R. White, UMass-Lowell (updated Dec. 2007)
%
                 
%
%   getting started
      clear all,  close all,  nfig = 0;
      global A b mu
%  
%   problem data                                                         
      alf = .0001;     % m^2/s  thermal diffusivity                            
      ul = 100;        % C      left (fixed) temp of bar                        
      ur = 400;        % C      initial temp on right                           
      L = 1;           % m      length of bar 
      imu = menu('Select decay rate of right side temp', ...
                 'Use slow decay rate (mu = 0.001 1/s)', ...
                 'Use rapid decay rate (mu = 0.005 1/s)');
      if imu == 1,   mu = .001;   end    
      if imu == 2,   mu = .005;   end    
%
%   set dx increments and time domain parameters
      Nx = 24;   dx = L/(Nx+1);   x = dx:dx:L-dx;  % spatial grid
      to = 0;    tf = 4800;                        % time interval (sec)
%
%   set edit parameters  and initial temperature distribution     
      Ne = 6;                   % number of edit intervals (including initial time)
      DT = (tf-to)/Ne;          % interval width
      t = linspace(to,tf,Ne);   % times of interest for spatial snapshots
      f = ((ur-ul)*x/L + ul)';  % initial temperature distribution
%
%   setup coeff matrix
      A = zeros(Nx,Nx);   b = zeros(Nx,1);    aa = alf/(dx*dx);
      b(1) = aa*ul;   A(1,1) = -2*aa;  A(1,2) = aa;       % left boundary 
      for i = 2:Nx-1
        A(i,i-1) = aa;  A(i,i) = -2*aa;  A(i,i+1) = aa;   % interior nodes 
      end
      b(Nx) = aa*ur;  A(Nx,Nx-1) = aa;  A(Nx,Nx) = -2*aa; % right boundary
%
%   now use ode23 to solve system of ODEs -- gives spatial profile at each time point
      [tt,u] = ode23('ht1d_ss3a',t,f);
      u = u';  % store solution as u(x,t1), u(x,t2), ... where i^th column is for t_i
%
%   add boundary points to solution for ploting
      xp = [0 x L];   up = zeros(Nx+2,Ne);
      up(:,1) = [ul; u(:,1); ur];
      for j = 2:Ne;   up(:,j) = [ul; u(:,j); ur*exp(-mu*t(j))];   end
%
%   set color and marker code for creating plots
      Ncm = 6;  scm = ['r-';'g:';'b-';'m:';'c-';'k:']; 
%
%   plot curves of up for various times
      st = zeros(Ne,9);
      nfig = nfig+1;   figure(nfig);   clf reset;                        
      jj = 0;
      for j = 1:Ne
        jj = jj+1;  if jj > 6, jj = 1; end
        plot(xp,up(:,j),scm(jj,:),'LineWidth',2); hold on 
        st(j,:) = sprintf('%5.1f min',t(j)/60);
      end
      axis([0 L 0 400]); 
      hold off                          
      title('HT1D\_SS3:  Bar Temperature Profiles at Several Times (u_R = e^{-\mut})')                  
      grid,xlabel('Distance (m)'),ylabel('Temperature (C)')       
      legend(char(st),'Location','NorthWest')
      text(0.03,20,['\mu = ',num2str(mu),' s^{-1}']) 
%
%   end of problem
