diff options
| -rw-r--r-- | useful_functions.m | 87 |
1 files changed, 55 insertions, 32 deletions
diff --git a/useful_functions.m b/useful_functions.m index 1765dd9..38841a4 100644 --- a/useful_functions.m +++ b/useful_functions.m @@ -59,19 +59,24 @@ function L=Liouville_operator_matrix( modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ) %------------------------- - L=zeros(N); % NxN matrix useful_constants; + L=zeros(N); % NxN matrix Nfreq=length(modulation_freq); - for p=1:N - for s=1:N - j=rhoLiouville_r(p); - k=rhoLiouville_c(p); - m=rhoLiouville_r(s); - n=rhoLiouville_c(s); + % Lets be supper smart and speed up L matrix construction + % 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; + 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(p); - w2i=rhoLiouville_w(s); + w1i=rhoLiouville_w(p0); + w2i=rhoLiouville_w(s0); w_jk=modulation_freq(w1i); w_mn=modulation_freq(w2i); % thus we know L matrix element frequency which we need to match @@ -81,33 +86,51 @@ function L=Liouville_operator_matrix( decay_part=0; Lt=0; for w3i=1:Nfreq + % at most we should have only one matching frequency w_iner=modulation_freq(w3i); - decay_part=0; if ((w_iner == w_l)) - %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; + % 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 endif endfor - if ((p == s)) - Lt+=-im_one*w_jk; - endif - Lps=Lt; - L(p,s)=Lps; endfor endfor endfunction |