%
%   PLANEWALL_1.M     Heat Transfer 1-D Planar Wall (Copper Plate)
%
%   The program implements the analytical solution using the Separation of Variables
%   (SOV) method to a relatively simple heat transfer problem.  The result of the
%   SOV method (see notes) is an infinite series expansion, and the goal here is 
%   to evaluate and plot the T(x,t) distribution given by the infinite series.
%
%   The physical problem models 1-D space-time heat transfer in a large planar wall.
%   The wall is initially at a constant temperature of 100 C.  At time t = 0, the 
%   left face is exposed to an ice bath (0 C) and the right surface is perfectly 
%   insulated (no heat transfer or T' = 0).  The object of the simulation is to 
%   compute and visualize the space-time temperature profile within the wall (in 
%   this case, a copper plate 10 cm thick).
%
%   This represents another example of evaluating infinite series...
%
%   File prepared by J. R. White, UMass-Lowell (last update: Oct. 2017)
%

      clear all,   close all,   nfig = 0;
%
%   problem data
      alf = 1.15;             % cm^2/s thermal diffusivity for copper
      L = 10;                 % cm     thickness of plate
      Ti = 100;               % C      initial temp throughout plate
      Nx = 101;               % number of x values
      x = linspace(0,L,Nx)';  % vector of points to evaluate function
      nmax = 50;              % max number of nonzero terms in expansion
      tol = 0.001;            % tolerance to stop series evaluation
%   
%   calc ln = (2n-1)*pi/2L  and   an = 4*Ti/(2n-1)*pi)  for n = 1,2,3,...
      ln = zeros(1,nmax);   an =  zeros(1,nmax);
      d = pi/(2*L);    c = 4*Ti/pi;    
      for n = 1:nmax
        ln(n) = (2*n-1)*d;    an(n) = c/(2*n-1);
      end
%
%   choose type of plot to be produced
      ic = menu('Choose type of plot', ...
                'Show T(x) for several different times (2-D plot)  ', ...
                'Show T(x,t) in a variety of surface plots (3-D plots) ');
%
%   now evaluate series expansion for several different times 
      if ic == 1
        tt = [0 1 5 15 30 60];        % use this for simple 2-D plot
      else
        tt = linspace(0,120,49);      % use this for 3-D plots (more time points)
      end
      Nt = length(tt);   T = zeros(Nx,Nt);   
%
%   since the initial temperature profile is known, let's just set it here
      T(:,1) = Ti*ones(size(x));
%
%   now loop over remaining times
      for j = 2:Nt
        t = tt(j);   cc = -alf*t;   mrerr = 1.0;   n = 0;   Tx = zeros(size(x));
        while mrerr > tol   &&   n < nmax
          n = n+1;  Txn = an(n)*exp(cc*ln(n)*ln(n))*sin(ln(n)*x);
          Tx = Tx + Txn;   
          i = find(Tx);  rerr = Txn(i)./Tx(i);   mrerr = max(abs(rerr));
        end
        T(:,j) = Tx;
        disp([' Needed ',num2str(n),' terms for convergence at t = ',num2str(t),' s'])
      end
%
%   *****  Create 2-D plot  *****
      if ic == 1
%
%   set color and marker code for creating 2-D plots
      Ncm = 6;  scm = ['r-';'g:';'b-';'m:';'c-';'k:']; 
      h = zeros(Nt,1);  st = zeros(Nt,9);
%
%   plot curves of T(x,t) for various times
      nfig = nfig+1;  figure(nfig)
      for j = 1:Nt
        h(j) = plot(x,T(:,j),scm(j,:),'LineWidth',2); hold on 
        st(j,:) = sprintf('%5.1f sec',tt(j));
      end 
      title('PlaneWall\_1:  Plate Temperature Profile at Various Times')
      grid,xlabel('Distance (cm)'),ylabel('Temperature (C)')
      legend(h,char(st),'Location','SouthEast')
%
      end
%
%   *****  Create various 3-D plots  *****
      if ic == 2
%
      nfig = nfig+1;  figure(nfig);  colormap(jet) 
      surf(tt,x,T), shading interp; colorbar, view(60,30), hold on
      [cc,hh] = contour3(tt,x,T, [80 60 40 20 10 5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis('tight')
      title('PlaneWall\_1:  T(x,t) in Copper Plate')
      ylabel('x location (cm)'),   xlabel('time (sec)')
      zlabel('Temperature (C)'),   hold off                                 
%
      nfig = nfig+1;  figure(nfig);  colormap(jet)
      surf(tt,x,T), shading interp; colorbar, view(2), hold on
      [cc,hh] = contour3(tt,x,T, [80 60 40 20 10 5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis('tight')
      title('PlaneWall\_1:  T(x,t) in Copper Plate')
      ylabel('x location (cm)'),   xlabel('time (sec)')
      zlabel('Temperature (C)'),  gtext('Temperature (C)'),  hold off 
%
      end
%
%   end of example
