copy psr to faraday

1 files changed, 337 insertions, 0 deletions

diff --git a/faraday/useful_functions.m b/faraday/useful_functions.m
new file mode 100644
index 0000000..eb3e880
--- /dev/null
+++ b/faraday/useful_functions.m
@@ -0,0 +1,337 @@
+1;
+
+function ret=decay_total(g_decay,i)
+% calculate total decay for particular level taking in account all branches
+	ret=sum(g_decay(i,:));
+endfunction
+
+function ret=kron_delta(i,j)
+% kroneker delta symbol
+	if ((i==j))
+		ret=1;
+	else
+		ret=0;
+	endif
+endfunction
+
+function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels)
+% this function create from Liouville density vector 
+% the density matrix with given modulation frequency
+	rho=zeros(Nlevels);
+	rho(:)=rhoLiouville((freqIndex-1)*Nlevels^2+1:(freqIndex)*Nlevels^2);
+	rho=rho.';
+endfunction
+
+function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
+% unwrap density matrix to Liouville density vector and assign all possible
+% modulation frequencies as well
+% resulting vector should be Nlevels x Nlevels x length(modulation_freq)
+	N  = Nfreq*Nlevels*Nlevels;
+	rho_size = Nlevels*Nlevels;
+	rhoLiouville_w=zeros(N,1);
+	rhoLiouville_r=zeros(N,1);
+	rhoLiouville_c=zeros(N,1);
+	
+	w=1:Nfreq;
+	w_tmplate=(repmat(w,rho_size,1))(:);
+	rhoLiouville_w=w_tmplate;
+	r=1:Nlevels;
+	r_tmplate=(repmat(r,Nlevels,1))(:);
+	rhoLiouville_r=(repmat(r_tmplate,Nfreq,1))(:)';
+	c=(1:Nlevels)';% hold column value of rho_rc
+	rhoLiouville_c=repmat(c,Nfreq*Nlevels,1);
+endfunction	
+
+
+function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
+		Nlevels, ...
+		H0, g_decay, g_dephasing, dipole_elements ...
+		)
+% create (Nlevels*Nlevels)x*(Nlevels*Nlevels)
+% sub matrices of Liouville operator 
+% which repeat themselves for each modulation frequency
+% based on recipe from Eugeniy Mikhailov thesis
+	%-------------------------
+	useful_constants;
+	rho_size=Nlevels*Nlevels;
+
+	% now we create Liouville indexes list
+	[Ndummy, rhoLiouville_w_notused, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,1);
+
+	kron_delta_m=eye(Nlevels);
+	% note that L0 and decay parts depend only on combination of indexes
+	% jk,mn but repeats itself for every frequency
+	L0m=zeros(rho_size); % (Nlevels^2)x(Nlevels^2) matrix
+	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.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
+		j=rhoLiouville_r(p);
+		k=rhoLiouville_c(p);
+		for s=1:rho_size
+			% s= m*Nlevels+n 
+			m=rhoLiouville_r(s);
+			n=rhoLiouville_c(s);
+
+			% calculate unperturbed part (Hamiltonian without EM field)
+			L0m(p,s)=H0(j,m)*kron_delta_m(k,n)-H0(n,k)*kron_delta_m(j,m);
+			decay_part_m(p,s)= ...
+				( ...
+					  decay_total(g_decay,k)/2 ...
+					+ decay_total(g_decay,j)/2 ...
+					+ g_dephasing(j,k) ...
+				)* kron_delta_m(j,m)*kron_delta_m(k,n) ...
+				- kron_delta_m(m,n)*kron_delta_m(j,k)*g_decay(m,j) ;
+			polarizability_m.linear(p,s)= ( dipole_elements.linear(j,m)*kron_delta_m(k,n)-dipole_elements.linear(n,k)*kron_delta_m(j,m) );
+			polarizability_m.left(p,s)= ( dipole_elements.left(j,m)*kron_delta_m(k,n)-dipole_elements.left(n,k)*kron_delta_m(j,m) );
+			polarizability_m.right(p,s)= ( dipole_elements.right(j,m)*kron_delta_m(k,n)-dipole_elements.right(n,k)*kron_delta_m(j,m) );
+		endfor
+	endfor
+	L0m=-im_one/hbar*L0m - decay_part_m; 
+endfunction
+
+function L=Liouville_operator_matrix( ...
+		N, ... 
+		L0m, polarizability_m, ...
+		E_field, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
+		)
+% Liouville operator matrix construction
+% based on recipe from Eugeniy Mikhailov thesis
+	%-------------------------
+	useful_constants;
+	L=zeros(N); % NxN matrix
+	Nfreq=length(modulation_freq);
+
+	% Lets be supper smart and speed up L matrix construction
+	% since it has a lot of voids.
+	% By creation of rhoLiouville we know that there are
+	% consequent chunks of rho_ij modulated with same frequency
+	% this means that rhoLiouville is split in to Nfreq chunks 
+	% with length Nlevels*Nlevels=N/Nfreq
+	rho_size=N/Nfreq;
+	
+	% creating building blocks of L by rho_size * rho_size
+	for w3i=1:Nfreq
+		w_iner=modulation_freq(w3i);
+		if ((w_iner == 0))
+			% calculate unperturbed part (Hamiltonian without EM field)
+			L_sub{w3i}=L0m;
+		else
+			% calculate perturbed part (Hamiltonian with EM field)
+			% in other word interactive part of Hamiltonian 
+			L_sub{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
+	endfor
+
+	% Liouville matrix operator has Nlevels*Nlevels blocks 
+	% which governed by the same modulation frequency
+	for p_freq_cntr=1:Nfreq
+		p0=1+(p_freq_cntr-1)*rho_size;
+		% we guaranteed to know frequency of final and initial rhoLiouville
+		w1i=rhoLiouville_w(p0); % final
+		w_jk=modulation_freq(w1i);
+		for s_freq_cntr=1:Nfreq
+			s0=1+(s_freq_cntr-1)*rho_size;
+			w2i=rhoLiouville_w(s0); % initial
+			w_mn=modulation_freq(w2i);
+			% thus we know L matrix element frequency which we need to match
+			w_l=w_jk-w_mn;
+			% lets search this frequency in the list of available frequencies
+			% but since we not guaranteed to find it lets assign temporary 0 to Liouville matrix element
+			w3i=(w_l == modulation_freq);
+			if (any(w3i))	
+				% yey, requested modulation frequency exist
+				% lets do  L sub matrix filling
+				% at most we should have only one matching frequency
+				w_iner=modulation_freq(w3i);
+				L(p0:p0+rho_size-1,s0:s0+rho_size-1) = L_sub{w3i};
+
+			endif
+		endfor
+		% diagonal elements are self modulated
+		% due to rotating wave approximation
+		L(p0:p0+rho_size-1,p0:p0+rho_size-1)+= -im_one*w_jk*eye(rho_size);
+	endfor
+endfunction
+
+
+function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ...
+		N, L, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
+% now generally rhoL_dot=0=L*rhoL has infinite number of solutions
+% since we always can resclale rho vector with arbitrary constant
+% lets constrain our density matrix with some physical meaning
+% 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
+	for i=1:N
+		w2i=rhoLiouville_w(i);
+		m=rhoLiouville_r(i);
+		n=rhoLiouville_c(i);
+		w=modulation_freq(w2i);
+		if ((w==0) & (m==n))
+			L(1,i)=1;
+		else
+			L(1,i)=0;
+		endif
+	endfor
+	rhoLiouville_dot= zeros(N,1);
+	% 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
+	rhoLiouville_dot(1)=1;
+endfunction
+
+function kappa=susceptibility(wi, rhoLiouville, dipole_elements, E_field)
+% calculate susceptibility for the field at given frequency index
+	Nlevels=( size(dipole_elements.linear)(1) );
+	rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
+	kappa.linear=0;
+	kappa.left=0;
+	kappa.right=0;
+
+	kappa.linear = sum( sum( transpose(dipole_elements.linear) .* rho ) );
+	kappa.left   = sum( sum( transpose(dipole_elements.left)   .* rho ) );
+	kappa.right  = sum( sum( transpose(dipole_elements.right)  .* rho ) );
+
+	kappa.linear /= E_field.linear( wi ); 
+	kappa.left   /= E_field.left( wi );
+	kappa.right  /= E_field.right( wi ); 
+endfunction
+
+function index=freq2index(freq, modulation_freq) 
+% convert modulation freq to its index in the modulation_freq vector
+	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_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;
+	L0m                = atom_properties.L0m              ; 
+	polarizability_m   = atom_properties.polarizability_m ; 
+	dipole_elements    = atom_properties.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, E_field);
+	xi_linear = xi.linear;
+	xi_right  = xi.right;
+	xi_left   = xi.left;
+endfunction		
+
+function [dEdz_coef_linear, dEdz_coef_left, dEdz_coef_right] =dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field)
+% complex absorption coefficient 
+% at given E_field frequency component
+	Nlevels=( size(dipole_elements.linear)(1) );
+	rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
+
+	dEdz_coef_linear = sum( sum( transpose(dipole_elements.linear) .* rho ) );
+	dEdz_coef_left   = sum( sum( transpose(dipole_elements.left)   .* rho ) );
+	dEdz_coef_right  = sum( sum( transpose(dipole_elements.right)  .* rho ) );
+
+endfunction 
+
+function [dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem)
+% complex absorption coefficient for each field
+	global atom_properties;
+	L0m                = atom_properties.L0m              ; 
+	polarizability_m   = atom_properties.polarizability_m ; 
+	dipole_elements    = atom_properties.dipole_elements  ; 
+
+	E_field            = atom_field_problem.E_field          ; 
+	modulation_freq    = atom_field_problem.modulation_freq  ; 
+	Nfreq=length(modulation_freq);
+
+	rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
+	for wi=1:Nfreq
+	 [dEdz_linear(wi), dEdz_left(wi), dEdz_right(wi)] = dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field);
+	endfor
+endfunction		
+
+function relative_transmission=total_relative_transmission(atom_field_problem)
+% summed across all frequencies field absorption coefficient
+	global atom_properties;
+	[dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem);
+	%dEdz_total= dEdz_linear .* conj(dEdz_linear) ...
+									 %+dEdz_left   .* conj(dEdz_left) ...
+									 %+dEdz_right  .* conj(dEdz_right);
+	%total_absorption = sum(dEdz_total);
+
+	%total_absorption = |E|^2 - | E - dEdz | ^2		
+	E_field.linear = atom_field_problem.E_field.linear;
+	E_field.right  = atom_field_problem.E_field.right;
+	E_field.left   = atom_field_problem.E_field.left;
+	
+	modulation_freq    = atom_field_problem.modulation_freq  ; 
+	pos_freq_indexes=(modulation_freq>0);
+	% WARNING INTRODUCED HERE BUT MUST BE DEFINE IN ATOM PROPERTIES FILE
+	%n_atoms - proportional to a number of interacting atoms
+	n_atoms=500;
+	transmited_intensities_vector = ...	
+		 abs(E_field.linear + n_atoms*(1i)*dEdz_linear).^2 ...
+		+abs(E_field.right  + n_atoms*(1i)*dEdz_right).^2 ...
+		+abs(E_field.left   + n_atoms*(1i)*dEdz_left).^2;
+	transmited_intensities_vector = transmited_intensities_vector(pos_freq_indexes);
+	input_intensities_vector =  ...	
+		 abs(E_field.linear).^2   ...
+		+abs(E_field.right).^2    ...
+		+abs(E_field.left).^2;
+	input_intensities_vector = input_intensities_vector(pos_freq_indexes);
+		
+	relative_transmission = sum(transmited_intensities_vector) / sum(input_intensities_vector);
+endfunction
+
+% create full list of atom modulation frequencies from positive light frequency amplitudes
+function [modulation_freq, E_field_amplitudes] = ...
+light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(...
+	positive_light_frequencies, positive_light_field_amplitudes )
+	% we should add 0 frequency as first element of our frequencies list
+	modulation_freq    = cat(2, 0, positive_light_frequencies, -positive_light_frequencies);
+	% negative frequencies have complex conjugated light fields amplitudes
+	E_field_amplitudes.left = cat(2, 0, positive_light_field_amplitudes.left, conj(positive_light_field_amplitudes.left) );
+	E_field_amplitudes.right = cat(2, 0, positive_light_field_amplitudes.right, conj(positive_light_field_amplitudes.right) );
+	E_field_amplitudes.linear = cat(2, 0, positive_light_field_amplitudes.linear, conj(positive_light_field_amplitudes.linear) );
+endfunction	
+
+
+
+% vim: ts=2:sw=2:fdm=indent