diff options
Diffstat (limited to 'liouville.m')
| -rw-r--r-- | liouville.m | 98 |
1 files changed, 26 insertions, 72 deletions
diff --git a/liouville.m b/liouville.m index baff45a..67dde6b 100644 --- a/liouville.m +++ b/liouville.m @@ -1,77 +1,21 @@ - -levels=1:3; - -% ----------- |1> -% / \ -% E_d / \ -% / \ E_p -% / \ -% -------- |3> \ -% \ -% ___________ |2> - - % some physical constants hbar=1; im_one=0+1i; -Nlevels=3; -w12=1e9; -w_hpf=6800; -w13=w12-w_hpf - - -% unperturbed Hamiltonian energy levels -levels_energy=[ w12, 0, w_hpf]; -levels_energy=levels_energy*hbar; -H0=zeros(Nlevels); -H0=diag(levels_energy); -%for i=1:Nlevels - %H0(i,i)=levels_energy(i); -%endfor - -% decay matrix g(i,j) correspnds to decay from i-->j -gamma=6; -gamma_13=.01 -g_decay=zeros(Nlevels); -g_decay(1,2)=gamma; %upper level decay -g_decay(1,3)=gamma; %upper level decay -g_decay(3,2)=gamma_13; % lower levels mixing -g_decay(2,3)=gamma_13; % lower levels mixing - -%defasing matris -g_deph=0; -g_dephasing=zeros(Nlevels); -g_dephasing(1,2)=g_deph; -g_dephasing(2,1)=g_dephasing(1,2); -g_dephasing(1,3)=g_deph; -g_dephasing(3,1)=g_dephasing(1,3); - - - -% dipole matrix -dipole_elements=zeros(Nlevels); -dipole_elements(1,2)=1; -dipole_elements(2,1)=dipole_elements(1,2); -dipole_elements(3,1)=1; -dipole_elements(1,3)=dipole_elements(3,1); - - -%EM field definition -Ep=100.0; -Epc=conj(Ep); -Ed=0; -Edc=conj(Ed); -wd=w13; -wp=w12; -modulation_freq=[0, wp, wd, -wp, -wd, wp-wd, wd-wp]; -E_field =[0, Ep, Ed, Epc, Edc, 0, 0 ]; +%three_levels; +two_levels; + Nfreq=length(modulation_freq); % now we create Liouville indexes list % we unwrap density matrix and assign all posible % frequencies as well % resulting vector should be Nlevels x Nlevels x length(modulation_freq) +N=length(modulation_freq)*Nlevels*Nlevels; +rhoLiouville=zeros(N,1); +rhoLiouville_w=rhoLiouville; +rhoLiouville_r=rhoLiouville; +rhoLiouville_c=rhoLiouville; i=0; for w=1:length(modulation_freq) for r=1:Nlevels @@ -85,9 +29,7 @@ for w=1:length(modulation_freq) endfor endfor -\ -N=length(rhoLiouville); % Liouville operator matrix L=zeros(N); % NxN matrix Li=zeros(N); % NxN Liouville interactive @@ -101,7 +43,11 @@ function ret=decay_total(g_decay,i) endfunction function ret=kron_delta(i,j) - ret=((i==j)); + if ((i==j)) + ret=1; + else + ret=0; + endif endfunction for p=1:N @@ -130,17 +76,21 @@ for p=1:N ( decay_total(g_decay,k)/2 + decay_total(g_decay,j)/2 + g_dephasing(j,k) )* kron_delta(j,m)*kron_delta(k,n) \ - kron_delta(m,n)*kron_delta(j,k)*g_decay(m,j) ; decay_part=decay_part*hbar/im_one; - L0=decay_part; + L0+=decay_part; Lt=L0; else - Li= ( dipole_elements(j,m)*kron_delta(k,n)+dipole_elements(n,k)*kron_delta(j,m) )*E_field(w3i); + Li= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) )*E_field(w3i); Lt=Li; endif - Lt=-im_one/hbar*Lt*kron_delta(w_jk-w_iner,w_mn); % above if should be done only if kron_delta is not zero - Lt+=-im_one*w_jk*kron_delta(w_iner,w_jk); - L(p,s)=Lt; + %Lt=-im_one/hbar*Lt*kron_delta(w_jk-w_iner,w_mn); % above if should be done only if kron_delta is not zero + % no need for above kron_delta since the same conditon checked in the outer if statement + Lt=-im_one/hbar*Lt; endif endfor + if ((p == s)) + Lt+=-im_one*w_jk; + endif + L(p,s)=Lt; endfor endfor @@ -158,6 +108,8 @@ for i=1:N w=modulation_freq(w2i); if ((w==0) & (m==n)) L(1,i)=1; + else + L(1,i)=0; endif endfor @@ -184,6 +136,8 @@ endfunction rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq) +rho_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq) +rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq) %rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq) %L*rhoLiouville |