extra work moved to function, code for sweeping laser is added

2 files changed, 154 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));
 
diff --git a/useful_functions.m b/useful_functions.m
index 748a677..5522afb 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -28,3 +28,116 @@ function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq)
 	endfor
 endfunction
 
+% we unwrap density matrix to Liouville density vector and assign all possible
+% modulation frequencies as well
+% resulting vector should be Nlevels x Nlevels x length(modulation_freq)
+function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
+	N=Nfreq*Nlevels*Nlevels;
+	rhoLiouville_w=zeros(N,1);
+	rhoLiouville_r=zeros(N,1);
+	rhoLiouville_c=zeros(N,1);
+	i=0;
+	for w=1:Nfreq;
+		for r=1:Nlevels
+			for c=1:Nlevels
+				i+=1;
+				rhoLiouville(i)=0;
+				rhoLiouville_w(i)=w; % hold frequency modulation index
+				rhoLiouville_r(i)=r; % hold row value of rho_rc
+				rhoLiouville_c(i)=c; % hold column value of rho_rc
+			endfor
+		endfor
+	endfor
+endfunction	
+
+
+% Liouville operator matrix construction
+function L=Liouville_operator_matrix( 
+		N, 
+		H0, g_decay, g_dephasing, dipole_elements,
+		E_field, 
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
+		)
+	%-------------------------
+	useful_constants;
+	L=zeros(N); % NxN matrix
+
+	Nfreq=length(modulation_freq);
+
+	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))
+						% calculate unperturbed part (Hamiltonian without EM field)
+						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
+						% calculate perturbed part (Hamiltonian with EM field)
+						% in othe word interactive part of Hamiltonian 
+						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
+endfunction
+
+
+% 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
+function [rhoLiouville_dot, L]=constran_rho_and_match_L(
+		N, L,
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
+	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
+	rhoLiouville_dot= zeros(N,1);
+	% 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;
+endfunction
+