diff options
| -rw-r--r-- | liouville.m | 79 | ||||
| -rw-r--r-- | useful_functions.m | 19 |
2 files changed, 50 insertions, 48 deletions
diff --git a/liouville.m b/liouville.m index 204c1a3..4c78557 100644 --- a/liouville.m +++ b/liouville.m @@ -29,53 +29,60 @@ detuning_freq=zeros(1,N_detun_steps+1); 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; -for detuning_p_cntr=1:N_detun_steps+1; - -wp0=w12; -detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1); -wp=wp0+detuning_p; -wm=wd-(wp-wd); -%modulation_freq=[0, wp, wd, wm, -wp, -wd, -wm, wp-wd, wd-wp]; -%E_field =[0, Ep, Ed, Em, Epc, Edc, Emc, 0, 0 ]; -modulation_freq=[0, wp, wd, -wp, -wd, wp-wd, wd-wp]; -E_field =[0, Ep, Ed, Epc, Edc, 0, 0 ]; -Nfreq=length(modulation_freq); -% now we create Liouville indexes list -[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq); -rhoLiouville=zeros(N,1); -% Liouville operator matrix construction -L=Liouville_operator_matrix( - N, - H0, g_decay, g_dephasing, dipole_elements, - E_field, - modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c +% 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 ... ); +for detuning_p_cntr=1:N_detun_steps+1; + wp0=w12; + detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1); + wp=wp0+detuning_p; + wm=wd-(wp-wd); + %modulation_freq=[0, wp, wd, wm, -wp, -wd, -wm, wp-wd, wd-wp]; + %E_field =[0, Ep, Ed, Em, Epc, Edc, Emc, 0, 0 ]; + modulation_freq=[0, wp, wd, -wp, -wd, wp-wd, wd-wp]; + E_field =[0, Ep, Ed, Epc, Edc, 0, 0 ]; + 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]=constran_rho_and_match_L( - N, L, - 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; + %solving for density matrix vector + rhoLiouville=L\rhoLiouville_dot; -%rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq); % 0 frequency, -%rho_p=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq); % probe frequency -%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); -detuning_freq(detuning_p_cntr)=detuning_p; + %rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq); % 0 frequency, + %rho_p=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq); % probe frequency + %rho_d=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq); % drive frequency + %rho_m=rhoOfFreq(rhoLiouville, 4, Nlevels, Nfreq); % opposite sideband frequency -%kappa_p_re=real(kappa_p); -%kappa_p_im=imag(kappa_p); + kappa_p(detuning_p_cntr)=sucseptibility(2, rhoLiouville, dipole_elements, Nlevels, Nfreq); + %kappa_m(detuning_p_cntr)=sucseptibility(4, rhoLiouville, dipole_elements, Nlevels, Nfreq); + detuning_freq(detuning_p_cntr)=detuning_p; + %kappa_p_re=real(kappa_p); + %kappa_p_im=imag(kappa_p); endfor figure(1); plot(detuning_freq, real(kappa_p)); title("probe dispersion"); figure(2); plot(detuning_freq, imag(kappa_p)); title("probe absorption"); @@ -83,3 +90,5 @@ figure(2); plot(detuning_freq, imag(kappa_p)); title("probe absorption"); %figure(4); plot(detuning_freq, imag(kappa_m)); title("off resonant absorption"); elapsed_time = etime (clock (), t0) + +% vim: ts=2:sw=2:fdm=indent diff --git a/useful_functions.m b/useful_functions.m index 56989df..a1135f7 100644 --- a/useful_functions.m +++ b/useful_functions.m @@ -55,7 +55,7 @@ endfunction % 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( ... +function [L0m, polarizability_m]=L0_and_polarization_submatrices( ... rho_size, ... H0, g_decay, g_dephasing, dipole_elements, ... E_field, ... @@ -69,7 +69,7 @@ function [L0m, polarization_m]=L0_and_polarization_submatrices( ... 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 + polarizability_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 @@ -85,7 +85,7 @@ function [L0m, polarization_m]=L0_and_polarization_submatrices( ... 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) ); + polarizability_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; @@ -95,7 +95,7 @@ endfunction % based on recipe from Eugeniy Mikhailov thesis function L=Liouville_operator_matrix( N, - H0, g_decay, g_dephasing, dipole_elements, + L0m, polarizability_m, E_field, modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ) @@ -112,13 +112,6 @@ function L=Liouville_operator_matrix( % 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 @@ -150,7 +143,7 @@ function L=Liouville_operator_matrix( % 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); + -im_one/hbar*polarizability_m* E_field(w3i); endif % diagonal elements are self modulated % due to rotating wave approximation @@ -170,7 +163,7 @@ endfunction % 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]=constran_rho_and_match_L( +function [rhoLiouville_dot, L]=constrain_rho_and_match_L( N, L, modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c) for i=1:N |