%  Explore nature of transient diffusion response
%  (Patankar Ch 4, p. 54).
%   This is simplest possible problem, to demonstrate fundamentals
%   of transient diffusion with constant kappa. 
%
%  d phi/dt = kappa d^2 phi /dx^2                        (1)
%
%   Boundary Condition
%   Neigbors at +- dx are held at 
%       phi_E = 0 and phi_W = 0                          (2)
%   
%   Initial condition
%   At time t=0, phi_P = phi_0, and phi is piecewise linear
%    i.e. a "tent" between x_W and x_E
%
%       phi(x,0) = phi_0 ( 1 + x/dx )   -dx <= x < 0
%                = phi_0 ( 1 - x/dx )     0 <= x <= dx   (3) 
%                = 0                  x < -dx or x > dx 
%
%  This is a classic problem that can be solved by Fourier analysis
%  if we assume that the "tent" is periodic with period 4 dx, and
%  goes to -phi_0 at +- 2 dx.
%
%  This code explores impact of discretizing space and time by
%  comparing solutions after time (or a sinfgle time step) dt
%  in units of 
%      tau = dx^2/(2 kappa)                              (4)
%  for 3 cases
%     (1) continuous in x and t
%     (2) discretized space dx, continuous in time
%     (3) discretized both space and time
%
%
% (1) When space and time are continuous, solution is given by
%   Fourier series -
%  -----------------------------------------------------------
%   
%  phi(x,t) =  [ phi_0/(4 dx) ] * sum_{m=-inf}^inf {
%        8 dx/(m pi)^2 (1 - cos(m pi) ) cos(m pi x /(2 dx) )
%      * exp( -kappa ( (m pi)/(2 dx) )^2 t )       }
%
%   Because cosine is symmetric function, all the negative m's are
%   equivalent to the positive m's, and m=0 term is zero (no DC)
%   so we can replace the 
%     sum_{m=-inf}^inf  with  2 sum_{m=1}^inf
%   
%      phi(x,t) = phi_0 * sum_{m=1}^inf  {
%        4/(m pi)^2 (1 - cos(m pi) ) cos(m pi x /(2 dx) )
%      * exp( -kappa ( (m pi)/(2 dx) )^2 t )  }       (5)
%
%  At center x_P
%      phi_P(t) = phi_0 * sum_{m=-inf}^inf {
%              4/(m pi)^2 (1 - cos(m pi) ) 
%            * exp( -kappa ( (m pi)/(2 dx) )^2 t )  }    (6)
%
%  (2) With spatial discretization,
%  --------------------------------
%    d phi_P/dt = kappa ( phi_E - 2 phi_P + phi_W )/dx^2
%               = - (2 kappa/dx^2) phi_P                 (7)
%
%   because phi_E = phi_W = 0.  
%
%   Solution of (7) for continuous time is
%       phi_P(t) = phi_0 * exp( -2 kappa t / dx^2 )      (8)
%
%
%   (3) With both spatial and temporal discretization,
%   ---------------------------------------------------
%   ( phi_P - phi_P^0 ) = ( kappa dt / dx^2 ) * {
%             theta ( phi_E - 2 phi_P + phi_W )          (9)     
%           + (1 - theta) ( phi_E^0 - 2 phi_P^0 + phi_W^0 )  }
%
%  BC sets phi_E = phi_W = phi_E^0 = phi_w^0 = 0
%  so that for time step dt, (9) reduces to
%
%   phi_P(dt) = phi_P^0 [ ( 1 - (1 - theta)(dt/tau) ) /
%                     ( 1 +  theta (dt/tau) )    ]
%
%
%  -------------------------------------------------------------------
%
%  Set spatial node separation
    dx = 1;
%
%  set up detailed space vector for plotting exact solution
    nx = 101;
    x_vec = linspace( -dx, dx, nx );
%
%  diffusivity
    kappa = 1;
%
%  Value of phi at x_P at t=0
    phi_P_0 = 1;
%
%  Characteristic time  tau = 2 kappa / dx^2
    tau = dx^2 / (2 * kappa );
%
%
%  set time step dt in numerical solution 
%  in units of tau = dx^2/(2 kappa)
    dt = 0.05;
%
%  set up time vector tau_vec in units of tau
    tau_length = 4;
    n_times = fix( tau_length/dt )  + 1;
    tau_vec = 0: dt: tau_length;
%
%  set up t_vec in units of time
    t_vec = tau_vec * tau;
%
%
%
%  (1) Solution with continuous space and time
%  ============================================
    m_max = 50;   %  highest term in cosine series
%
%  make text string about m_max for plot
    m_max_words = 'Number of cosine terms m_{max} =   ';
    m_max_str    = num2str(m_max);
    m_max_message  = [ m_max_words  m_max_str ];
%
    fig_no = 1;
    figure( fig_no );
    hold on;grid on;box on;
    set(gca, 'ylim', [0, 1] );
    xlabel('Distance  x (in units of future grid dx)')
    ylabel('\phi_0(x)')
    title('Construction of initial profile \phi_0(x) from cosine sum')
    text( -0.5, 0.32, 'Blue - m even; Red - m odd');
    text( -0.5, 0.12, m_max_message );
%
    phi_0 = zeros(size(x_vec) );
%
% find current phi_0 profile m after m terms
%    Get plot color, and plot profile
    for m = 1:m_max   
               phi_0 = phi_0 + phi_P_0 ...
               * (  ( 4 /(m*pi)^2)  * (1 - cos(m*pi) ) ...
               * cos(m*pi*x_vec/(2*dx) ) );
       if( rem(m,2) == 0)
         color = 'b';
       else
         color='r';
       end     %  Get red or blue color ...
%
       plot( x_vec, phi_0, 'color', color );
%
       if(m==1)
         cosine_instruction_message =  ...
               'Hit any key to continue building cosine series';
         text( -0.5, 0.04, cosine_instruction_message );  
         disp( cosine_instruction_message );
   %
       end   %  if(m == 1 ) ...
   %
   %
   %  show current m
        current_m_message = ['Current m =  ' num2str(m) ];
        patch( [-0.5 -0.5 -0.05 -0.05 ], [0.2 0.26 0.26 0.2 ], 'w' )
        text( -0.5, 0.22, current_m_message )
   %
   %
       pause
   %
    end    %   for m = 1:m_max ...
%
%
%  now find solution at later times
%  --------------------------------
  %
  %  make strings to report dt and max time
     dt_words = 'Plot interval (in units of \tau) =  ';
     dt_string = num2str(dt);
     time_step_message = [ dt_words dt_string ];
  %
     max_time_words = 'Maximum time =  ';
     max_time_string = num2str( tau_length );
     max_time_message = [ max_time_words  max_time_string ];
  %
     time_step_instruction_message =  ...
         'Hit any key to show next time step';

  %
  % start figure to plot evolution of profiles phi_t(x,t)
    figure( fig_no + 1 );
    hold on;grid on;box on;
    plot( x_vec, phi_0, 'k' )
    xlabel('Distance  x  (in units of grid dx)')
    ylabel('\phi(x,t)')
    title('Evolution \phi(x,t)')
    text( 0.1, 0.94, time_step_message )
    text( 0.25, 0.83, max_time_message );
    text( -0.9, 0.84, time_step_instruction_message );  
%
    disp( time_step_instruction_message );
%
%
% loop over times
%  --------------
    for i_time = 1:n_times
  %
  % sum over attenuated cosine components at step i_time
     phi_x_t = zeros(size(phi_0) );
   %
     for m = 1:m_max
       phi_x_t = phi_x_t + phi_P_0 ...
            * ( 4 /(m*pi)^2) * ( 1 - cos(m*pi) ) ...
            * cos(m*pi*x_vec/(2*dx) ) ...
            * exp( -kappa * (m*pi/(2*dx) )^2 * t_vec(i_time) );
     end     %  for m = 1:m_max
  %
  % plot this time step
     plot(x_vec, phi_x_t, 'r' )
  %
  %
  % show current time
     current_time_message = ...
          ['Current time =  ' num2str( (i_time-1)*dt)];
     patch( [ 0.33 0.33 0.9 0.9 ], [0.7 0.76 0.76 0.7 ], 'w' )
     text( 0.35, 0.72, current_time_message )
   %
     pause
   %
    end      %  for i_time = 1:n_times
%
%  recalculate and plot temporal response at x_P
%  ---------------------------------------------
%
    figure( fig_no + 2 )
    hold on;grid on;box on;
    xlabel('Time step dt  (in units of \tau = dx^2/(2 \kappa)')
    ylabel('\phi_P(t)')
    title('Evolution \phi_P(t). k-exact; r-discrete x; others-discrete x &< t, various \theta')
%
%  initialize solution vector at x_P
    phi_cs_ct = zeros( size(t_vec) );
%
%  loop on times
    for i_time = 1:n_times
  %
  %  loop on wavenumber summation
      for m = 1:m_max
       phi_cs_ct(i_time) = phi_cs_ct(i_time) +  ...
              phi_P_0 * (4/(m*pi)^2) * (1 - cos(m*pi) ) ...
            * exp( -kappa * ( m*pi/(2*dx) )^2 * t_vec(i_time) );
  %
      end  %  for m = 1:m_max
% 
    end  %  for i_time - 1:n_times
%
%  Plot it as function of time in tau units
    plot( tau_vec, phi_cs_ct, 'k' )
% 
      pause
%
%
%  (2) Solution phi(t) spatial discrete space, continuous time
%  ============================================================
%
    phi_ds_ct = phi_P_0 * exp( - t_vec / tau );
%
%  Plot it
    plot( tau_vec, phi_ds_ct, 'r' )
%
    pause
%
%
%  (3) Solution with discrete time and space
%  ==========================================
%

%   initialize theta and curve-coun ter i_curve
    theta = 1;
    i_curve = 0;
%
%  Loop on values of time-splitting parameter theta
    while( ~isempty( theta ) )
  %
  %  augment curve counter
      i_curve = i_curve+1;
  %
  %   find current level to write theta & color on plot
      ylevel_1 = 1 - i_curve/10;
      ylevel_2 = ylevel_1 + 0.08;
  %
  %  Get time-step splitting ratio theta
      disp('Just hit RTN if you don''t want another plot.');
      theta = input('What value for time-splitting parameter theta?   ');
  %
      if( ~isempty( theta ) )
    %
    %  Get plotting color
        color = input('What color for this curve?  ');
        phi_ds_dt = phi_P_0 * ...
          ( ( 1 - (1-theta)* t_vec/tau ) ./ ...
            ( 1 + theta * t_vec/tau ) );
    %
    % Plot it
       plot( tau_vec, phi_ds_dt, 'color', color )
    %
    %  show theta and color for current curve
        current_curve_message = ['\theta = ' num2str(theta) ...
                 '   color '  color ];
        patch( [ 1  1 2 2 ], ...
             [ ylevel_1 ylevel_2 ylevel_2 ylevel_1 ], 'w' )
        text( 1.05, ylevel_1+0.04, current_curve_message )
   %

      end   %   if( ~isempty( theta ) ) ...
  %
    end  %  while( ~isempty( theta ) ...
%
%  
%
%