delt = .01;
endt = input('End time?');
delx = .03;
t = 0:delt:endt;
tmesh=size(t,2);
a = 5.1;
x = -a:delx:a;
s = input('Initial velocity?');
spread = input('How many samples?  The nth sample will have velocity v(n-1) + 1.') - 1;
m = 1;
k = s:(spread+s);
vr = input('What is the recoil velocity?  .6 for rubidium.');
xmesh = size(x,2);
j = 1:xmesh;
p = 1:xmesh;
ymin = -1;
ymax = 1;
d = ones(1,xmesh);
e = ones(xmesh,1);
%l = ymin:.2:ymax;
pot = input('Enter a potential function of x.  For a random potential, enter 999.');
if pot == 999
    pot = d;
    fluc = input('How large are the fluctuations?')
    while m<xmesh+1;
        pot(m) = pot(m)*2*rand*fluc;
    m = m+1;
end
pot = pot-fluc;
m = 1;
end
%for j = 1:xmesh;
%    if x(j) == 0;
%pot(j)=0;
%end
e = ones(xmesh,1);
K = spdiags([-e (2+(delx^2)*pot').*e -e], -1:1, xmesh, xmesh);
c=delt/delx^2;
sigma=.3;
K(xmesh,1) = -1;
K(1,xmesh) = -1;
A = zeros(xmesh,(spread+1));
nv = ones((spread+1),1);
fc = zeros((spread+1),1);
figure
while m<spread+2;
    z = abs(t-(a/(2*k(m))));
[O,W]=min(z);
    n=1;
    u = (2/(pi*sigma^2))^.25.*exp(-(x).^2/(sigma^2)).*(exp(i*(k(m)+vr).*x)+exp(i*(k(m)-vr).*x));
    %v = 2*cos(k*x).*exp(i*vr*x);
    %u = exp(i*2.*x)+exp(i*-2.*x)+exp(i*3.*x)+exp(i*-3.*x)
while n<tmesh;
    y = (abs(u)).^2;
    subplot(3,1,1); plot(x,[real(u); imag(u)])
    axis([-a a ymin ymax])
    set(gca,'xtick',[-10:2:10])
    subplot(3,1,2); plot(x,y)
    set(gca,'xtick',[-10:2:10])
    axis([-a a 0 2*ymax])
    u = ((speye(xmesh) + .5i*c.*K)*((speye(xmesh) - .5i*c.*K)\(u')))';
   % v = ((speye(xmesh) + .5i*c.*K)*((speye(xmesh) - .5i*c.*K)\(v')))';
    n = n + 1;
    text(0,1.5,sprintf('t=%g',t(n)));
    drawnow
    if n == W
        q = u;
    end
end
yprime=y';
f = round((xmesh/2) - (pi/(k(m)*delx)));
g = round((xmesh/2) + (pi/(k(m)*delx)));
cut = q(f:g);
A(j,m)=A(j,m)+yprime(j);
avg = (A*nv)./m;
Imax = max(cut);
Imin = min(cut);
fc(m,1) = (abs(Imax-Imin))/(abs(Imax+Imin));
subplot(3,1,3);plot(x,avg)
axis([-a a 0 2*ymax])
text(-1,1.5,sprintf('fringe contrast around x = 0 at t = a/(2k) = %g',fc(m)));
m=m+1;
pause
end