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



% This program calculates the FVH (Feynman Vernon Hellwarth) evolution for photon locking.
% 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 photon locking scheme the atom starts from the ground state.
% The field vector initially points to the negative 'x' axis: (-E, 0, 0).
% The Bloch vector evolves regularly (as in the simple FVH evolution), but
% when it reaches the negative 'y' axis the field phase is changed instantaneously so the field 
% vector points to the negative 'y' axis: (0, -E, 0). In this situation the Bloch vector, and the
% field vector are co-alligned. Thus there is no further evolution of the Bloch vector. If the 
% field phase is changed when the Bloch vector does not point exactly to
% the negative 'y' axis, the Bloch vector will process around the field
% vector with a small angle.

clear all;

% Initial parameters
E=1; % Field amplitude
D=0; % 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

% 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');
set(hplt(1),'linewidth',4);
set(hplt(1),'xdata',[0 -E]/a,'ydata',[0 0],'zdata',[0 D]/a);
set(hplt(2),'linewidth',4);


Mindex=1;
% FVH solutions
for j=1:27
    t=(j-1)*dt;
    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)));
    set(hplt(2),'xdata',[0 r1],'ydata',[0 r2],'zdata',[0 r3]);
    M(Mindex)=getframe(gcf);
    Mindex=Mindex+1;
    %pause(.1);
end;



% Photon locking
r=[r1;r2;r3];

% Change of phase of the Field vector
set(hplt(1),'xdata',[0 0]/a,'ydata',[0 -E]/a,'zdata',[0 D]/a);

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

% FVH solutions after change of phase
for j=1:207
    t=(j-1)*dt;
    r1=(c(2)*(i*a*exp(i*a*t))+c(3)*(-i*a*exp(-i*a*t)));
    r2=(c(1)*(-E)+c(2)*(D*exp(i*a*t))+c(3)*(D*exp(-i*a*t)));
    r3=(c(1)*(D)+c(2)*(E*exp(i*a*t))+c(3)*(E*exp(-i*a*t)));
    set(hplt(2),'xdata',[0 r1],'ydata',[0 r2],'zdata',[0 r3]);
    M(Mindex)=getframe(gcf);
    Mindex=Mindex+1;
    %pause(.1);
end;


