now we treat atom as polarization aware one

1 files changed, 21 insertions, 7 deletions

diff --git a/useful_functions.m b/useful_functions.m
index 74dd9a5..4b239a5 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -72,7 +72,9 @@ function [L0m, polarizability_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
-	polarizability_m=zeros(rho_size); % (NxN)x(NxN) matrix
+	polarizability_m.linear = zeros(rho_size); % (NxN)x(NxN) matrix
+	polarizability_m.left   = zeros(rho_size); % (NxN)x(NxN) matrix
+	polarizability_m.right  = 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
@@ -92,7 +94,9 @@ function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
 					+ g_dephasing(j,k) ...
 				)* kron_delta(j,m)*kron_delta(k,n) ...
 				- kron_delta(m,n)*kron_delta(j,k)*g_decay(m,j) ;
-			polarizability_m(p,s)= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) );
+			polarizability_m.linear(p,s)= ( dipole_elements.linear(j,m)*kron_delta(k,n)-dipole_elements.linear(n,k)*kron_delta(j,m) );
+			polarizability_m.left(p,s)= ( dipole_elements.left(j,m)*kron_delta(k,n)-dipole_elements.left(n,k)*kron_delta(j,m) );
+			polarizability_m.right(p,s)= ( dipole_elements.right(j,m)*kron_delta(k,n)-dipole_elements.right(n,k)*kron_delta(j,m) );
 		endfor
 	endfor
 	L0m=-im_one/hbar*L0m - decay_part_m; 
@@ -150,7 +154,10 @@ 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*polarizability_m* E_field(w3i);
+							-im_one/hbar*polarizability_m.linear * E_field.linear(w3i) ...
+							-im_one/hbar*polarizability_m.left   * E_field.left(w3i)   ...
+							-im_one/hbar*polarizability_m.right  * E_field.right(w3i)  ...
+							;
 					endif
 					% diagonal elements are self modulated
 					% due to rotating wave approximation
@@ -193,12 +200,16 @@ endfunction
 
 function kappa=susceptibility(wi, rhoLiouville, dipole_elements)
 % calculate susceptibility for the field at given frequency index
-	Nlevels=( size(dipole_elements)(1) );
+	Nlevels=( size(dipole_elements.linear)(1) );
 	rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
-	kappa=0;
+	kappa.linear=0;
+	kappa.left=0;
+	kappa.right=0;
 	for i=1:Nlevels
 		for j=1:Nlevels
-			kappa+=dipole_elements(j,i)*rho(i,j);
+			kappa.linear += dipole_elements.linear(j,i) * rho(i,j);
+			kappa.left   += dipole_elements.left(j,i)   * rho(i,j);
+			kappa.right  += dipole_elements.right(j,i)  * rho(i,j);
 		endfor
 	endfor
 endfunction
@@ -233,7 +244,7 @@ function rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field,
 	rhoLiouville=L\rhoLiouville_dot;
 endfunction
 
-function xi=susceptibility_steady_state_at_freq( atom_field_problem)
+function [xi_linear, xi_left, xi_right]=susceptibility_steady_state_at_freq( atom_field_problem)
 	% find steady state susceptibility at particular modulation frequency element
 	% at given E_field 
 	global atom_properties;
@@ -247,6 +258,9 @@ function xi=susceptibility_steady_state_at_freq( atom_field_problem)
 
 	rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
 	xi=susceptibility(freq_index, rhoLiouville, dipole_elements);
+	xi_linear = xi.linear;
+	xi_right  = xi.right;
+	xi_left   = xi.left;
 endfunction		
 
 % vim: ts=2:sw=2:fdm=indent