% Ploting a movie of a 2D RATIONAL KdV solution

if(exist('plotvel')~=1)
 error('Please first define plotvel, (plotvel=1 : show velocity, plotvel=0 : do not show velocity)');
end

x=-5:0.2:20;
y=-5:0.1:5;
t0=0; dt=0.05; t1=1.5;
[X,Y]=meshgrid(x,y);

xz=-3:0.01:2;
yz=-2:0.01:2;
[Xz,Yz]=meshgrid(xz,yz);

figure(1)
clf
set(gca,'FontSize',[16]);

p=2;
t=t0;

ln=2;
vi=[24 72];

xx=-1/p;
yy=0;
vx=3*p^2;
vy=0;
factor=5;
vx=vx/factor;
vy=vy/factor;

XX = X+(1/p)-3*p^2*t;
Z = 4*(p^2*Y.^2-XX.^2+(1/p^2))./(p^2*Y.^2+XX.^2+(1/p^2)).^2;
MZ = max(max(Z));
dZ = MZ/100;

XX = Xz+(1/p)-3*p^2*t;
Zz = 4*(p^2*Yz.^2-XX.^2+(1/p^2))./(p^2*Yz.^2+XX.^2+(1/p^2)).^2;
MZz = max(max(Zz));
dZz = MZ/100;

surfl(Xz,Yz,Zz);
shading interp
colormap(gray)
title(['2D soliton solution for the KdV; t=' num2str(t)])
ylabel('y')
xlabel('x')
zlabel('u(x,y)')
view([15 50])
axis('tight');
ax=axis;
if(plotvel)
 hold on
 plot3([xx xx+vx],[yy yy],MZz*[1 1],'r:','LineWidth',ln)
 plot3([xx xx],[yy yy+vy],MZz*[1 1],'r:','LineWidth',ln)
 plot3([xx xx+vx],[yy yy+vy],MZz*[1 1],'r-','LineWidth',ln)
 plot3([xx-10*vx xx+10*vx],[yy-10*vy yy+10*vy],MZz*[1 1],'r:','LineWidth',ln)
 plot3(xx,yy,MZz+dZz,'k.','MarkerSize',20)
 hold off
end


figure(2)
clf
set(gca,'FontSize',[16]);

for t=t0:dt:t1
 XX = X+(1/p)-3*p^2*t;
 Z = 4*(p^2*Y.^2-XX.^2+(1/p^2))./(p^2*Y.^2+XX.^2+(1/p^2)).^2;
 surfl(X,Y,Z);
 shading interp
 colormap(gray)
 title(['2D soliton solution for the KdV; t=' num2str(t)])
 ylabel('y')
 xlabel('x')
 zlabel('u(x,y)')
 view(vi)
 axis tight;
 if(plotvel)
  hold on
  plot3([xx xx+vx],[yy yy],MZ*[1 1],'r:','LineWidth',ln)
  plot3([xx xx],[yy yy+vy],MZ*[1 1],'r:','LineWidth',ln)
  plot3([xx xx+vx],[yy yy+vy],MZ*[1 1],'r-','LineWidth',ln)
  plot3([xx-10*vx xx+10*vx],[yy-10*vy yy+10*vy],MZ*[1 1],'r:','LineWidth',ln)
  plot3(xx,yy,MZ+dZ,'k.','MarkerSize',20)
  hold off
  xx=xx+dt*vx*factor;
 end
 pause(0.001)
 if(t==t0)
  fprintf('Press any key...\n');
  figure(2)
  set(gca,'FontSize',[16]);
  pause
 end
end