%Time-evolution of a wavepacket using split operator method
%with FFT
clear all;

n=1024;  %number of points in spatial grid Choose power of 2 for speed
xmax=20.0;
xmin=-1*xmax;
dx=(xmax-xmin)/(n-1);
for in=1:n
	 x(in)=xmin+dx*in;
pot(in)=0.0;  %No potential for now
%pot(in)=0.5*x(in)^2;  %Harmonic potential
end

%Set the grid in k-space
%Not the funny order imposed by the ordering in fast F-transform
deltak=2*pi/((n-1)*dx);
for in=1:n/2
	 kv(in)=(in-1)*deltak;
kv(n-in+1)=-1*in*deltak;
end

	 for in=1:n
		  Hk(in)=0.5*kv(in)^2;
end



dt=dx^2/4.0;  %need this timestep scaling for numerical stability
time_final=15;
samples=100;
tnum=floor(time_final/(dt));
jump=floor((tnum*1.0/samples));
i=complex(0,1);

%Initialize wavefunction
width=1;  
k0=2.0;
%k0=0.0;
pos=-5.0;
wf=exp(-0.5*(x-pos).^2/width^2).*exp(i*k0*x);
%wf=(x-pos).*exp(-0.5*(x-pos).^2/width^2).*exp(i*k0*x);
%wf=sin(2*x).*exp(-0.5*(x-pos).^2/width^2).*exp(i*k0*x);

wf=wf/sqrt(sum(abs(wf).^2)*dx);


    fignum=1;fs=20;
figure(fignum);clf;
th=0;
subplot(2,1,1);plot(x,abs(wf).^2,'LineWidth',2);grid on;
xlabel('x','FontSize',fs);
titletext=['t= ',num2str(th,'%5.2f')];
title(titletext,'FontSize',fs);
a=gca;set(a,'FontSize',fs);
axis([xmin xmax 0 0.5]);

subplot(2,1,2);
data=angle(wf);data=unwrap(data);
plot(x,data,'LineWidth',2);grid on;
xlabel('x','FontSize',fs);xlabel('\phi','FontSize',fs);
titletext=['t= ',num2str(th,'%5.2f')];
title(titletext,'FontSize',fs);
a=gca;set(a,'FontSize',fs);
pause(0.2);





%Time loop
for ti=1:tnum
	 th=dt*ti;

%Split operator time-step
	 wf=exp(-i*dt*pot).*wf;
wf=fft(wf)/sqrt(n*1.0);
wf=exp(-i*dt*Hk).*wf;
  wf=ifft(wf)*sqrt(n);

%Plot every now and then
  if (mod(ti,jump)==0)
    fignum=1;fs=16;
figure(fignum);clf;

subplot(2,1,1);plot(x,abs(wf).^2,'LineWidth',2);grid on;
xlabel('x','FontSize',fs);
titletext=['Prob. density: t= ',num2str(th,'%5.2f')];
title(titletext,'FontSize',fs);
a=gca;set(a,'FontSize',fs);
axis([xmin xmax 0 0.7]);

subplot(2,1,2);
data=angle(wf);data=unwrap(data);
plot(x,data,'LineWidth',2);grid on;
xlabel('x','FontSize',fs);ylabel('\phi','FontSize',fs);
titletext=['angle: t= ',num2str(th,'%5.2f')];
title(titletext,'FontSize',fs);
a=gca;set(a,'FontSize',fs);


time=th
  pause(0.2);
end

end