diff options
| -rw-r--r-- | liouville.m | 14 | ||||
| -rw-r--r-- | useful_functions.m | 42 |
2 files changed, 30 insertions, 26 deletions
diff --git a/liouville.m b/liouville.m index 4c78557..4ce714c 100644 --- a/liouville.m +++ b/liouville.m @@ -1,4 +1,5 @@ 1; +clear all; t0 = clock (); % we will use this latter to calculate elapsed time @@ -30,14 +31,15 @@ kappa_p =zeros(1,N_detun_steps+1); kappa_m =zeros(1,N_detun_steps+1); detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps; +% now we create Liouville indexes list +[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq); +rhoLiouville=zeros(N,1); % calculate E_field independent properties of athe atom % to be used as sub matrix templates for Liouville operator matrix [L0m, polarizability_m]=L0_and_polarization_submatrices( ... - Nlevels*Nlevels, ... - H0, g_decay, g_dephasing, dipole_elements, ... - E_field, ... - modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ... + Nlevels, ... + H0, g_decay, g_dephasing, dipole_elements ... ); for detuning_p_cntr=1:N_detun_steps+1; wp0=w12; @@ -77,8 +79,8 @@ for detuning_p_cntr=1:N_detun_steps+1; %rho_d=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq); % drive frequency %rho_m=rhoOfFreq(rhoLiouville, 4, Nlevels, Nfreq); % opposite sideband frequency - kappa_p(detuning_p_cntr)=sucseptibility(2, rhoLiouville, dipole_elements, Nlevels, Nfreq); - %kappa_m(detuning_p_cntr)=sucseptibility(4, rhoLiouville, dipole_elements, Nlevels, Nfreq); + kappa_p(detuning_p_cntr)=susceptibility(2, 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); diff --git a/useful_functions.m b/useful_functions.m index 606fe49..c0fdd7e 100644 --- a/useful_functions.m +++ b/useful_functions.m @@ -1,15 +1,15 @@ 1; -% calculate total decay for particular level taking in account all branches function ret=decay_total(g_decay,i) +% calculate total decay for particular level taking in account all branches ret=0; for k=1:size(g_decay)(1) ret=ret+g_decay(i,k); endfor endfunction -% kroneker delta symbol function ret=kron_delta(i,j) +% kroneker delta symbol if ((i==j)) ret=1; else @@ -17,9 +17,9 @@ function ret=kron_delta(i,j) endif endfunction +function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq) % this function create from Liouville density vector % the density matrix with given modulation frequency -function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq) rho=zeros(Nlevels); for r=1:Nlevels for c=1:Nlevels @@ -28,10 +28,10 @@ function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq) endfor endfunction -% we unwrap density matrix to Liouville density vector and assign all possible +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) -function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq) N=Nfreq*Nlevels*Nlevels; rhoLiouville_w=zeros(N,1); rhoLiouville_r=zeros(N,1); @@ -51,21 +51,24 @@ function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matr 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 -function [L0m, polarizability_m]=L0_and_polarization_submatrices( ... - rho_size, ... - H0, g_decay, g_dephasing, dipole_elements, ... - E_field, ... - modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ... - ) %------------------------- 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); + % note that L0 and decay parts depend only on combination of indexes % jk,mn but repeats itself for every frequency - L0m=zeros(rho_size); % (NxN)x(NxN) matrix + 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 @@ -95,14 +98,14 @@ function [L0m, polarizability_m]=L0_and_polarization_submatrices( ... L0m=-im_one/hbar*L0m - decay_part_m; endfunction -% Liouville operator matrix construction -% based on recipe from Eugeniy Mikhailov thesis 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 @@ -161,15 +164,15 @@ function L=Liouville_operator_matrix( ... 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 -function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ... - N, L, ... - modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c) for i=1:N w2i=rhoLiouville_w(i); m=rhoLiouville_r(i); @@ -188,9 +191,8 @@ function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ... rhoLiouville_dot(1)=1; endfunction -% calculate sucseptibility for the field at given frequency index -function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq) - +function kappa=susceptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq) +% calculate susceptibility for the field at given frequency index rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq); kappa=0; for i=1:Nlevels |