% This program computes the average position and momentum, the standard
% deviation for each of these, and produces some plots to help you
% visualize what is going on.  

% The program is written the Matlab language, but should also run on 
% the free alternatives: % Octave (www.octave.org) and Scilab (www.scilab.org).
% I actually use Octave. 

% I would strongly encourage you to learn some using one of these programs.
% If you continue on in any field of science, engineering, or mathematics
% they will be useful.  If you Googl "Matlab tutorial", you will find plenty
% of on-line introductions.


% Here are the inputs to the program.  Change these and see what happens.

% Range from xmin to xmax.  Computers do not work well with an infinite
% range.  We would like the wave function to be contained in this range
% so xmin and xmax may need to be adjusted.

xmin = -5;
xmax =  5;

% The number of points should be chosen so that wave function varies
% smoothly or slowly between adjacent points.

Npts = 1000;

% Now we are ready to define the wave function.  For this example
% I am choosing a Gaussian wavepacket which has 3 parameters,
% "width", "k", and "x0".  These are the things you probably 
% want to change first.

%%%%%%%%%%%%%%%%%% Part of Code to Play With %%%%%%%%%%%%%%%%%%%%%%%%%
width = 0.25;
k = 6;
x0 = 0;

dx = (xmax - xmin)/(Npts -1);
x = xmin + dx*(0:(Npts - 1));

% psi = exp(- (x - x0).^2/width^2) .* exp(i*k*x);
psi = exp(- (x - x0).^4/width^4) .* exp(i*k*x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% This wavefunction, psi, is not normalized.  Do so now.

psisquared = abs(psi).^2;
area = trapz(x,psisquared);
psi = psi/sqrt(area);

% Our first plot is of the |psi|^2 and the real and imaginary parts
% of psi.

figure(1)
plot(x,abs(psi).^2 ,x,real(psi),x,imag(psi));
legend('|psi|^2','real(psi)','imag(psi)');
xlabel('x')

% For our second plot look at x|psi|^2 and x^2 |psi|^2,
% and also compute <x>, <x^2>, sigmax .

figure(2)
plot(x,abs(psi).^2 ,x,x.*abs(psi).^2 ,x,(x.^2).*abs(psi).^2 )
legend('|psi|^2','x|psi|^2','x^2|psi|^2')
xlabel('x')

expvalx = trapz(x,x.*abs(psi).^2)
expvalxsq = trapz(x,(x.^2).*abs(psi).^2)
sigmax = sqrt(expvalxsq - expvalx.^2)

% Finally we computer the expectation values for the momentum,
% the momentum squared, and sigmap.  We set hbar equal to one.
% There are few subtleties here.  So this is not as straight
% forward as the position.

% start with p
momentumpsi = diff(psi)/i/dx;
psiinbetween = 0.5*(psi(1:(Npts -1)) + psi(2:Npts));
xinbetween = 0.5*(x(1:(Npts -1)) + x(2:Npts));
conjpsippsi = conj(psiinbetween).*momentumpsi;

expvalp = sum(conjpsippsi)*dx

% repeat for p^2
momentum2psi = diff(momentumpsi)/i/dx;
psiinbetween2 = 0.5*(psiinbetween(1:(Npts -2)) + psiinbetween(2:(Npts -1)));
xinbetween2 = 0.5*(x(1:(Npts -2)) + x(2:(Npts -1)));
conjpsip2psi = conj(psiinbetween2).*momentum2psi;

expvalp2 = sum(conjpsip2psi)*dx

sigmap = sqrt(real(expvalp2) - real(expvalp)^2)

figure(3)
plot(x,abs(psi).^2,xinbetween,real(conjpsippsi),xinbetween2,real(conjpsip2psi));
legend('|psi|^2','real(psi^*p psi)','real(psi^* p^2 psi)')

% Finally check the uncertainty principle

uncertaintyprinciple = sigmax * sigmap