Ressonância entre um pente de frequências e um átomo de dois níveis - programa no Matlab

Voltar

Programa em C++ para o estudo da ressonância da interação entre um pente de frequências e um átomo de dois níveis.

Função bloch - Equações de Bloch:

function [pop] = bloch(a1,a2,a12,b12,I,Omega,alpha,gama,a10,q22,dd)
if I==1 
    pop =  2*Omega*(b12*cos(alpha)-a12*sin(alpha))  +q22*a2 -(a1-a10)*gama;
elseif I==2 
    pop = -2*Omega*(b12*cos(alpha)-a12*sin(alpha))  -q22*a2       -a2*gama;
elseif I==3 
    pop = -dd*b12-Omega*(a2-a1)*sin(alpha)                    -0.5*q22*a12  -a12*gama;
elseif I==4 
    pop =  dd*a12+Omega*(a2-a1)*cos(alpha)                    -0.5*q22*b12  -b12*gama;
end

Arquivo do programa principal:

% 23-03-2012

clear;clc;

    gama=0.0025*1e9*0;
    PassoRKfento=10;
    q22=(2*pi)*5e6;
    q12=0.5*q22;
    fr=100e6;
    Tp=100e-15;            
    Pulsos=100;            
    Omega=0.01*q22/(fr*Tp)
    a10=1;                 
    a20=0;                 
    w2=(2*pi)*400e12+(2*pi)*20e6*0;
    wL=(2*pi)*400e12+(2*pi)*40e6*0;
    phi=(2*pi)*0.4*0;
               
    a(1)=a10;
    a(2)=a20;
    a(3)=0;
    a(4)=0;   
    t=0;N=-1;
    dd=w2-wL;
    Tr=1/fr;
    
    for I = 1:2*Pulsos+1
        if I/2 - ceil(I/2) == 0
            N = N + 1
            g = PassoRKfento;
            h = (1/g)*Tp;
            
            alpha = -N*wL*Tr + N*phi;
            for k = 1:g
                t = t + h;
                
                for j = 1:4
                    k1(j)=bloch(a(1),a(2),a(3),a(4),j,Omega,alpha,gama,a10,q22,dd);
                end
                                                            
                for j = 1:4
                    k2(j)=bloch(a(1)+k1(1)*h/2,a(2)+k1(2)*h/2,a(3)+k1(3)*h/2,a(4)+k1(4)*h/2,j,Omega,alpha,gama,a10,q22,dd);
                end
                   
                for j = 1:4
                    k3(j)=bloch(a(1)+k2(1)*h/2,a(2)+k2(2)*h/2,a(3)+k2(3)*h/2,a(4)+k2(4)*h/2,j,Omega,alpha,gama,a10,q22,dd);
                end
                   
                for j = 1:4
                    k4(j)=bloch(a(1)+k3(1)*h,a(2)+k3(2)*h,a(3)+k3(3)*h,a(4)+k3(4)*h,j,Omega,alpha,gama,a10,q22,dd);
                end
 
                for j = 1:4
                    b(j)=a(j)+h*(k1(j)/6+k2(j)/3+k3(j)/3+k4(j)/6); 
                end
                   
                for m = 1:4
                    a(m)=b(m);
                end
            
            end
            
        elseif I/2 - ceil(I/2) == -0.5
            t=t+(Tr-Tp);
            x=dd*(Tr-Tp);
            y=(Tr-Tp)*q22;
                                          
            b(1)=a(1)+a(2)*(1-exp(-y));
            b(2)=a(2)*exp(-y);
            b(3)=(a(3)*cos(x)-a(4)*sin(x))*exp(-0.5*y);
            b(4)=(a(3)*sin(x)+a(4)*cos(x))*exp(-0.5*y);
               
            for m = 1:4
                a(m)=b(m);
            end
            
        end
        
        w = b(1) + b(2);
        A(I,1) = t;
        A(I,2) = b(2);
        csvwrite('C:\Users\Marco Polo\Desktop\dados.dat', A);
        
    end