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

1 files changed, 44 insertions, 35 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