diff options
Diffstat (limited to 'psr/useful_functions.m')
| -rw-r--r-- | psr/useful_functions.m | 337 |
1 files changed, 337 insertions, 0 deletions
diff --git a/psr/useful_functions.m b/psr/useful_functions.m new file mode 100644 index 0000000..eb3e880 --- /dev/null +++ b/psr/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 |