%
%  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
%
%  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')
fs=16;

mycase=1; % 1:orig, 2:HW#9.2, 3:1 KdV soliton

%%%%%%%%%%%%%%%%%%%%%
% Equation Parameters
%%%%%%%%%%%%%%%%%%%%%
switch mycase
 case 1 %%Orig
  c =  4.0;  b =  0.0;  g =  0.75; d =  1.0;  

 case 2 %%HW#9.2
  c =  4.0;  b =  0.0;  g =  0.25; d =  1.0;  

 case 3 %KdV 1-sol
  c =  0.0;  b =  6.0;  g =  0.0; d =  1.0;  
end

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

%%%%%%%
% Setup
%%%%%%%
switch mycase
 case 1 %%Orig
  k = 1.0;			% wave number
  Nperiods=4;			% how many periods is the domain
  L=Nperiods*2*pi/k;		% Domain length
  N = 128;			% Number of grid points
  tfin=4.0;			% Final time
  dt = 0.01;			% Time step (if run becomes untable reduce dt)

 case 2 %%HW#9.2
  k = 1.0;			% wave number
  Nperiods=3;			% how many periods is the domain
  L=Nperiods*2*pi/k;		% Domain length
  tfin=6.0;			% Final time
  dt = 0.015;			% Time step (if run becomes untable reduce dt)

 case 3 %%KdV 1-sol
  k = 1.0;                        % wave number
  Nperiods=4;                     % how many periods is the domain
  L=Nperiods*2*pi/k;              % Domain length
  N = 128;                        % Number of grid points
  tfin=10.0;			% Final time
  dt = 0.0025;			% Time step (if run becomes untable reduce dt)
end


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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch mycase
 case 1 %%Orig
  A = 1;				% Amplitude
  ug  = A*cos(k*xg);		% Initial condition
  pmin = -1.5; pmax = 1.5;	% Y-axis display window
  %%add more modes
  %kk=[1:10];
  %AA=1*rand(1,10);
  %for ii=1:5
  % ug = ug+AA(ii)*cos(kk(ii)*xg+2*pi*rand(1));		% Initial condition
  %end

 case 2 %%HW#9.2
  ug  = 3*cos(xg).^2+sin(xg);
  pmin = -1.25; pmax = 3.5;

 case 3 %%KdV 1-sol
  v = 2;
  ug  = (v/2)*sech((xg-L/2)*(sqrt(v)/2)).^2; 
  pmin = -0.10; pmax = 1.10;
  %%add pertrurbation
  %ug = ug+0.1*max(ug)*(rand(N,1)-1*0.5);
end

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;
set(gca,'FontSize',fs)
plot(x,u,'r',x,u,'k.','markersize',7);
axis([0 L pmin pmax]);
title(blurb);  
xlabel('\bf x-axis','FontSize',fs);  ylabel('\bf u-axis','FontSize',fs)
text(min(min(x))+L/20,pmax*0.92,disp_param,'FontSize',fs);
text(max(max(x))-L/10,pmax*0.90,['t=' num2str(t)],'FontSize',fs);
pause(0.1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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-spectral + 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;
   set(gca,'FontSize',fs)
   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','FontSize',fs);  ylabel('\bf u-axis','FontSize',fs)
   text(min(min(x))+L/20,pmax*0.92,disp_param,'FontSize',fs);
   text(max(max(x))-L/10,pmax*0.90,['t=' num2str(t)],'FontSize',fs);
   pause(0.1)
  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')
hold off
hleg=legend(['  t=',num2str(tfin),'  '],'  t=0  ');
set(hleg,'FontSize',fs);

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


