%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Solution to inviscid Burger's Equation
%         Using the Spectral Transform Method
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;


t_stop = 0.20;    % end time for integration
nskip  = 10;      % number of integration steps between pause/printout 

% define the spatial domain

nx =  64;         % number of x-grid points, also 
                  % number of retained Fourier modes
xmin = 0;         % left edge of domain
xmax = 1;         % right edge of domain
dt = .005;         % time step

% *** Define and compute misc. numerical paramters

dx = (xmax-xmin)/(nx);
xx=linspace(xmin,xmax-dx,nx);  % points on x-axis
zeta=xx;                       % x points to be used as initial positions
im = sqrt(-1);   % avoid confusion about i=sqrt(-1)
max_step = round(t_stop/dt);

% *** Initialize the approximate and exact solutions

psi = sin(2*pi.*xx);
exact = psi;


% *** Fourier transform the initial field

fpsi = fft(psi);

% *** Extra stuff for the spectral transform method

%   pad Fourier coefficients with zeros to obtain an array of 
%   length 2*nx for use with the transform method
%   Note: it would only be necessary to pad to a length of
%   1.5*nx, but padding to 2*nx makes the coding simpler

pad_fac=2;   % factor to account for the change in the normalization
             % of ifft(fpsi_2) by 2*nx since the vector length is now
             % doubled

fpsi_2 = pad_fac*[fpsi(1:nx/2) fpsi(nx/2+1)/2 zeros(1,nx-1) ...
                   fpsi(nx/2+1)/2. fpsi(nx/2+2:nx)];

% *** create arrays containing the wavenumber of each mode
%     Maximum wavenumber is 2*pi/(2*dx) = pi/dx = pi*nx/(xmax-xmin)

kscale = 2*pi/(xmax-xmin);

%           wavenumbers on double resolution mesh
knum_2 = [0:nx/2 zeros(1,nx-1) -nx/2:-1]*kscale;

% *** Initialize the time integration arrays

%     initialize the spectral method
sp_new = fpsi_2; sp_mid = fpsi_2; sp_old = fpsi_2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Main time integration loop

for step = 1:max_step;

    time=step*dt;
%
%  ****** spectral approximation implemented via the transform method *** 
%
   sp_old = sp_mid;
   sp_mid = sp_new;

%    Compute the spatial derivative of the concentration field
   fderiv_sp = im*sp_mid.*knum_2;

%    Compute forcing on the fine-spaced physical grid

   forcing_2 = real(ifft(fderiv_sp)).*real(ifft(sp_mid));

%     Transform back to wavenumber space and zero the
%     coefficients for the waves not retained in the original
%     truncation

   sp_forcing_2 = fft(forcing_2);

   sp_forcing = [sp_forcing_2(1:nx/2) 2*sp_forcing_2(nx/2+1) zeros(1,nx), ...
              sp_forcing_2(3*nx/2+2:2*nx)];

%   step forward: use an initial forward step, then leapfrog

   if (step == 1)
     sp_new = sp_old - dt*sp_forcing;
   else
     sp_new = sp_old - 2*dt*sp_forcing;
   end

%  **************** exact solution **************************

   for i = 1:nx,
%        for each i point, 
%            solve 0 = xx - zeta -sin(2*pi*zeta)*time for zeta
%            use zeta from the previous time step as the first guess
     ztemp = fzero(@(x) xx(i)-x-sin(2*pi*x)*time , zeta(i));
     zeta(i)=ztemp;
     exact(i) = sin(2*pi*zeta(i));
   end % for i

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%             Plot solutions

  sp=[sp_new(1:nx/2) 2*sp_new(nx/2+1) sp_new(3*nx/2+2:2*nx)]/pad_fac;
  spec=real(ifft(sp));

  figure(1)
  plot(xx,exact,'b-');
  hold on
  plot(xx,spec,'r--');
  axis([xmin xmax -1.2 1.2])
  xlabel(strcat('M = ',num2str(nx)));
  title(strcat('Time = ',num2str(time)));
  legend('exact','spectral');
  hold off
  pause(0.1)


  if ( (mod(step,nskip) == 0) || (step==max_step) )
    fprintf('At time %5.3f',time)
    fprintf('  sum Spectral**2 = %7.3f',sum(spec.*spec) )
    fprintf('\n Press any key...\n');
    pause
  end 

end % main integration step


