% PIB: Particle in a box
a =  1e-9 ;% one nanometer  length of the box
hbar = 1.e-34  ;%  Joule Sec
me= 1e-30  % mass of the electron
n= [ 1 2  3 4];
Npts=100
k = pi*n/a;
E = (hbar*k).^2/(2*me);
x = linspace(0,a,Npts);
dx = x(2)-x(1);
Amp = sqrt(2/a);
Psi = Amp*sin((pi*n')*(x/a));
SF = 1e24   % This is the E to Amplitude Scale Factor  A/E
SF = Amp/((hbar*pi/a).^2/(2*me));
figure(1)
hold on
for k=1:4
    plot(x,Psi(k,:)+E(k)*SF)  % ,x,Psi(2,:)+E(2)*SF,x,Psi(3,:)+SF*E(3))
end
xlabel('position [meters]');, ylabel('Wave Amplitude')
disp('paused after Fig 1')
pause
% Look at orthogonality pick two functions mult together and do a cumsum
npick = [ 1 2 ]
Ppick =  sqrt(2/a)*sin((pi*npick')*(x/a));
Pprod =  prod(Ppick,1)
PC = cumsum(Pprod);
figure(2)
plot(x,Ppick(1,:),x,Ppick(2,:),x,Pprod/Amp,x,PC/Amp/sqrt(Npts))
xlabel('position [meters]');, ylabel('Wave Amplitude')
NormPsi = (Ppick*Ppick')*dx



% put some time dependence onto the wave function.  
% Use the ones chosen in npick
kv = pi*npick/a;
E = (hbar*kv).^2/(2*me);

% Need the momentum operated on wave functions:  P*Phi(n)  Mom-Psi
MC = -i*hbar*(kv'*ones(size(x)));
MomPsi = Amp*(MC.*cos((pi*npick')*(x/a)));

Wfreq = (E-E(1))/hbar;
RelAmp = [ 1  1]
RelAmp = RelAmp/sqrt(RelAmp*RelAmp')
t=0
Psitime = (exp(-i*Wfreq*t).*RelAmp)*Ppick
Probtime = conj(Psitime).*Psitime;
Meanx = dx*(x*Probtime');
SqX =dx*((x.^2)*Probtime')
Rmsx = sqrt(SqX);
Var = SqX - Meanx^2;
% convert Meanx and Rmsx into a gaussian curve to see
y = (x-Meanx).^2/(2*Var);
gx = Amp*exp(-y);
MeanMom = dx*((exp(-i*Wfreq*t).*RelAmp)*MomPsi)*Psitime';
figure(3)

hp=plot(x,Psitime,x,Probtime/(Amp),Meanx,Amp,'k*',Rmsx,Amp/2,'ks',x,gx,'g--')
xlabel('position [meters]');, ylabel('Wave Amplitude')
Tmax = 4*2*pi/Wfreq(2);
set(gca,'ylim',[-9 9]*1e4)
Ntime=1000
disp('paused after Fig 3 before movie')
pause
PosTime =[];
MomTime =[];
t = linspace(0,Tmax,Ntime)
for k=1:Ntime
    % where is the mean position x  xrms is the root mean square position
    Psitime = (exp(-i*Wfreq*t(k)).*RelAmp)*Ppick;
    Probtime = conj(Psitime).*Psitime;
    Ptime = Probtime/(Amp);
    Meanx = dx*(x*Probtime');
    SqX =dx*((x.^2)*Probtime')
    Rmsx = sqrt(SqX);
    Var = SqX - Meanx^2;
    y = (x-Meanx).^2/(2*Var);
    gx = Amp*exp(-y);
    % retain mean position and momentum as a function of time:
    MeanMom = dx*((exp(-i*Wfreq*t(k)).*RelAmp)*MomPsi)*Psitime';
    PosTime =[PosTime Meanx];
    MomTime =[MomTime MeanMom];

    %plot(x,Psitime,x,Probtime/(Amp),Meanx,Amp,'k*',Rmsx,Amp/2,'ks')
    set(hp(1),'ydata',Ptime);, set(hp(2),'ydata',Psitime);
    set(hp(3),'xdata',Meanx);, set(hp(4),'xdata',Rmsx);, set(hp(5),'ydata',gx)
    pause(.03)
end

% now can plot mean and momentum vs time and parameterically vs each other
figure(4)
subplot(3,1,1), plot(t,PosTime)
xlabel('time [sec]'),ylabel('position [meters]')
subplot(3,1,2), plot(t,MomTime/me)
xlabel('time [sec]'),ylabel('velocity [meters/sec]')
subplot(3,1,3), plot(PosTime, MomTime/me)
xlabel('position [meters]'),ylabel('velocity [meters/sec]')