%                                                                               
%   HT1D_SOV3.M   Heat Transfer in a 1-D Finite Bar (Example 10.3 in Class Notes)  
%                                                                               
%   Analytical Solution using Separation of Variables (SOV) to the following problem:
%      ut(x,t) = alf*uxx(x,t)    
%   where
%      u(L,t) = g(t) = ur*exp(-mu*t)    --- decaying temperature on right side
%      u(0,t) = ul                      --- fixed temp on left   
%      u(x,0) = f(x) = (ur-ul)x/L + ul  --- linear profile  
%   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 where the RHS endpoint temperature is time
%   dependent.  The Eigenfunction Expansion Method is used to approximate the 
%   v(x,t) solution (see class notes).                          
%                
%   File generated 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                             
      L = 1;            % m      length of bar                                   
      ul = 100;         % C      left (fixed) temp of bar                        
      ur = 400;         % C      initial temp on right                           
      maxt =  80;       % max terms in expansion
      tol = 0.001;      % tolerance used to stop series expansion
      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    
%
%   vector of points to evaluate function 
      Nx = 101;   x = linspace(0,L,Nx)';
%
%   w(x) is the steady state solution
      w = ul*(L-x)/L;                         
%
%   calc ln, ano and cn
      for n = 1:maxt                                 
        ln(n) = n*pi/L;    ano(n) = (2*ur/(pi*n))*(-1)^(n+1);
        cn(n) = (2*alf*ur*ln(n)/L)*(-1)^(n+1);
      end             
%
%   define time points and initialize space-time distribution
      Nt = 5;    tt = [0 600 1200 2400 4800];   % time in seconds
      u = zeros(Nx,Nt);
%
%   start loop over time points
      for i = 1:Nt                                                     
        t = tt(i);   bb = -mu*t;  cc = -alf*t;  v = zeros(Nx,1);
        n = 0;   mrerr = 1.0;
          while mrerr > tol   &   n < maxt
            n = n+1;  k1 = exp(bb);  k2 = exp(cc*ln(n)*ln(n));  k3 = alf*ln(n)*ln(n);
            vn = (ano(n)*k2 + cn(n)*(k1-k2)/(k3-mu))*sin(ln(n)*x);
            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
%
%   set color and marker code for creating plots
      Ncm = 6;
      scm = ['r-';      % red solid
             'g:';      % green dotted
             'b-';      % blue solid
             'm:';      % magenta dotted
             'c-';      % cyan solid
             'y:'];     % yellow dotted
%
%   plot curves of u for various times
      st = zeros(Nt,9);
      if Nt < 7 
        nfig = nfig+1;   figure(nfig)                        
        for i = 1:Nt
          h(i) = plot(x,u(:,i),scm(i,:),'LineWidth',2); hold on 
          st(i,:) = sprintf('%5.1f min',tt(i)/60);
        end 
        hold off 
        axis([0 1 0 400]);                        
        title('HT1D\_SOV3:  Bar Temp Profile with g(t) = ur*exp(-mu*t) (Example 10.3)')  
        grid,xlabel('Distance (m)'),ylabel('Temperature (C)')                      
        legend(h,char(st))
        text(0.04,20,['mu = ',num2str(mu),' 1/s']) 
      end
%
%   end of problem
