%  
%   HT1D_SOV2.M   Heat Transfer in a 1-D Finite Bar  (Example 10.2 in Class Notes)
%  
%   Analytical Solution using Separation of Variables (SOV) to the following problem: 
%      ut(x,t) = alf*uxx(x,t)   with  u(0,t) = ul   u(L,t) = ur   u(x,0) = u0 
%   has solution given by
%      u(x,t) = v(x,t) + w(x)     
%   where v(x,t) is the transient solution and w(x) is the steady state solution                          
%
%   The goal here is to evaluate and visualize the temperature profile in the 
%   1-D finite bar of length L.  Two options are provided:
%     1. Create snapshots of u(x) for several times.  This is an easy
%        way to visualize the full space-time solution, u(x,t).
%     2. Use Matlab's animation capability to create a movie that
%        shows the spatial temperature profile changing in time.
%
%   File prepared by J. R. White, UMass-Lowell  (Aug. 2003)
%

%
%   getting started
      clear all,   close all,   nfig = 0;
%  
%   problem data                                                         
      alf = 0.0001;    % m^2/s  thermal diffusivity                            
      ul = 100;        % C      fixed temp at left end                         
      ur = 400;        % C      fixed temp at right end                        
      u0 = 70;         % C      initial uniform temp of bar                    
      L = 1;           % m      length of bar                                 
      maxt = 50;       % max number of terms in expansion
      tol = 0.001;     % tolerance used to stop series expansion
%
%   the steady state solution w(x)
      Nx = 101;   x = linspace(0,L,Nx)';   w = (ur-ul)*x/L + ul;                                                   
%
%   now v(x,t) is given as an infinite series                                 
%   calc terms for n = 1,2,3,4,5,...,max                                   
      cc1 = 600*L/pi;   cc2 = 60/pi;  
      for n = 1:maxt
        lam(n) = n*pi/L;   c(n) = (1/n)*(cc1*(-1)^n + cc2*((-1)^n-1));
      end      
%
%   define time points and initialize space-time distribution
      opt = menu('Output plot option?','Create profiles at 5 time points', ...
                 'Make a movie with Nt frames');
      switch opt
        case 1     % use for plotting snapshots
          Nt = 5;   tt = [0 50 250 1000 5000];
        case 2     % use for making smooth movie
          Nt = input(' Input number of frames for movie (5-100): ');
          if Nt < 5,    Nt = 5;    end
          if Nt > 100,  Nt = 100;  end
          tt = linspace(0,5000,Nt); 
      end
%
%   start loop over time points
      u = zeros(Nx,Nt);
      for i = 1:Nt                                                     
        t = tt(i);   cc = -alf*t;   mrerr = 1.0;   n = 0;   v = zeros(size(x));
          while mrerr > tol   &   n < maxt
            n = n+1;  vn = c(n)*sin(lam(n)*x).*exp(cc*lam(n)*lam(n));  v = v + vn; 
            j = find(v);                    % finds indices of nonzero values of v(x)
            mrerr = max(abs(vn(j)./v(j)));  % compute max relative error
          end
        u(:,i) = v + w;
        disp([' Needed ',num2str(n),' terms for convergence at t = ',num2str(t),' s'])
      end
%
%   plot curves of u for various times (if opt = 1)
      switch opt
        case 1
          nfig = nfig+1;   figure(nfig)                        
          plot(x,u,'LineWidth',2)
          axis([0 1 0 400]);                        
          title('HT1D\_SOV2:  Bar Temperature for Several Times (Example 10.2)')                  
          grid,xlabel('Distance (m)'),ylabel('Temperature (C)')       
          for i = 1:Nt,   gtext(['t = ',num2str(tt(i)),' s']),   end
%
%   let's make a movie of the time-dependent profile (if opt = 2)
        case 2
          nfig = nfig+1;     figure(nfig) 
          plot(x,u(:,1),'LineWidth',2);    axis([0 L 0 400]); 
          M = moviein(Nt); 
          for i = 1:Nt
            plot(x,u(:,i),'LineWidth',2);  axis([0 L 0 400]);
            M(:,i) = getframe;
          end
%
%   play the movie one more time and show initial and final profiles
          movie(M,0);   hold on;
          plot(x,u(:,1),'LineWidth',2);    grid on;  
          title('HT1D\_SOV2:  Bar Temperature for Initial and Final Times (Example 10.2)');                  
          xlabel('Distance (m)'),ylabel('Temperature (C)');      
          gtext('initial'),gtext('final')
          hold off
      end
%
%   end of problem


