%
%  This code solves the nonlinear, dispersive, diffusive wave equation:
%
%		u  + c u + b u u  - g u  + d u    = 0
%		 t      x       x      xx     xxx
%
%  For c =  0.0;  b =  6.0;  g =  0.0; d =  1.0; i.e. the KdV equation:
%
%		u  + 6 u u  +  u    = 0
%		 t        x     xxx
%
%  using a spectral method over the interval [0:L]
%  Ricardo Carretero (30/10/02) [based on David Muraki's code]
%
%  Note: The initial condition HAS TO BE PERIODIC on the box [0:L]
%

global ik
clear('i')

%%%%%%%%%%%%%%%%%%%%%
% Equation Parameters
%%%%%%%%%%%%%%%%%%%%%
c =  0.0;  b =  6.0;  g =  0.0; d =  1.0;  

par = [c,b,g,d];

%%%%%%%
% Setup
%%%%%%%
%KdV 1-sol
L=30;				% Domain length
N = 1*128;			% Number of grid points
tfin=6.0;			% Final time

%%KdV 2-sol
%L=50;				% Domain length
%N = 2*128;			% Number of grid points
%tfin=10.0;			% Final time

dt = 0.002;			% Time step (if run becomes untable reduce dt)
Nsnaps = 50; 			% Number of snapshots to display on screen
Nsaves = 200;			% Number of solutions to save on allu
lw=2;				% linewidth for display

tsnaps = tfin/Nsnaps;		% Time between screen snapshots
tsaves = tfin/Nsaves;		% Time between saved solutions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% makes sure that dt divides tsnaps and tsaves
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(tsnaps/dt~=fix(tsnaps/dt))
 fprintf('dt=%f does not divide tsnaps=%f',dt,tsnaps);
 error('incompatible value for dt');
end
if(tsaves/dt~=fix(tsaves/dt))
 fprintf('dt=%f does not divide tsaves=%f',dt,tsaves);
 error('incompatible value for dt');
end

blurb = '\bf Wave equation: u_t + c u_x + \beta u u_x - \gamma u_{xx} + \delta u_{xxx} = 0';
disp_param = ['c=' num2str(c) ', \beta=' num2str(b) ', \gamma=' num2str(g) ', \delta=' num2str(d) ];


%%%%%%%%%%%%%%%%%%%
%  set grid vectors
%%%%%%%%%%%%%%%%%%%
ik = i*(2*pi/L)*(-N/2:N/2-1)';
fact = exp(dt*(c*ik-g*ik.^2+d*ik.^3));
dx = L/N;  xg = (0:dx:L-dx)';  x = [xg ; L];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial conditions and plotting window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%KdV 1-sol
% use v=1 and then v=2, v=4, v=8, ...
v = 1;
x0 = 3+8/(v+1);
ug  = (v/2)*sech((xg-x0)*(sqrt(v)/2)).^2;
pmin = -0.10; pmax = v/2+0.1;
pmax=4.1;

%%KdV 2-sol
%v1 = 4;
%% use v2=2 and then v2=3
%v2 = 2;
%x01 = 5;
%x02 = 15;
%ug1  = (v1/2)*sech((xg-x01)*(sqrt(v1)/2)).^2;
%ug2  = (v2/2)*sech((xg-x02)*(sqrt(v2)/2)).^2;
%ug = ug1 + ug2;
%pmin = -0.10; pmax = 2.10;

t = 0; 				% initial time
vg = fftshift(fft(ug))/N;	% FFT(initial condition)
allu = ug;			% this vector stores the saved snapshots
allt = 0;			% this vector stores the saving times
tsa = 0;			% saving loop-time
tsn = 0;			% displaying loop-time
u = [ug ; ug(1)];

%%%%%%%%%%%%%%%%%%%%%%%%
% Plot initial condition
%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
clf;
plot(x,u,'r',x,u,'k.','markersize',7);
axis([0 L pmin pmax]);
title(blurb);  
xlabel('\bf x-axis');  ylabel('\bf u-axis')
text(min(min(x))+L/20,pmax*0.92,disp_param);
text(max(max(x))-L/10,pmax*0.92,['t=' num2str(t)]);
pause(0.001)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quick check that initial condition is indeed periodic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this routine detects if the mas change of u happens at the boundary
if(abs(ug(1)-ug(N))>max(abs(diff(ug))))
 fprintf('I think your IC is not PERIODIC on the box [0:L] => Check plot!!!');
 error('Non-periodic IC');
end

input('press any key to continue');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  timestep loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for t = dt:dt:tfin

  [vg,tj] = wave_integrator_rk(vg,0,dt,par);	% pseudo-spectrol + Runge-Kutta method
  vg = vg./fact;
  ug = real(ifft(fftshift(vg))*N);

  tsa = tsa + dt;
  tsn = tsn + dt;
  if(tsa+dt/2>tsaves)		% it is time to save solution
   tsa=0;
   allu = [allu,ug];
   allt = [allt,t];
  end
  if(tsn+dt/2>tsnaps)		% it is time to display solution
   tsn=0;
   figure(1); clf;
   plot(x,[ug ; ug(1)],'b','LineWidth',lw);%  plot(x,[ug ; ug(1)],'k.');  
   axis([0 L pmin pmax]);
   title(blurb);  xlabel('\bf x-axis');  ylabel('\bf u-axis')
   text(min(min(x))+L/20,pmax*0.92,disp_param);
   text(max(max(x))-L/10,pmax*0.92,['t=' num2str(t)]);
   pause(0.001)
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of timestep loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot final condition together with initial
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on
plot(x,[allu(:,1);allu(1:1)],'r--');  
plot(x,[allu(:,1);allu(1:1)],'k.','markersize',7);
title(blurb);  xlabel('\bf x-axis');  ylabel('\bf u-axis')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calls the ploting subroutine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
input('press ANOTHER key to continue');
wave_plotsurf(allt,xg,allu,L,1,Nsaves+1,'no')