TAG v1.1 with major speedups

1 files changed, 22 insertions, 20 deletions

diff --git a/useful_functions.m b/useful_functions.m
index 606fe49..c0fdd7e 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -1,15 +1,15 @@
 1;
 
-% calculate total decay for particular level taking in account all branches
 function ret=decay_total(g_decay,i)
+% calculate total decay for particular level taking in account all branches
 	ret=0;
 	for k=1:size(g_decay)(1)
 		ret=ret+g_decay(i,k);
 	endfor
 endfunction
 
-% kroneker delta symbol
 function ret=kron_delta(i,j)
+% kroneker delta symbol
 	if ((i==j))
 		ret=1;
 	else
@@ -17,9 +17,9 @@ function ret=kron_delta(i,j)
 	endif
 endfunction
 
+function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq)
 % this function create from Liouville density vector 
 % the density matrix with given modulation frequency
-function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq)
 	rho=zeros(Nlevels);
 	for r=1:Nlevels
 		for c=1:Nlevels
@@ -28,10 +28,10 @@ function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq)
 	endfor
 endfunction
 
-% we unwrap density matrix to Liouville density vector and assign all possible
+function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
+% 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);
@@ -51,21 +51,24 @@ function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matr
 endfunction	
 
 
+function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
+		Nlevels, ...
+		H0, g_decay, g_dephasing, dipole_elements ...
+		)
 % 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, polarizability_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;
+	rho_size=Nlevels*Nlevels;
+
+	% now we create Liouville indexes list
+	[Ndummy, rhoLiouville_w_notused, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,1);
+
 	% 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
+	L0m=zeros(rho_size); % (Nlevels^2)x(Nlevels^2) 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
@@ -95,14 +98,14 @@ function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
 	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, ... 
 		L0m, polarizability_m, ...
 		E_field, ...
 		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
 		)
+% Liouville operator matrix construction
+% based on recipe from Eugeniy Mikhailov thesis
 	%-------------------------
 	useful_constants;
 	L=zeros(N); % NxN matrix
@@ -161,15 +164,15 @@ function L=Liouville_operator_matrix( ...
 endfunction
 
 
+function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ...
+		N, L, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
 % 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]=constrain_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);
@@ -188,9 +191,8 @@ function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ...
 	rhoLiouville_dot(1)=1;
 endfunction
 
-% calculate sucseptibility for the field at given frequency index
-function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)
-
+function kappa=susceptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)
+% calculate susceptibility for the field at given frequency index
 	rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq);
 	kappa=0;
 	for i=1:Nlevels