speed up by another factor of 4 by calculatiing field independent atom property matrices

2 files changed, 50 insertions, 48 deletions

diff --git a/liouville.m b/liouville.m
index 204c1a3..4c78557 100644
--- a/liouville.m
+++ b/liouville.m
@@ -29,53 +29,60 @@ detuning_freq=zeros(1,N_detun_steps+1);
 kappa_p      =zeros(1,N_detun_steps+1);
 kappa_m      =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;
-wm=wd-(wp-wd);
-%modulation_freq=[0,  wp, wd, wm,   -wp, -wd, -wm,  wp-wd, wd-wp];
-%E_field        =[0,  Ep, Ed, Em,   Epc, Edc, Emc,  0,     0    ];
-modulation_freq=[0,  wp, wd, -wp, -wd,  wp-wd, wd-wp];
-E_field        =[0,  Ep, Ed,  Epc, Edc, 0,     0    ];
-Nfreq=length(modulation_freq);
 
-% 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
+% calculate E_field independent properties of athe atom
+% to be used as sub matrix templates for Liouville operator matrix
+[L0m, polarizability_m]=L0_and_polarization_submatrices( ...
+		Nlevels*Nlevels, ...
+		H0, g_decay, g_dephasing, dipole_elements, ...
+		E_field, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
 		);
+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;
+	wm=wd-(wp-wd);
+	%modulation_freq=[0,  wp, wd, wm,   -wp, -wd, -wm,  wp-wd, wd-wp];
+	%E_field        =[0,  Ep, Ed, Em,   Epc, Edc, Emc,  0,     0    ];
+	modulation_freq=[0,  wp, wd, -wp, -wd,  wp-wd, wd-wp];
+	E_field        =[0,  Ep, Ed,  Epc, Edc, 0,     0    ];
+	Nfreq=length(modulation_freq);
+
+	% now we create Liouville indexes list
+	[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
 
+	% Liouville operator matrix construction
+	L=Liouville_operator_matrix( 
+			N, 
+			L0m, polarizability_m,
+			E_field, 
+			modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
+			);
 
-%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);
 
+	%use the fact that sum(rho_ii)=1 to constrain solution
+	[rhoLiouville_dot, L]=constrain_rho_and_match_L(
+			N, L,
+			modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);
 
-%solving for density matrix vector
-rhoLiouville=L\rhoLiouville_dot;
 
+	%solving for density matrix vector
+	rhoLiouville=L\rhoLiouville_dot;
 
-%rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq); % 0 frequency, 
-%rho_p=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq); % probe frequency
-%rho_d=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq); % drive frequency
-%rho_m=rhoOfFreq(rhoLiouville, 4, Nlevels, Nfreq); % opposite sideband frequency
 
-kappa_p(detuning_p_cntr)=sucseptibility(2, rhoLiouville, dipole_elements, Nlevels, Nfreq);
-%kappa_m(detuning_p_cntr)=sucseptibility(4, rhoLiouville, dipole_elements, Nlevels, Nfreq);
-detuning_freq(detuning_p_cntr)=detuning_p;
+	%rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq); % 0 frequency, 
+	%rho_p=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq); % probe frequency
+	%rho_d=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq); % drive frequency
+	%rho_m=rhoOfFreq(rhoLiouville, 4, Nlevels, Nfreq); % opposite sideband frequency
 
-%kappa_p_re=real(kappa_p);
-%kappa_p_im=imag(kappa_p);
+	kappa_p(detuning_p_cntr)=sucseptibility(2, rhoLiouville, dipole_elements, Nlevels, Nfreq);
+	%kappa_m(detuning_p_cntr)=sucseptibility(4, rhoLiouville, dipole_elements, Nlevels, Nfreq);
+	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)); title("probe dispersion");
 figure(2); plot(detuning_freq, imag(kappa_p)); title("probe absorption");
@@ -83,3 +90,5 @@ figure(2); plot(detuning_freq, imag(kappa_p)); title("probe absorption");
 %figure(4); plot(detuning_freq, imag(kappa_m)); title("off resonant absorption");
 
 elapsed_time = etime (clock (), t0)
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/useful_functions.m b/useful_functions.m
index 56989df..a1135f7 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -55,7 +55,7 @@ endfunction
 % 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( ...
+function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
 		rho_size, ...
 		H0, g_decay, g_dephasing, dipole_elements, ...
 		E_field, ...
@@ -69,7 +69,7 @@ function [L0m, polarization_m]=L0_and_polarization_submatrices( ...
 	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
+	polarizability_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
@@ -85,7 +85,7 @@ function [L0m, polarization_m]=L0_and_polarization_submatrices( ...
 			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) );
+			polarizability_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; 
@@ -95,7 +95,7 @@ endfunction
 % based on recipe from Eugeniy Mikhailov thesis
 function L=Liouville_operator_matrix( 
 		N, 
-		H0, g_decay, g_dephasing, dipole_elements,
+		L0m, polarizability_m,
 		E_field, 
 		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
 		)
@@ -112,13 +112,6 @@ function L=Liouville_operator_matrix(
 	% 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
@@ -150,7 +143,7 @@ function L=Liouville_operator_matrix(
 						% 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);
+							-im_one/hbar*polarizability_m* E_field(w3i);
 					endif
 					% diagonal elements are self modulated
 					% due to rotating wave approximation
@@ -170,7 +163,7 @@ endfunction
 % 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(
+function [rhoLiouville_dot, L]=constrain_rho_and_match_L(
 		N, L,
 		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
 	for i=1:N