%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Monte-carlo wavefunction simulation for a two-level atom [J. Dalibard et al., PRL68,0580 (1992)].
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MonteCarlo
    figure(1);						clf;

    G =2*pi*6e6;
    W =     5*G;
    D =    .5*G;
    dt=   .5e-9;
    t =[0:400]*dt;

    M =MonteCarlo_Bloch(G,W,D);
    rho0=[1 0 0 0]';
    for jj=1:length(t)
        rho=expm(M*t(jj))*rho0;
        rh1(jj)=real(rho(2));
    end

    H =MonteCarlo_Hamilton(G,W,D);
    nn1=zeros(size(t));
    for jj=1:200
        Y =[1 0]';
        for j=1:length(t)
            n1(j)=Y(2)'*Y(2);  

            Y=expm(-1i*H*dt)*Y;		%Y=(eye(2)-i*H*dt)*Y;
            Z=rand;
            if Y'*Y>Z
                Y=Y/(Y'*Y);
            else
                Y=[1 0]';
            end
        end
        nn1=nn1+n1;

        subplot(211),plot(t*1e9,rh1,'k',t*1e9,n1,'g','LineWidth',1.4)
        set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'FontSize',16,'LineWidth',1.4)
        ylabel('\fontsize{16}{\it\rho}_{11}')
%         title(['Run = ',num2str(jj)]);
        subplot(212),plot(t*1e9,rh1,'k',t*1e9,nn1/jj,'b','LineWidth',1.4)
        set(gca,'xLim',[0 max(t)*1e9],'yLim',[0 1],'FontSize',16,'LineWidth',1.4)
        xlabel('\fontsize{16}{\itt}   (ns)');
        ylabel('\fontsize{16}{\it\rho}_{11}')
        pause(.01);
    end
end

function M = MonteCarlo_Bloch(G,W,D)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bloch equation matrix.

    M=[	0		G		1i/2*W		-1i/2*W
        0		-G		-1i/2*W		1i/2*W
        1i/2*W	-1i/2*W	-1i*D-G/2	0
       -1i/2*W	1i/2*W	0			1i*D-G/2	];
end

function H = MonteCarlo_Hamilton(G,W,D)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hamilton matrix.

    H=[	0 	W/2
        W/2	-D-1i*G/1	];
end
