%
%   SOV_EX2.M     Heat Transfer 1-D Finite Bar (HT Example #2)
%
%   Analytical Solutionn using Separation of Variables (SOV) to the following 
%   problem:
%           ut(x,t) = alf*uxx(x,t)   
%   with    u(0,t) = uL   and   -kux(L,t) = h(u(L,t)-uinf)   
%   where   u(x,0) = f(x) = uo = 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 a cold liquid bath at 10 C and the right end is exposed to a 
%   convective environment characterized by h/k and uinf.
%   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;
      global hok L ALPN uL uo CC1
%
%   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
      uL = 10;                % C      left boundary temperature
      uinf = 30;              % C      ambient temperature
      hok1 = [2.5e-4 2.5e-2]; % 1/cm   h/k for copper (with h = [10 1000] W/m^2-s)
      stl1 = [' LOW';'HIGH']; % titles to plots depending on choice of cooling rate
      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
%
%   choose high or low cooling rate at right surface
      ic = menu('Choose high or low cooling rate at right surface', ...
                   '         LOW Cooling Rate (h = 10 W/m^2-C)         ', ...
                   '        HIGH Cooling Rate (h = 1000 W/m^2-C)');
      hok = hok1(ic);   stl = stl1(ic,:);
%
%   compute steady state solution
      CC1 = hok*(uinf-uL)/(1+hok*L);   w = CC1*x + uL;
%   
%   calc eigenvalues
      ln = zeros(Nmax,1);
      for n = 1:Nmax
        g1 = (2*n-1)*pi/2;  g2 = (2*n+1)*pi/2;
        pn = fzero('sov_ex2a',[g1 g2]);
        ln(n) = pn/L;
      end
%
%   evaluate the norm of the eigenfunctions
      norm = zeros(Nmax,1);
      for n = 1:Nmax
        ALPN = ln(n);
        norm(n) = quadl('sov_ex2b',0,L);
      end
%
%   find expansion coeffs 
      an = zeros(Nmax,1);  
      for n = 1:Nmax
        ALPN = ln(n);
        an(n) = quadl('sov_ex2c',0,L)/norm(n);
      end
%
%   now evaluate transient solution expansion for several different times
      tt = [0 30 120 300 600 1800];   Nt = length(tt);   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 < Nmax
          n = n+1;  vn = an(n)*exp(cc*ln(n)*ln(n))*sin(ln(n)*x);
          v = v + vn;   rerr = vn(2:Nx-1)./v(2:Nx-1);   mrerr = max(abs(rerr));
        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-';'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,u(:,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 #2 --  Bar Temperature Profile at Various Times (', ...
             stl,' Cooling Rate)'])
      grid,xlabel('Distance (cm)'),ylabel('Temperature (C)')
      legend(h,char(st))
      hold off
%
%   end of example
