diff options
| -rw-r--r-- | useful_functions.m | 118 |
1 files changed, 74 insertions, 44 deletions
diff --git a/useful_functions.m b/useful_functions.m index 38841a4..56989df 100644 --- a/useful_functions.m +++ b/useful_functions.m @@ -51,7 +51,48 @@ function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matr endfunction +% 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, polarization_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; + % 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 + 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 + polarization_m=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(k,n)-H0(n,k)*kron_delta(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(j,m)*kron_delta(k,n) \ + - kron_delta(m,n)*kron_delta(j,k)*g_decay(m,j) ; + polarization_m(p,s)= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) ); + endfor + endfor + 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, H0, g_decay, g_dephasing, dipole_elements, @@ -64,71 +105,58 @@ function L=Liouville_operator_matrix( Nfreq=length(modulation_freq); % Lets be supper smart and speed up L matrix construction - % it has a lot of voids. + % 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; + + [L0m, polarization_m]=L0_and_polarization_submatrices( ... + rho_size, ... + H0, g_decay, g_dephasing, dipole_elements, ... + E_field, ... + modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ... + ); + + + % Liouville matrix operator has Nlevels*Nlevels blocks + % which governed by the same modulation frequency for p_freq_cntr=1:Nfreq for s_freq_cntr=1:Nfreq p0=1+(p_freq_cntr-1)*rho_size; s0=1+(s_freq_cntr-1)*rho_size; % we guaranteed to know frequency of final and initial rhoLiouville - w1i=rhoLiouville_w(p0); - w2i=rhoLiouville_w(s0); + w1i=rhoLiouville_w(p0); % final w_jk=modulation_freq(w1i); + 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 - decay_part=0; - Lt=0; for w3i=1:Nfreq % at most we should have only one matching frequency w_iner=modulation_freq(w3i); if ((w_iner == w_l)) % yey, requested modulation frequency exist - % lets do L submatrix filling - for p_cntr=0:(rho_size-1) - p=p0+p_cntr; - for s_cntr=0:(rho_size-1) - s=s0+s_cntr; - - decay_part=0; - Lt=0; - - j=rhoLiouville_r(p); - k=rhoLiouville_c(p); - m=rhoLiouville_r(s); - n=rhoLiouville_c(s); - - %such frequency exist in the list of modulation frequencies - if ((w_iner == 0)) - % calculate unperturbed part (Hamiltonian without EM field) - L0=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m); - decay_part=\ - ( decay_total(g_decay,k)/2 + decay_total(g_decay,j)/2 + g_dephasing(j,k) )* kron_delta(j,m)*kron_delta(k,n) \ - - kron_delta(m,n)*kron_delta(j,k)*g_decay(m,j) ; - Lt=L0; - else - % calculate perturbed part (Hamiltonian with EM field) - % in othe word interactive part of Hamiltonian - Li= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) )*E_field(w3i); - Lt=Li; - endif - %Lt=-im_one/hbar*Lt*kron_delta(w_jk-w_iner,w_mn); % above if should be done only if kron_delta is not zero - % no need for above kron_delta since the same conditon checked in the outer if statement - Lt=-im_one/hbar*Lt - decay_part; - % diagonal elements are self modulated - if ((p == s)) - Lt+=-im_one*w_jk; - endif - Lps=Lt; - L(p,s)=Lps; - endfor - endfor + % lets do L sub matrix filling + + %such frequency exist in the list of modulation frequencies + if ((w_iner == 0)) + % calculate unperturbed part (Hamiltonian without EM field) + L(p0:p0+rho_size-1,s0:s0+rho_size-1)=L0m; + else + % 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*polarization_m* E_field(w3i); + endif + % diagonal elements are self modulated + % due to rotating wave approximation + if ((p0 == s0)) + L(p0:p0+rho_size-1,s0:s0+rho_size-1)+= -im_one*w_jk*eye(rho_size); + endif endif endfor endfor @@ -174,3 +202,5 @@ function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq) endfor endfor endfunction + +% vim: ts=2:sw=2:fdm=indent |