1 files changed, 41 insertions, 99 deletions
diff --git a/liouville.m b/liouville.m
index bc9d524..b232a56 100644
--- a/liouville.m
+++ b/liouville.m
@@ -3,9 +3,7 @@
 useful_functions;
 
 % some physical constants
-hbar=1;
-im_one=0+1i;
-
+useful_constants;
 
 % load atom energy levels and decay description
 three_levels;
@@ -16,113 +14,57 @@ field_description;
 
 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
-		for c=1:Nlevels
-			i+=1;
-			rhoLiouville(i)=0;
-			rhoLiouville_w(i)=w;
-			rhoLiouville_r(i)=r;
-			rhoLiouville_c(i)=c;
-		endfor
-	endfor
-endfor
 
 
-% Liouville operator matrix 
-L=zeros(N); % NxN matrix
-Li=zeros(N); % NxN Liouville interactive
-L0=zeros(N); % NxN Liouville from unperturbed hamiltonian
-
-for p=1:N
-	for s=1:N
-		j=rhoLiouville_r(p);
-		k=rhoLiouville_c(p);
-		m=rhoLiouville_r(s);
-		n=rhoLiouville_c(s);
-		% we garanted to know frequency of final and initial rhoLiouville
-		w1i=rhoLiouville_w(p);
-		w2i=rhoLiouville_w(s);
-		w_jk=modulation_freq(w1i);
-		w_mn=modulation_freq(w2i);
-		% thus we know L matrix element frequency which we need to match
-		w_l=w_jk-w_mn;
-		% lets search this wrequency in the list of available frequencyes
-		% but since we not garanteed to find it lets assign temporary 0 to Liouville matrix element
-		L(p,s)=0;
-		decay_part=0;
-		Lt=0;
-		for w3i=1:Nfreq
-			w_iner=modulation_freq(w3i);
-			decay_part=0;
-			if ((w_iner == w_l))
-				%such frequency exist in the list of modulation frequencies
-				if ((w_iner == 0))
-					L0=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m);
-					decay_part=\
-						( 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) ;
-					Lt=L0;
-				else
-					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
-				% no need for above kron_delta since the same conditon checked in the outer if statement
-				Lt=-im_one/hbar*Lt - decay_part; 
-			endif
-		endfor
-		if ((p == s))
-			Lt+=-im_one*w_jk;
-		endif
-		L(p,s)=Lt;
-	endfor
-endfor
+%tune probe frequency
+detuning_p=0;
+N_detun_steps=15;
+detuning_p_min=-15;
+detuning_p_max=-detuning_p_min;
+detuning_freq=zeros(1,N_detun_steps+1);
+kappa_p      =zeros(1,N_detun_steps+1);
+detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps;
+for detuning_p_cntr=1:N_detun_steps+1;
 
+wp0=w12;
+detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1);
+wp=wp0+detuning_p;
+modulation_freq=[0,  wp, wd, -wp, -wd, wp-wd, wd-wp];
 
-% now generally rhoL_dot=0=L*rhoL has infinite number of solutions
-% since we always can resclale rho vector with arbitrary constant
-% lets constrain our density matrix with some physical meaning
-% sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1
-% we will replace first row of Liouville operator with this condition
-% thus rhoLiouville_dot(1)=1
-for i=1:N
-	w2i=rhoLiouville_w(i);
-	m=rhoLiouville_r(i);
-	n=rhoLiouville_c(i);
-	w=modulation_freq(w2i);
-	if ((w==0) & (m==n))
-		L(1,i)=1;
-	else
-		L(1,i)=0;
-	endif
-endfor
+% now we create Liouville indexes list
+[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
+rhoLiouville=zeros(N,1);
+
+% Liouville operator matrix construction
+L=Liouville_operator_matrix( 
+		N, 
+		H0, g_decay, g_dephasing, dipole_elements,
+		E_field, 
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
+		);
 
-rhoLiouville_dot=rhoLiouville*0;
-% sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1
-% we will replace first row of Liouville operator with this condition
-% thus rhoLiouville_dot(1)=1
-rhoLiouville_dot(1)=1;
 
+%use the fact that sum(rho_ii)=1 to constrain solution
+[rhoLiouville_dot, L]=constran_rho_and_match_L(
+		N, L,
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);
 
 %solving for density matrix vector
 rhoLiouville=L\rhoLiouville_dot;
 
 
-rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq)
-rho_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq)
-rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq)
+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
+kappa_p(detuning_p_cntr)=sum(sum(rho_1));
+detuning_freq(detuning_p_cntr)=detuning_p;
+
+%kappa_p_re=real(kappa_p);
+%kappa_p_im=imag(kappa_p);
+
+endfor
+figure(1); plot(detuning_freq, real(kappa_p));
+figure(2); plot(detuning_freq, imag(kappa_p));