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

1 files changed, 113 insertions, 0 deletions

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
+