function WavePacketDoppler

%Normal Doppler

xo = 0.0; %	xo = wave-packet center [def.=0]
Dx = 10.; %	Dx = standard deviation from xo of the probability distribution [def.=0.9 A]
Dx = 3.; %	Dx = standard deviation from xo of the probability distribution [def.=0.9 A]
Ko = -15.0;%	Ko = averaged momentum/h_bar (wave-number) [def.= 0 A^-1]
Ko = -15.0;%	Ko = averaged momentum/h_bar (wave-number) [def.= 0 A^-1]
relMass = 1.0; %	relMass = particle mass/electron mass ratio [def.=1].	
xWidth = .10; %	xWidth = amplitude of the x-axis interval [def.=40 A]
pNum = 2001; %	pNum = (odd) number of points covering the x-axis interval [def.=1001].
Wo = .20; % Hz-1
t = 1.;

gcf = figure('Units','inches','Position',[-.25 0.5 11 3.5],'Color',[1 1 1])
%axes('fontsize',axissize,'position',[.15 .15 .8 .73]);
hold off;
cla; % clear the screen
hold on;
axis([-100 350 -.25 .25]);
axis off;

if Ko==0,
	lambda= 1e16; % wave length in Angstrom
else
	lambda= 2*pi/Ko;
end 
N=(pNum-1)/2;
xWidth = pNum;
x=xWidth*(-N:N)/pNum;
x=x(:);

A = 0.7;
A2 = A^2;
aw0 = 0.005;
ak0 = 0.05;
aDw = -0.01;
aDk = 0.85;
t = 0.0;
aDw1 = aDw;
avs = 0.03;
%Dk = Dk*(Dw/(Dw+vs))
%k0 = k0*(Dw/(Dw+vs));

tburst1 = 2000;
tburst2 = 4000;
tburst3 = 6000;
tburst4 = 8000;
first1 = 1;
first2 = 1;
first3 = 1;
first4 = 1;

mov = moviein(500);
framecount = 1;

for t = 1:20:10000
    xxx = avs*t;
   cg = abs(Wo/Ko);
aDk1 = aDk*(cg/(cg-avs));
aDk2 = aDk*(cg/(cg+avs));
%reverse
%   f1 =     A2*sin((aw0-2.0*aDw1)*t - (ak0+2.0*aDk1)*x);
%   f2 =     A2*sin((aw0+2.0*aDw1)*t + (ak0+2.0*aDk2)*x);
%forward direction
   f1 =     A2*sin((aw0-2.0*aDw1)*t + (ak0+2.0*aDk1)*x);
   f2 =     A2*sin((aw0+2.0*aDw1)*t - (ak0+2.0*aDk2)*x);

   xo = (cg+avs)*t;
   x2=(x-xo).^2;
   Dx2=Dx*Dx; %+0.01*t;
   %A=1.0/sqrt(1+t^2/Dx2^2);
   Psi1 = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi1.*Psi1));
   Psi1=Psi1/sqrt(Norm); % normalized wave function
   Psi1 = Psi1.*f1;
   
   xo = (cg)*t;
   x2=(x+xo).^2;
   Dx2=Dx*Dx;
   Psi2 = exp(-x2/(4*Dx2));
   

   iK=sqrt(-1)*Ko;	% otherways Psi is real
%iW=sqrt(-1)*Wo;
%Psi = Psi.*exp(iK*x);
   
   Norm = sum(abs(Psi2.*Psi2));
   Psi2=Psi2/sqrt(Norm); % normalized wave function
Psi2 = Psi2.*f2;
   
if t > tburst1
    if (first1 == 1)
        xoa1 = xxx;
        ta1 = t;
        first1 = 0;
    end
   xcenter = xoa1 + (cg+avs)*(t-ta1); 
   x2=(x-xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi1a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi1a.*Psi1a));
   Psi1a=Psi1a/sqrt(Norm); % normalized wave function
   Psi1a = Psi1a.*f1;
   Psi1 = Psi1+Psi1a;  
   
   xcenter = xoa1 - (avs-cg)*(t-ta1); 
   x2=(-x+xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi2a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi2a.*Psi2a));
   Psi2a=Psi2a/sqrt(Norm); % normalized wave function
   Psi2a = Psi2a.*f2;
   Psi2 = Psi2+Psi2a;  
   
   
end
if t > tburst2
    if (first2 == 1)
        xoa2 = xxx;
        ta2 = t;
        first2 = 0;
    end
   xcenter = xoa2 + (cg+avs)*(t-ta2); 
   x2=(x-xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi1a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi1a.*Psi1a));
   Psi1a=Psi1a/sqrt(Norm); % normalized wave function
   Psi1a = Psi1a.*f1;
   Psi1 = Psi1+Psi1a;  
   
   xcenter = xoa2 - (avs-cg)*(t-ta2); 
   x2=(-x+xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi2a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi2a.*Psi2a));
   Psi2a=Psi2a/sqrt(Norm); % normalized wave function
   Psi2a = Psi2a.*f2;
   Psi2 = Psi2+Psi2a;  
   
end
if t > tburst3
    if (first3 == 1)
        xoa3 = xxx;
        ta3 = t;
        first3 = 0;
    end
   xcenter = xoa3 + (cg+avs)*(t-ta3); 
   x2=(x-xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi1a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi1a.*Psi1a));
   Psi1a=Psi1a/sqrt(Norm); % normalized wave function
   Psi1a = Psi1a.*f1;
   Psi1 = Psi1+Psi1a;  
   
   xcenter = xoa3 - (avs-cg)*(t-ta3); 
   x2=(-x+xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi2a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi2a.*Psi2a));
   Psi2a=Psi2a/sqrt(Norm); % normalized wave function
   Psi2a = Psi2a.*f2;
   Psi2 = Psi2+Psi2a;  
   
end
if t > tburst4
    if (first4 == 1)
        xoa4 = xxx;
        ta4 = t;
        first4 = 0;
    end
   xcenter = xoa4 + (cg+avs)*(t-ta4); 
   x2=(x-xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi1a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi1a.*Psi1a));
   Psi1a=Psi1a/sqrt(Norm); % normalized wave function
   Psi1a = Psi1a.*f1;
   Psi1 = Psi1+Psi1a;  
   
   xcenter = xoa4 - (avs-cg)*(t-ta4); 
   x2=(-x+xcenter).^2;
   Dx2=Dx*Dx; %+0.01*t;
   Psi2a = exp(-x2/(4*Dx2));
   Norm = sum(abs(Psi2a.*Psi2a));
   Psi2a=Psi2a/sqrt(Norm); % normalized wave function
   Psi2a = Psi2a.*f2;
   Psi2 = Psi2+Psi2a;  
   
end


   if t == 1
      hd1=plot(x,Psi1,'LineWidth',2);
      hd2=plot(x,Psi2,'color',[0 .7 0],'LineWidth',2);
      hd3=plot(xxx,0,'.r','MarkerSize',25);
   else
      set(gcf,'DoubleBuffer','on');
      set(hd1,'ydata',real(Psi1));
      set(hd2,'ydata',real(Psi2));
      set(hd3,'xdata',xxx);
      drawnow    
   end
   
   mov(framecount) = getframe(gcf);
   framecount = framecount + 1;
   
end

movie2avi(mov, 'Doppler_R.avi','fps', 30,'compression','None','quality', 85 );
