
function rs =  eig_double_well(instruct,outfile)

% written by Einat Frishman, Quantum Molecular Dynamics, Weizmann, 1996.
% This program accompanies Ch. 9 of Introduction to Quantum Mechanics:
% A Time Dependent Perspective, by David Tannor (University Science Books, 2001).

% Uses the adiabatic switching method to calculate eigenfunctions for the double
% well, starting from the known analytical form of the eigenfunctions of the harmonic oscillator.
% This macro uses functions hpsi, poten, eigf
% parameters of the calculation

%clear

%omega = eval(instruct.omega); %--------------------------------------


xdim=1;                    % range of x dimension
nop=128;                   % no of points /2
xrange=nop/xdim;           %
mu=0.5*1836;               % the reduced mass
dt=2;                    % time step in each iteration
noiter=500;                % no of iteration
omega=0.1;            %--------------------------------------------
tinit=0.0;
potmax=3;
levelmax=12;
c=-i;

% total switching time is dt*noiter.
% dt shouldn't be too large, and dt*noiter should be as large as possible.

finaltime=dt*noiter;

v0=0.5*mu*omega*omega;
fprintf('No. of time steps: \t\t\t%g.\n',noiter)
x=[-1*nop:nop-1]/xrange;  % x range vector
n=length(x);

% x range vector:

x=[-1*nop:nop-1]/xrange;

%k-vector

k=[-1*nop:nop-1]*pi/xdim;
kine=1.*k.*k/(2.*mu);

alpha=mu*omega;
delk=sqrt(alpha)*5;
delx=7/sqrt(alpha);




for j=1:levelmax;
psi(j,:)=eigf(j,x,alpha);
norm=psi(j,:)*psi(j,:)'/xrange;
psi(j,:)=psi(j,:)/norm^0.5;

    pote=poten(mu,omega,x,tinit,finaltime);
%    [hpsinit(j,:),bb]=hpsi(psi(j,:),pote,kine);
%    einit=psi(j,:)*hpsinit(j,:)'/xrange;
%    fprintf('einit: \t\t%g.\n',einit)

psinit(j,:)=psi(j,:);

% PRINTING INITIAL WAVEFUNCTIONS:

f1 = figure(1);
subplot(3,4,j);
hold off;
plot(x,psinit(j,:));
xlabel('x (a.u.)')
ylabel('')
psisqr(j,:)=psi(j,:).*conj(psi(j,:));
maxsqr=max(psisqr(1,:));
maxx=sqrt(maxsqr)*1.2;
minx1=-delx;
maxx1=delx;
axis([minx1 maxx1 -maxx maxx]);
phshift=1;
psi(j,:)=psi(j,:).*exp(-i*pi*[-1*nop:nop-1]*phshift);
end;

%    pote1=pote(1,250)
%    psi1=psi(1,250)
%    hpsinit1=hpsinit(1,250)
%    fprintf('einit: \t\t%g.\n',einit)

text(-7,23,'Initial Wavefunctions')
%print single.ps;
tfin=tinit+dt*noiter;

% first step:
for j=1:levelmax;
  psi(j,:)=psi(j,:).*exp(-i*poten(mu,omega,x,0,finaltime)*dt/2);
fprintf('operator: \t\t%g%g.\n',exp(-i*poten(mu,omega,x(1),0,finaltime)*dt/2))
end;

% middle steps:

for iter=1:noiter
  t=iter*dt+tinit;

  for j=1:levelmax;
    m=fft(psi(j,:));


%  if (j==1);
%  subplot(3,3,8);
%   plot(k,m.*conj(m));
%  axis([-delk delk 0 max(m.*conj(m))]);
%  end;

    m=m.*exp(-i*kine*dt);
    b=ifft(m);
    pote=poten(mu,omega,x,t,finaltime);
    psi(j,:)=b.*exp(-i*pote*dt);
    [hpsif(j,:),bb]=hpsi(psi(j,:),pote,kine);
    E(j,iter)=psi(j,:)*hpsif(j,:)'/xrange;
%    fprintf('E=: %g.\n',E(j,iter))
%    plot(x,pote,'y-');
  end;

  if (iter/50==fix(iter/50));
  fprintf('iteration: %g.\n',iter);
  end;

end;


pause(1/100);

fprintf ('energy levels: \n')
for j=1:levelmax;
fprintf('level no. %g,%g.\n',j,E(j,noiter))
end

%last step:



for j=1:levelmax;
m=fft(psi(j,:));
m=m.*exp(-i*kine*dt/2);
psi(j,:)=ifft(m);
psisqr(j,:)=psi(j,:).*conj(psi(j,:));

% moving back psi in momentum space:

psi(j,:)=psi(j,:).*exp(i*pi*[-1*nop:nop-1]*phshift);


% PRINTING FINAL WAVEFUNCTIONS:
f2 = figure(2);
subplot(3,4,j);
hold off;

ata=atan(imag(psi(j,120))./real(psi(j,120)));
psifin(j,:)=real(psi(j,:)*exp(c*ata));

if (psifin(j,50)/psinit(j,50))<0
  psifin(j,:)=-psifin(j,:);
end;

norm=psi(j,:)*ones(n,1)/xrange;
plot(x,psifin(j,:),'r-');
axis([minx1 maxx1 -maxx maxx]);
xlabel('x (a.u.)')

end;
text(-7,23,'Final Wavefunctions')

end
%print double.ps;



f3 = figure(3);

hold off;
maxheight=poten(mu,omega,0,finaltime,finaltime)
vec=[1:noiter]*switch2(t,finaltime)*maxheight/noiter;
for i=1:levelmax;
plot(vec,real(E(i,:)));
hold on;
end;
xlabel('Height of Barrier (a.u.)')
ylabel('Energy (a.u.)')
title('Energy Levels')
%print energy.ps;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f4 = figure(4);

text(-3.2,13.5,'Final Potential and Energy Levels')
subplot(1,2,1);
hold off;
plot(x,poten(mu,omega,x,tinit,finaltime),'g-');
hold on;
for j=1:levelmax;
plot(x,real(E(j,1))*ones(1,n),'r-');
end
axis([minx1,maxx1,0,potmax])
xlabel('x (a.u.)')
ylabel('Initial Potential (a.u.)')
subplot(1,2,2);
hold off;
plot(x,poten(mu,omega,x,tfin,finaltime),'g-');
hold on;
for j=1:levelmax;
plot(x,real(E(j,noiter))*ones(1,n),'r-');
end
axis([minx1,maxx1,0,potmax])
xlabel('x (a.u.)')
ylabel('Final Potential (a.u.)')
%print potentials.ps;

% estimating energy difference:

phi1=psifin(1,:);
phi2=psifin(2,:);
phi_right=(phi1-phi2)/sqrt(2);
phi_left=(phi1+phi2)/sqrt(2);

overlap=phi_right*phi_left'/xrange

[y,ym]=min(poten(mu,omega,x,tfin,finaltime));

% Distance between minima:
dx=2*abs(x(ym));

% Height of maximum:

v=1.0-y(1);
mom=sqrt(2*mu*v);

lambda=dx*mom;
enerdiff=0.05*exp(-lambda)
d12=E(2,noiter)-E(1,noiter)
d34=E(4,noiter)-E(3,noiter)
d56=E(6,noiter)-E(5,noiter)

d56/d34

  return;
end;



