now one function solves the steady state susceptibility problem, thus we can define cellarray with initial conditions and then call cellfun to calculate them at once

1 files changed, 14 insertions, 30 deletions

diff --git a/liouville.m b/liouville.m
index 2f6a28e..fd1c840 100644
--- a/liouville.m
+++ b/liouville.m
@@ -35,7 +35,7 @@ detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps;
 [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
 rhoLiouville=zeros(N,1);
 
-% calculate E_field independent properties of athe atom
+% calculate E_field independent properties of the atom
 % to be used as sub matrix templates for Liouville operator matrix
 [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
 		Nlevels, ...
@@ -50,42 +50,26 @@ for detuning_p_cntr=1:N_detun_steps+1;
 	%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);
+	freq_index=freq2index(wp,modulation_freq);
 
-	% now we create Liouville indexes list
-	[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
+	atom_field_problem.L0m              = L0m;
+	atom_field_problem.polarizability_m = polarizability_m;
+	atom_field_problem.dipole_elements  = dipole_elements;
+	atom_field_problem.E_field          = E_field;
+	atom_field_problem.modulation_freq  = modulation_freq;
+	atom_field_problem.freq_index       = freq_index;
 
-	% Liouville operator matrix construction
-	L=Liouville_operator_matrix( 
-			N, 
-			L0m, polarizability_m,
-			E_field, 
-			modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
-			);
+	problems_cell_array{detuning_p_cntr}=atom_field_problem;
 
 
-	%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;
-
+	%kappa_p(detuning_p_cntr)=susceptibility_steady_state_at_freq( atom_field_problem);
+	detuning_freq(detuning_p_cntr)=detuning_p;
+endfor
 
-	%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
+% once we define all problems the main job is done here
+kappa_p=cellfun( @susceptibility_steady_state_at_freq, problems_cell_array);
 
-	kappa_p(detuning_p_cntr)=susceptibility(freq2index(wp,modulation_freq), rhoLiouville, dipole_elements, Nlevels, Nfreq);
-	%kappa_m(detuning_p_cntr)=susceptibility(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");
 %figure(3); plot(detuning_freq, real(kappa_m)); title("off resonant sideband dispersion");