use of L operator matrix property to obtain speed up of factor of 2

1 files changed, 55 insertions, 32 deletions

diff --git a/useful_functions.m b/useful_functions.m
index 1765dd9..38841a4 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -59,19 +59,24 @@ function L=Liouville_operator_matrix(
 		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
 		)
 	%-------------------------
-	L=zeros(N); % NxN matrix
 	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);
+	% Lets be supper smart and speed up L matrix construction
+	% it has a lot of voids.
+	% By creation of rhoLiouville we know that there are
+	% consequent chunks of rho_ij modulated with same frequency
+	% this means that rhoLiouville is split in to Nfreq chunks 
+	% with length Nlevels*Nlevels=N/Nfreq
+	rho_size=N/Nfreq;
+	for p_freq_cntr=1:Nfreq
+		for s_freq_cntr=1:Nfreq
+			p0=1+(p_freq_cntr-1)*rho_size;
+			s0=1+(s_freq_cntr-1)*rho_size;
 			% we guaranteed to know frequency of final and initial rhoLiouville
-			w1i=rhoLiouville_w(p);
-			w2i=rhoLiouville_w(s);
+			w1i=rhoLiouville_w(p0);
+			w2i=rhoLiouville_w(s0);
 			w_jk=modulation_freq(w1i);
 			w_mn=modulation_freq(w2i);
 			% thus we know L matrix element frequency which we need to match
@@ -81,33 +86,51 @@ function L=Liouville_operator_matrix(
 			decay_part=0;
 			Lt=0;
 			for w3i=1:Nfreq
+				% at most we should have only one matching frequency
 				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; 
+					% yey, requested modulation frequency exist
+					% lets do  L submatrix filling
+					for p_cntr=0:(rho_size-1)
+						p=p0+p_cntr;
+						for s_cntr=0:(rho_size-1)
+							s=s0+s_cntr;
+
+							decay_part=0;
+							Lt=0;
+							
+							j=rhoLiouville_r(p);
+							k=rhoLiouville_c(p);
+							m=rhoLiouville_r(s);
+							n=rhoLiouville_c(s);
+
+							%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; 
+							% diagonal elements are self modulated
+							if ((p == s))
+								Lt+=-im_one*w_jk;
+							endif
+							Lps=Lt;
+							L(p,s)=Lps;
+						endfor
+					endfor
 				endif
 			endfor
-			if ((p == s))
-				Lt+=-im_one*w_jk;
-			endif
-			Lps=Lt;
-			L(p,s)=Lps;
 		endfor
 	endfor
 endfunction