seems too start outputing reasonable values at least for 2 levels system

3 files changed, 147 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
diff --git a/three_levels.m b/three_levels.m
new file mode 100644
index 0000000..beb70b3
--- /dev/null
+++ b/three_levels.m
@@ -0,0 +1,67 @@
+levels=1:3;
+
+%      -----------  |1>
+%         /     \
+%    E_d /       \ 
+%       /         \ E_p
+%      /           \
+%    -------- |3>   \
+%                    \
+%               ___________ |2>
+
+
+
+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_23=.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_23; % lower levels mixing
+g_decay(2,3)=gamma_23; % 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=10.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    ];
+modulation_freq=[0,  wp, -wp];
+E_field        =[0,  Ep,  Epc];
+
diff --git a/two_levels.m b/two_levels.m
new file mode 100644
index 0000000..18860cb
--- /dev/null
+++ b/two_levels.m
@@ -0,0 +1,54 @@
+1;
+
+%      -----------  |1>
+%               \
+%                \ 
+%                 \ E_p
+%                  \
+%                   \
+%                    \
+%               ___________ |2>
+
+
+
+Nlevels=2;
+w12=1e4;
+w_hpf=68;
+w13=w12-w_hpf;
+
+
+% unperturbed Hamiltonian energy levels
+levels_energy=[ w12, 0];
+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_23=.01;
+g_decay=zeros(Nlevels);
+g_decay(1,2)=gamma; %upper level decay
+
+%defasing matris
+g_deph=0;
+g_dephasing=zeros(Nlevels);
+g_dephasing(1,2)=g_deph;
+g_dephasing(2,1)=g_dephasing(1,2);
+
+
+% dipole matrix
+dipole_elements=zeros(Nlevels);
+dipole_elements(1,2)=1;
+dipole_elements(2,1)=dipole_elements(1,2);
+
+
+%EM field definition
+Ep=1.0;
+Epc=conj(Ep);
+wp=w12;
+modulation_freq=[0,  wp, -wp];
+E_field        =[0,  Ep,  Epc];
+