%
%   SOV_EX1.M     Heat Transfer 1-D Finite Bar (HT Example #1)
%
%   Analytical Solution using Separation of Variables (SOV) to the following 
%   problem:
%      ut(x,t) = alf*uxx(x,t)   with  u(0,t) = 0   ux(L,t) = 0   u(x,0) = f(x)
%      where    u(x,0) = f(x) = 100 C  (constant initial temperature)  
%
%   This problem models a 1-D laterally insulated bar which is initially at
%   a constant temperature of 100 C.  At time t = 0, the left end is placed
%   in an ice bath (0 C) and the right end is perfectly insulated (u' = 0).
%   The object of the simulation is to compute and visualize the space-time
%   temperature profile within the bar.
%
%   File prepared by J. R. White, UMass-Lowell (Nov. 2004)
%

%
%   getting started
      clear all,   close all,   nfig = 0;
%
%   problem data
      alf = 1.15;             % cm^2/s thermal diffusivity for copper
      L = 50;                 % cm     length of bar
      uo = 100;               % C      initial temp throughout bar
      Nx = 101;               % number of x values
      x = linspace(0,L,Nx)';  % vector of points to evaluate function
      nmax = 50;              % max number of nonzero terms
      tol = 0.001;            % tolerance to stop series evaluation
%   
%   calc ln = (2n-1)*pi/2L  and   an = 4uo/(2n-1)*pi)  for n = 1,2,3,...
      ln = zeros(1,nmax);   an =  zeros(1,nmax);
      d = pi/(2*L);    c = 4*uo/pi;    
      for n = 1:nmax
        ln(n) = (2*n-1)*d;    an(n) = c/(2*n-1);
      end
%
%   now evaluate series expansion for several different times
      tt = [0 30 120 300 600 1800];   Nt = length(tt);   ut = zeros(Nx,Nt);   
%
      for i = 1:Nt
        t = tt(i);   cc = -alf*t;   mrerr = 1.0;   n = 0;   u = zeros(size(x));
        while mrerr > tol   &   n < nmax
          n = n+1;  un = an(n)*exp(cc*ln(n)*ln(n))*sin(ln(n)*x);
          u = u + un;   rerr = un(2:Nx-1)./u(2:Nx-1);   mrerr = max(abs(rerr));
        end
        ut(:,i) = u;
        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-';'g:';'b-';'m:';'c-';'k:'];
%
%   plot curves of u(x,t) for various times
      nfig = nfig+1;  figure(nfig)
      for j = 1:Nt
        h(j) = plot(x,ut(:,j),scm(j,:),'LineWidth',2); hold on 
        st(j,:) = sprintf('%5.1f min',tt(j)/60);
      end 
      plot(x,uo*ones(size(x)),'r--','LineWidth',2)    % actual initial profile
      title('SOV Example #1 --  Bar Temperature Profile at Various Times')
      grid,xlabel('Distance (cm)'),ylabel('Temperature (C)')
      legend(h,char(st))
      hold off
%
%   end of example
