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, 43 insertions, 3 deletions

diff --git a/useful_functions.m b/useful_functions.m
index dda3d1d..e0c9e8d 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -17,7 +17,7 @@ function ret=kron_delta(i,j)
 	endif
 endfunction
 
-function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq)
+function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels)
 % this function create from Liouville density vector 
 % the density matrix with given modulation frequency
 	rho=zeros(Nlevels);
@@ -191,9 +191,10 @@ function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ...
 	rhoLiouville_dot(1)=1;
 endfunction
 
-function kappa=susceptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)
+function kappa=susceptibility(wi, rhoLiouville, dipole_elements)
 % calculate susceptibility for the field at given frequency index
-	rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq);
+	Nlevels=( size(dipole_elements)(1) );
+	rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
 	kappa=0;
 	for i=1:Nlevels
 		for j=1:Nlevels
@@ -207,4 +208,43 @@ function index=freq2index(freq, modulation_freq)
 	index=[1:length(modulation_freq)](modulation_freq==freq);
 endfunction
 
+function rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq)
+% calculates rhoLiouville vector assuming steady state situation and normalization of rho_ii to 1
+	Nlevels=sqrt( size(L0m)(1) );
+	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]=constrain_rho_and_match_L(
+			N, L,
+			modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);
+
+	%solving for density matrix vector
+	rhoLiouville=L\rhoLiouville_dot;
+endfunction
+
+function xi=susceptibility_steady_state_at_freq( atom_field_problem)
+	% find steady state susceptibility at particular modulation frequency element
+	% at given E_field 
+	L0m                = atom_field_problem.L0m              ; 
+	polarizability_m   = atom_field_problem.polarizability_m ; 
+	dipole_elements    = atom_field_problem.dipole_elements  ; 
+	E_field            = atom_field_problem.E_field          ; 
+	modulation_freq    = atom_field_problem.modulation_freq  ; 
+	freq_index         = atom_field_problem.freq_index       ; 
+
+	rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
+	xi=susceptibility(freq_index, rhoLiouville, dipole_elements);
+endfunction		
+
 % vim: ts=2:sw=2:fdm=indent