%  This file will compute the time-dependent Schrodinger equation in 2-D starting with an initial gaussian wavepacket input.
%  This file contains realistic numbers for the constants.

clear all;

format long;


%  Initialize all of the parameters


i_imag = sqrt(-1);

N = 200;

L = 100*10^-9;	% in meters

h = L/(N-1);

x = h*(0:N-1) - L/2;

y = h*(0:N-1) - L/2;

h_bar = 1.055*10^-34;

m_gaas = .067*9.11*10^-31;

m_algaas = .09*9.11*10^-31;

tau = 2e-16;

V0 = .1*1.6*10^-19;


%  Set up the Hameltonian matrix

ham = zeros(N);

coeff = -h_bar^2/(2*m_gaas*h^2);

coeff2 = -h_bar^2/(2*m_algaas*h^2);

% Define the x-directed Hamiltonian

for ii=2:1:N-1
    
    hamx(ii,ii-1) = coeff;
    
    hamx(ii,ii)= -2*coeff;
    
    hamx(ii,ii+1) = coeff;
    
end

%  Define the first and last points

hamx(1,1) = -2*coeff;

hamx(1,2) = coeff;

hamx(N,N) = -2*coeff;

hamx(N,N-1) = coeff;

hamx(N,1)=coeff;

hamx(1,N)=coeff;

%Left side

% Define the y-directed Hamiltonian

for ii=2:1:N-1
    
    hamy(ii,ii-1) = coeff;
    
    hamy(ii,ii)= -2*coeff;
    
    hamy(ii,ii+1) = coeff;
    
    if((ii>95)&(ii<105))
        
        hamy(ii,ii-1) = coeff;
    
        hamy(ii,ii)= -2*coeff+V0;
    
        hamy(ii,ii+1) = coeff;
        
    end
    
end

%  Define the first and last points

hamy(1,1) = -2*coeff+V0;

hamy(1,2) = coeff;

hamy(N,N) = -2*coeff+V0;

hamy(N,N-1) = coeff;

hamy(N,1)=coeff;

hamy(1,N)=coeff;

%  Set up first and last rows for periodic boundary conditions  for finite well

hamy(1,N)=coeff;

hamy(1,1)=(-2*coeff)+V0;

hamy(1,2)=coeff;

hamy(N,N-1)=coeff;

hamy(N,N)=(-2*coeff)+V0;

hamy(N,1)=coeff;

%  Compute the Crank-Nicholson matrix

dCNx =(inv(eye(N)+.5*i_imag*tau/h_bar*hamx)*(eye(N)-.5*i_imag*tau/h_bar*hamx));
dCNy =(inv(eye(N)+.5*i_imag*tau/h_bar*hamy)*(eye(N)-.5*i_imag*tau/h_bar*hamy));

%  Initialize the wavefunction  in the x-direction

x0 = -25e-9;

velocity = 7.24*10^5;    	% in m/s

k0 = m_gaas*velocity/h_bar;

sigma0=L/10;

psi=exp(i_imag*k0*x').*exp(-(x'-x0).^2/(2*sigma0^2));

Norm_psi=sum(psi.*conj(psi));

psi = (1/sqrt(Norm_psi))*psi;

%  Plot the initial wavefunction

figure(1);

clf;

plot(x,real(psi),'-',x,imag(psi),'--');

title('Initial wave function');

xlabel('x');

ylabel('\psi(x)');

legend('Real','Imag');

drawnow;

pause(1);


%  Initialize loop and plot variables

%max_iter=L/(velocity*tau);

max_iter=50000;

plot_iter=max_iter/500;

p_plot(:,1)=psi.*conj(psi);

iplot=1;

figure(2);

clf;

axisV=[-L/2 L/2 0 max(p_plot)];


%  Loop over the desired number of steps

for iter=1:max_iter
   
   psi=dCNx*psi;
   
   %  Record plotting values
   
   if (rem(iter,plot_iter)<1)
      
      iplot=iplot+1;
      
      p_plot(:,iplot)=psi.*conj(psi);
      
      plot(x,p_plot(:,iplot));
      
      xlabel('x');
      
      ylabel('P(x,t)');
      
      title(sprintf('Finished %g of %g iterations',iter,max_iter));
      
      axis(axisV);
      
      drawnow;
      
   end
   
end

%  Plot probability versus position at various times

%pxFinal=psi.*conj(psi);

%plot(x,p_plot(:,1:2:iplot),x,pxFinal);

%xlabel('x');

%ylabel('P(x,t)');

%title('Probability density at various times');


%  Now for the y-direction

y0 = 25e-9;

velocity = 0*10^6;		% in m/s

k0 = m_gaas*velocity/h_bar;

sigma0=L/10;

psi_y=exp(i_imag*k0*y').*exp(-(y'-y0).^2/(2*sigma0^2));

Norm_psi_y=sum(psi_y.*conj(psi_y));

psi_y = (1/sqrt(Norm_psi_y))*psi_y;



%  Plot the initial wavefunction

figure(3);

clf;

plot(y,real(psi_y),'-',y,imag(psi_y),'--');

title('Initial wave function');

xlabel('y');

ylabel('\psi(y)');

legend('Real','Imag');

drawnow;

pause(3);


%  Initialize loop and plot variables

%max_iter=L/(velocity*tau);

max_iter=50000;

plot_iter=max_iter/500;

py_plot(:,1)=psi_y.*conj(psi_y);

iplot=1;

figure(4);

clf;

axisV=[-L/2 L/2 0 max(py_plot)];


%  Loop over the desired number of steps

for iter=1:max_iter
   
   psi_y=dCNy*psi_y;
   
   %  Record plotting values
   
   if (rem(iter,plot_iter)<1)
      
      iplot=iplot+1;
      
      py_plot(:,iplot)=psi_y.*conj(psi_y);
      
      plot(y,py_plot(:,iplot));
      
      xlabel('y');
      
      ylabel('P(y,t)');
      
      title(sprintf('Finished %g of %g iterations',iter,max_iter));
      
      axis(axisV);
      
      drawnow;
      
   end
   
end

%  Plot probability versus position at various times

%pyFinal=psi_y.*conj(psi_y);

%plot(y,py_plot(:,1:2:iplot),y,pyFinal);

%xlabel('y');

%ylabel('P(y,t)');

%title('Probability density at various times');




%  Find the final wavepacket

for i=1:1:iplot
   
   for j=1:1:N
      
      for k=1:1:N
   
   pFinal(k,j,i) = p_plot(j,i)*py_plot(k,i);
   
end

end

end

%  Plot the probability versus position 

figure(5);

smax = max(max(pFinal(:,:,1)));

axis([-L/2 L/2 -L/2 L/2 0 smax]);

for i=1:1:iplot
   
   	surf(x/1e-9,y/1e-9,pFinal(:,:,i));
   
   	shading interp;
      
      xlabel('x');
      
      ylabel('y');
      
      zlabel('P(y,t)');
      
      title(sprintf('Finished %g of %g iterations',i,iplot));
      
      drawnow;
      
   end
   
   
   figure(6);

axis([-L/2 L/2 -L/2 L/2]);

caxis([0 smax]);

for i=1:1:iplot
   
   	pcolor(x/1e-9,y/1e-9,pFinal(:,:,i));
    
    caxis([0 smax]);
   
	   shading interp;
   
   	colorbar('h');
      
      xlabel('x');
      
      ylabel('y');
           
      title(sprintf('Finished %g of %g iterations',i,iplot));
      
      drawnow;
      
   end










   
   
