% This program was written on 20/08/2003 by Erez Boukobza (Tannor's group).



% This program calculates the FVH (Feynman Vernon Hellwarth) evolution for adiabatic following.
% r1 & r2 represent coherences; r3 represents population differene between
% the excited & ground states (r3=-1 the atom is 100% in the ground state; 
% r3=1 the atom is 100% in the excited state).

% In the adiabatic following scheme the atom starts from the ground state.
% As the field vector changes adiabatically from: (-E, 0, -D) to (-E, 0, D) where D>>E 
% (E & D are the field amplitude and detuning respectively, both are
% positive here), the atom (Bloch vector) follows it by precessing around
% it constantly. Finally the population is almost completely inverted
% (almost all of it is in the excited state), and it remains there as long
% as the field is on.



clear all;

% Initial parameters
E=1; % Field amplitude
D=-50; % Detuning from resonance frequency
a=sqrt(D^2+E^2); % Rabi frequency
dt=2*pi/100; % time step for evolution 
z=[-1:0.01:1]; % coordinates for sphere environment
t=0;

% Initial condirions for the FVH evolution
r1_0=0;
r2_0=0;
r3_0=-1;
r_0=[r1_0;r2_0;r3_0];

% Matrix equation for obtaining the coefficients for the FVH solution
A=[-E D D;0 -i*a i*a;D E E];
c=A\r_0;

% Creation of sphere environment
[x,y,z]=sphere;
z(x<0&y<0)=inf;
surf(x,y,z);
xlabel('r_1''','position',[-.75 -2.3 -.35],'fontsize',12,'color','b');
ylabel('r_2''','position',[-2.2 -1.1 -.25],'fontsize',12,'color','b');
zlabel('r_3''','rotation',0,'position',[-2.4 -.25 .93],'fontsize',12,'color','b');
h=gcf;
set(h,'position',[50 50 350 300]);
set(h,'color','w');
h0=get(h,'Children');
h1=get(h0,'Children');
set(h1,'FaceAlpha',0.5);
set(h1,'FaceColor','interp');
set(h1,'EdgeAlpha',0.5);
hold on;

% Creation of Bloch vector, and Field vector
hplt=plot3([0 0],[0 0],[0 0],'b',[0 0],[0 0],[0 0],'r');%,0*z,0*z,z,'m',z,0*z,0*z,'m',0*z,z,0*z,'m'); 
%set(hplt(3),'linewidth',0.5);set(hplt(4),'linewidth',0.5);set(hplt(5),'linewidth',0.5);
set(hplt(1),'linewidth',4);
set(hplt(1),'xdata',[0 -E]/a,'ydata',[0 0],'zdata',[0 D]/a);
set(hplt(2),'linewidth',4);



% FVH solutions
indx1=1;
indx2=1;
for k=1:341
    set(hplt(1),'xdata',[0 -E]/a,'ydata',[0 0],'zdata',[0 D]/a);
    ch=k;
    for j=1:10
        t=dt*j;
        %t=10*(k-1)+dt*(j-1);
        r1=(c(1)*(-E)+c(2)*(D*exp(i*a*t))+c(3)*(D*exp(-i*a*t)));
        r2=(c(2)*(-i*a*exp(i*a*t))+c(3)*(i*a*exp(-i*a*t)));
        r3=(c(1)*(D)+c(2)*(E*exp(i*a*t))+c(3)*(E*exp(-i*a*t)));
        if (k>227)
            if mod(indx1,3)==0
                set(hplt(2),'xdata',[0 r1],'ydata',[0 r2],'zdata',[0 r3]);
                %drawnow;
                M3(indx2)=getframe(gcf);
                indx2=indx2+1;
            end;
        end;
        %set(hplt(2),'xdata',[0 r1],'ydata',[0 r2],'zdata',[0 r3]);
        %drawnow;
        %M(indx)=getframe;
        indx1=indx1+1;
    end;
        r=[r1;r2;r3]; % Initial conditions for the next FVH solution (new D)
    if D<-5  
        D=-50+k*1; % Detuning from resonance frequency
        a=sqrt(D^2+E^2); % Rabi frequency
        % Matrix equation for obtaining the coefficients for the FVH solution
        A=[-E D D;0 -i*a i*a;D E E];
        %[v,d] = eig(A);
        c=A\r;
    end;
    if (D>=-5 & D<5)
        D=-5+(k-45)*0.04; % Detuning from resonance frequency
        a=sqrt(D^2+E^2); % Rabi frequency
        % Matrix equation for obtaining the coefficients for the FVH solution
        A=[-E D D;0 -i*a i*a;D E E];
        %[v,d] = eig(A);
        c=A\r;
    end;
    if D>=5
        D=5+(k-295)*1; % Detuning from resonance frequency
        a=sqrt(D^2+E^2); % Rabi frequency
        % Matrix equation for obtaining the coefficients for the FVH solution
        A=[-E D D;0 -i*a i*a;D E E];
        %[v,d] = eig(A);
        c=A\r;
    end;
end;