major speed up by avoiding repetive calculations of the same submatrices

1 files changed, 74 insertions, 44 deletions

diff --git a/useful_functions.m b/useful_functions.m
index 38841a4..56989df 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -51,7 +51,48 @@ function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matr
 endfunction	
 
 
+% create (Nlevels*Nlevels)x*(Nlevels*Nlevels)
+% sub matrices of Liouville operator 
+% which repeat themselves for each modulation frequency
+% based on recipe from Eugeniy Mikhailov thesis
+function [L0m, polarization_m]=L0_and_polarization_submatrices( ...
+		rho_size, ...
+		H0, g_decay, g_dephasing, dipole_elements, ...
+		E_field, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
+		)
+	%-------------------------
+	useful_constants;
+	% note that L0 and decay parts depend only on combination of indexes
+	% jk,mn but repeats itself for every frequency
+	L0m=zeros(rho_size); % (NxN)x(NxN) matrix
+	decay_part_m=zeros(rho_size); % (NxN)x(NxN) matrix
+	% polarization matrix will be multiplied by field amplitude letter
+	% polarization is part of perturbation part of Hamiltonian
+	polarization_m=zeros(rho_size); % (NxN)x(NxN) matrix
+	for p=1:rho_size
+		% p= j*Nlevels+k 
+		% this might speed up stuff since less matrix passed back and force
+		j=rhoLiouville_r(p);
+		k=rhoLiouville_c(p);
+		for s=1:rho_size
+			% s= m*Nlevels+n 
+			m=rhoLiouville_r(s);
+			n=rhoLiouville_c(s);
+
+			% calculate unperturbed part (Hamiltonian without EM field)
+			L0m(p,s)=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m);
+			decay_part_m(p,s)=\
+				( 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) ;
+			polarization_m(p,s)= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) );
+		endfor
+	endfor
+	L0m=-im_one/hbar*L0m - decay_part_m; 
+endfunction
+
 % Liouville operator matrix construction
+% based on recipe from Eugeniy Mikhailov thesis
 function L=Liouville_operator_matrix( 
 		N, 
 		H0, g_decay, g_dephasing, dipole_elements,
@@ -64,71 +105,58 @@ function L=Liouville_operator_matrix(
 	Nfreq=length(modulation_freq);
 
 	% Lets be supper smart and speed up L matrix construction
-	% it has a lot of voids.
+	% since 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;
+	
+	[L0m, polarization_m]=L0_and_polarization_submatrices( ...
+		rho_size, ...
+		H0, g_decay, g_dephasing, dipole_elements, ...
+		E_field, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
+		);
+
+
+	% Liouville matrix operator has Nlevels*Nlevels blocks 
+	% which governed by the same modulation frequency
 	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(p0);
-			w2i=rhoLiouville_w(s0);
+			w1i=rhoLiouville_w(p0); % final
 			w_jk=modulation_freq(w1i);
+			w2i=rhoLiouville_w(s0); % initial
 			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 frequency in the list of available frequencies
 			% but since we not guaranteed to find it lets assign temporary 0 to Liouville matrix element
-			decay_part=0;
-			Lt=0;
 			for w3i=1:Nfreq
 				% at most we should have only one matching frequency
 				w_iner=modulation_freq(w3i);
 				if ((w_iner == w_l))
 					% 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
+					% lets do  L sub matrix filling
+
+					%such frequency exist in the list of modulation frequencies
+					if ((w_iner == 0))
+						% calculate unperturbed part (Hamiltonian without EM field)
+						L(p0:p0+rho_size-1,s0:s0+rho_size-1)=L0m;
+					else
+						% calculate perturbed part (Hamiltonian with EM field)
+						% in other word interactive part of Hamiltonian 
+						L(p0:p0+rho_size-1,s0:s0+rho_size-1)= ...
+							-im_one/hbar*polarization_m* E_field(w3i);
+					endif
+					% diagonal elements are self modulated
+					% due to rotating wave approximation
+					if ((p0 == s0))
+						L(p0:p0+rho_size-1,s0:s0+rho_size-1)+= -im_one*w_jk*eye(rho_size);
+					endif
 				endif
 			endfor
 		endfor
@@ -174,3 +202,5 @@ function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)
 		endfor
 	endfor
 endfunction
+
+% vim: ts=2:sw=2:fdm=indent