From f25ef4ac36b03b1c4e26ef0aa14905734c2db9fe Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Wed, 16 Nov 2011 16:33:49 -0500 Subject: clean up of unused files and link to the master copy --- faraday/basis_transformation.m | 58 +----- faraday/dipole_elementRb87D1line.m | 178 +----------------- faraday/ouput_psr_vs_detuning_combo.m | 36 ---- faraday/output_results.m | 31 ---- faraday/output_xi_results.m | 46 ----- faraday/propagation_problem.m | 79 ++++++++ faraday/run_me.m | 2 +- faraday/useful_constants.m | 7 +- faraday/useful_functions.m | 338 +--------------------------------- 9 files changed, 84 insertions(+), 691 deletions(-) mode change 100644 => 120000 faraday/basis_transformation.m mode change 100644 => 120000 faraday/dipole_elementRb87D1line.m delete mode 100644 faraday/ouput_psr_vs_detuning_combo.m delete mode 100644 faraday/output_results.m delete mode 100644 faraday/output_xi_results.m create mode 100644 faraday/propagation_problem.m mode change 100644 => 120000 faraday/useful_constants.m mode change 100644 => 120000 faraday/useful_functions.m diff --git a/faraday/basis_transformation.m b/faraday/basis_transformation.m deleted file mode 100644 index 8b8eb04..0000000 --- a/faraday/basis_transformation.m +++ /dev/null @@ -1,57 +0,0 @@ -1; - -% matrix of circular to linear transformation -% [x, y, z]' = lin2circ * [r, l, z]' -function transformation_matrix = circ2lin() - transformation_matrix = ... - [ ... - [ 1/sqrt(2), 1/sqrt(2), 0]; ... - [ 1i/sqrt(2),-1i/sqrt(2), 0]; ... - [ 0, 0, 1] ... - ]; -endfunction - - -% matrix of linear to circular transformation -% [r, l, z]' = lin2circ * [x, y, z]' -function transformation_matrix = lin2circ() - transformation_matrix = ... - [ ... - [ 1/sqrt(2), -1i/sqrt(2), 0]; ... - [ 1/sqrt(2), 1i/sqrt(2), 0]; ... - [ 0, 0, 1] ... - ]; -endfunction -% linear basis rotation -% x axis untouched -% z and y rotated by angle theta around 'x' axis -% [x_new, y_new, z_new]' = oldlin2newlin * [x_old, y_old, z_old]' -function oldlin2newlin_m = oldlin2newlin(theta) - oldlin2newlin_m = [ ... - [ 1, 0, 0]; ... - [ 0, cos(theta), -sin(theta)]; ... - [ 0, sin(theta), cos(theta)]... - ]; -endfunction - -% rotate x polarized light by angle phi around -% light propagation axis (Z) -function [E_field_x, E_field_y] = rotXpolarization(phi, E_field_linear) - % important negative frequency behave as they rotate in opposite direction - E_field_x=cos(phi)*E_field_linear; - E_field_y=sin(phi)*E_field_linear; -endfunction - -% transform x,y,z linearly polarized light in the lab/light system coordinate -% to left, right, linear along z atom system of coordinate -% atom magnetic field is along new axis Z wich is at angle theta with respect to -% light propagation direction -function E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq) - coord_transf_m = lin2circ() * oldlin2newlin(theta); - E_field_pos_freq.right = coord_transf_m(1,1)*E_field_lab_pos_freq.x + coord_transf_m(1,2)*E_field_lab_pos_freq.y + coord_transf_m(1,3)*E_field_lab_pos_freq.z; - E_field_pos_freq.left = coord_transf_m(2,1)*E_field_lab_pos_freq.x + coord_transf_m(2,2)*E_field_lab_pos_freq.y + coord_transf_m(2,3)*E_field_lab_pos_freq.z; - E_field_pos_freq.linear = coord_transf_m(3,1)*E_field_lab_pos_freq.x + coord_transf_m(3,2)*E_field_lab_pos_freq.y + coord_transf_m(3,3)*E_field_lab_pos_freq.z; -endfunction - - -% vim: ts=2:sw=2:fdm=indent diff --git a/faraday/basis_transformation.m b/faraday/basis_transformation.m new file mode 120000 index 0000000..28d8beb --- /dev/null +++ b/faraday/basis_transformation.m @@ -0,0 +1 @@ +../basis_transformation.m \ No newline at end of file diff --git a/faraday/dipole_elementRb87D1line.m b/faraday/dipole_elementRb87D1line.m deleted file mode 100644 index 9235274..0000000 --- a/faraday/dipole_elementRb87D1line.m +++ /dev/null @@ -1,177 +0,0 @@ -function d=dipole_elementRb87D1line(Fl,ml,Fu,mu) -% Fl, ml are F and m quantum numbers of lower state -% Fu, mu are F and m quantum numbers of upper state -% F is total momentum and m is projection - d.left = 0; %default return value - d.linear = 0; %default return value - d.right = 0; %default return value - if ( mu==(ml+1) ) - % sigma plus polarization - % ------ Fl=2 -> Fu=2 -------- - if ( (ml==-2) & (Fl==2) & (Fu==2) ) - d.right=sqrt(1/6); - endif - if ( (ml==-1) & (Fl==2) & (Fu==2) ) - d.right=sqrt(1/4); - endif - if ( (ml== 0) & (Fl==2) & (Fu==2) ) - d.right=sqrt(1/4); - endif - if ( (ml== 1) & (Fl==2) & (Fu==2) ) - d.right=sqrt(1/6); - endif - if ( (ml== 2) & (Fl==2) & (Fu==2) ) - d.right=0; - endif - % ------ Fl=2 -> Fu=1 -------- - if ( (ml==-2) & (Fl==2) & (Fu==1) ) - d.right=sqrt(1/2); - endif - if ( (ml==-1) & (Fl==2) & (Fu==1) ) - d.right=sqrt(1/4); - endif - if ( (ml== 0) & (Fl==2) & (Fu==1) ) - d.right=sqrt(1/12); - endif - if ( (ml== 1) & (Fl==2) & (Fu==1) ) - d.right=0; - endif - if ( (ml== 2) & (Fl==2) & (Fu==1) ) - d.right=0; - endif - % ------ Fl=1 -> Fu=2 -------- - if ( (ml==-1) & (Fl==1) & (Fu==2) ) - d.right = -sqrt(1/12); - endif - if ( (ml== 0) & (Fl==1) & (Fu==2) ) - d.right = -sqrt(1/4); - endif - if ( (ml== 1) & (Fl==1) & (Fu==2) ) - d.right = -sqrt(1/2); - endif - % ------ Fl=1 -> Fu=1 -------- - if ( (ml==-1) & (Fl==1) & (Fu==1) ) - d.right = -sqrt(1/12); - endif - if ( (ml== 0) & (Fl==1) & (Fu==1) ) - d.right = -sqrt(1/12); - endif - if ( (ml== 1) & (Fl==1) & (Fu==1) ) - d.right = 0; - endif - endif - if ( mu==(ml+0) ) - % pi polarization - % ------ Fl=2 -> Fu=2 -------- - if ( (ml==-2) & (Fl==2) & (Fu==2) ) - d.linear=-sqrt(1/3); - endif - if ( (ml==-1) & (Fl==2) & (Fu==2) ) - d.linear=-sqrt(1/12); - endif - if ( (ml== 0) & (Fl==2) & (Fu==2) ) - d.linear=0; - endif - if ( (ml== 1) & (Fl==2) & (Fu==2) ) - d.linear=sqrt(1/12); - endif - if ( (ml== 2) & (Fl==2) & (Fu==2) ) - d.linear=sqrt(1/3); - endif - % ------ Fl=2 -> Fu=1 -------- - if ( (ml==-2) & (Fl==2) & (Fu==1) ) - d.linear = 0; - endif - if ( (ml==-1) & (Fl==2) & (Fu==1) ) - d.linear=sqrt(1/4); - endif - if ( (ml== 0) & (Fl==2) & (Fu==1) ) - d.linear=sqrt(1/3); - endif - if ( (ml== 1) & (Fl==2) & (Fu==1) ) - d.linear=sqrt(1/4); - endif - if ( (ml== 2) & (Fl==2) & (Fu==1) ) - d.linear = 0; - endif - % ------ Fl=1 -> Fu=2 -------- - if ( (ml==-1) & (Fl==1) & (Fu==2) ) - d.linear = sqrt(1/4); - endif - if ( (ml== 0) & (Fl==1) & (Fu==2) ) - d.linear = sqrt(1/2); - endif - if ( (ml== 1) & (Fl==1) & (Fu==2) ) - d.linear = sqrt(1/4); - endif - % ------ Fl=1 -> Fu=1 -------- - if ( (ml==-1) & (Fl==1) & (Fu==1) ) - d.linear = sqrt(1/12); - endif - if ( (ml== 0) & (Fl==1) & (Fu==1) ) - d.linear = 0; - endif - if ( (ml== 1) & (Fl==1) & (Fu==1) ) - d.linear = -sqrt(1/12); - endif - endif - if ( mu==(ml-1) ) - % sigma minus polarization - % ------ Fl=2 -> Fu=2 -------- - if ( (ml==-2) & (Fl==2) & (Fu==2) ) - d.left = 0; - endif - if ( (ml==-1) & (Fl==2) & (Fu==2) ) - d.left = -sqrt(1/6); - endif - if ( (ml== 0) & (Fl==2) & (Fu==2) ) - d.left = -sqrt(1/4); - endif - if ( (ml== 1) & (Fl==2) & (Fu==2) ) - d.left = -sqrt(1/4); - endif - if ( (ml== 2) & (Fl==2) & (Fu==2) ) - d.left = -sqrt(1/6); - endif - % ------ Fl=2 -> Fu=1 -------- - if ( (ml==-2) & (Fl==2) & (Fu==1) ) - d.left = 0; - endif - if ( (ml==-1) & (Fl==2) & (Fu==1) ) - d.left = 0; - endif - if ( (ml== 0) & (Fl==2) & (Fu==1) ) - d.left = sqrt(1/12); - endif - if ( (ml== 1) & (Fl==2) & (Fu==1) ) - d.left = sqrt(1/4); - endif - if ( (ml== 2) & (Fl==2) & (Fu==1) ) - d.left = sqrt(1/2); - endif - % ------ Fl=1 -> Fu=2 -------- - if ( (ml==-1) & (Fl==1) & (Fu==2) ) - d.left = -sqrt(1/2); - endif - if ( (ml== 0) & (Fl==1) & (Fu==2) ) - d.left = -sqrt(1/4); - endif - if ( (ml== 1) & (Fl==1) & (Fu==2) ) - d.left = -sqrt(1/12); - endif - % ------ Fl=1 -> Fu=1 -------- - if ( (ml==-1) & (Fl==1) & (Fu==1) ) - d.left = 0; - endif - if ( (ml== 0) & (Fl==1) & (Fu==1) ) - d.left = sqrt(1/12); - endif - if ( (ml== 1) & (Fl==1) & (Fu==1) ) - d.left = sqrt(1/12); - endif - endif -endfunction - - - -% vim: ts=2:sw=2:fdm=indent diff --git a/faraday/dipole_elementRb87D1line.m b/faraday/dipole_elementRb87D1line.m new file mode 120000 index 0000000..6a11745 --- /dev/null +++ b/faraday/dipole_elementRb87D1line.m @@ -0,0 +1 @@ +../dipole_elementRb87D1line.m \ No newline at end of file diff --git a/faraday/ouput_psr_vs_detuning_combo.m b/faraday/ouput_psr_vs_detuning_combo.m deleted file mode 100644 index 154effc..0000000 --- a/faraday/ouput_psr_vs_detuning_combo.m +++ /dev/null @@ -1,36 +0,0 @@ -function ret=ouput_psr_vs_detuning_combo( ... - psr_rad_tnEp_pos_el, psr_rad_tnEp_neg_el, ... - psr_rad_smEp_pos_el, psr_rad_smEp_neg_el, ... - psr_rad_lgEp_pos_el, psr_rad_lgEp_neg_el, ... - psr_rad_grEp_pos_el, psr_rad_grEp_neg_el ... - , detuning_freq, B_field, theta, phi, psi_el ... - , output_dir) - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - figure(6); - tn_cl='blue'; - sm_cl='green'; - lg_cl='magenta'; - gt_cl='red'; - - plot ( ... - detuning_freq, psr_rad_tnEp_pos_el, '-;tn,pos el;', "linewidth",3, 'color',tn_cl - , detuning_freq, psr_rad_tnEp_neg_el, '-;tn,neg el;', "linewidth",1, 'color',tn_cl - , detuning_freq, psr_rad_smEp_pos_el, '-;sm,pos el;', "linewidth",3, 'color',sm_cl - , detuning_freq, psr_rad_smEp_neg_el, '-;sm,neg el;', "linewidth",1, 'color',sm_cl - , detuning_freq, psr_rad_lgEp_pos_el, '-;lg,pos el;', "linewidth",3, 'color',lg_cl - , detuning_freq, psr_rad_lgEp_neg_el, '-;lg,neg el;', "linewidth",1, 'color',lg_cl - , detuning_freq, psr_rad_grEp_pos_el, '-;gr,pos el;', "linewidth",3, 'color',gt_cl - , detuning_freq, psr_rad_grEp_neg_el, '-;gr,neg el;', "linewidth",1, 'color',gt_cl - ); - - str_title=sprintf("PSR. B field=%.5f Gauss, ellipticity=%.2f rad, theta=%.2f, phi=%.2f", B_field, psi_el, theta, phi); - title(str_title); - xlabel('detuning, MHz'); - ylabel('PSR, radians'); - fname=data_file_name(output_dir, 'PSR.','pdf', B_field, theta,phi,psi_el); - print(fname); - - ret=true; -endfunction - diff --git a/faraday/output_results.m b/faraday/output_results.m deleted file mode 100644 index eca6085..0000000 --- a/faraday/output_results.m +++ /dev/null @@ -1,31 +0,0 @@ -1; - - -load '/tmp/relative_transmission_vs_detuning.mat' ; - -figure(1); -hold off; -zoom_factor=1; -plot(detuning_freq, zoom_factor*(relative_transmission_vs_detuning), '-'); -title("relative transmission"); -xlabel("two photon detuning"); - -%load 'xi_vs_detuning.mat' ; - -%figure(1); - %hold off; - %plot(detuning_freq, imag(xi_linear), '-1;linear;'); - %hold on; - %plot(detuning_freq, imag(xi_left), '-2;left;'); - %plot(detuning_freq, imag(xi_right), '-3;right;'); - %title("probe absorption"); - %hold off; -%figure(2); - %hold off; - %plot(detuning_freq, real(xi_linear), '-1;linear;'); - %hold on; - %plot(detuning_freq, real(xi_left), '-2;left;'); - %plot(detuning_freq, real(xi_right), '-3;right;'); - %title("probe dispersion"); - %hold off; - diff --git a/faraday/output_xi_results.m b/faraday/output_xi_results.m deleted file mode 100644 index 7dc3855..0000000 --- a/faraday/output_xi_results.m +++ /dev/null @@ -1,46 +0,0 @@ -1; - - -load '/tmp/xi_vs_detuning.mat' ; - -figure(1); -hold off; -plot(detuning_freq, real(xi_left-xi_right), '-'); -title("differential real xi"); -xlabel("two photon detuning"); - -figure(2); -hold off; -plot(detuning_freq, imag(xi_left-xi_right), '-'); -title("differential imag xi"); -xlabel("two photon detuning"); - -figure(3); -hold off; -plot(detuning_freq, imag(xi_left), '-', detuning_freq, imag(xi_right), '-'); -title("imag xi"); -xlabel("two photon detuning"); - -figure(4); -hold off; -plot(detuning_freq, real(xi_left), '-', detuning_freq, real(xi_right), '-'); -title("real xi"); -xlabel("two photon detuning"); - -%figure(1); - %hold off; - %plot(detuning_freq, imag(xi_linear), '-1;linear;'); - %hold on; - %plot(detuning_freq, imag(xi_left), '-2;left;'); - %plot(detuning_freq, imag(xi_right), '-3;right;'); - %title("probe absorption"); - %hold off; -%figure(2); - %hold off; - %plot(detuning_freq, real(xi_linear), '-1;linear;'); - %hold on; - %plot(detuning_freq, real(xi_left), '-2;left;'); - %plot(detuning_freq, real(xi_right), '-3;right;'); - %title("probe dispersion"); - %hold off; - diff --git a/faraday/propagation_problem.m b/faraday/propagation_problem.m new file mode 100644 index 0000000..74b4471 --- /dev/null +++ b/faraday/propagation_problem.m @@ -0,0 +1,79 @@ +function [xi_linear, xi_left, xi_right]=propagation_problem(detuning_freq, Ep, psi_el, B_field, theta, phi) +% calculates transmission if light polarizations vs B field in the cell +% for given laser probe and B fields array +% Probe field defined by field strength (Ep) and ellipticity angle (pse_el) +% Magnetic field defined by magnitude (B_field) and angles theta and phi. +% +% Note: it is expensive to recalculate atom property for each given B_field strength +% so run as many calculation for constant magnetic field as possible + +t0 = clock (); % we will use this latter to calculate elapsed time + + +% load useful functions; +useful_functions; + +% some physical constants +useful_constants; + +basis_transformation; % load subroutines + + +% load atom energy levels and decay description +rb87_D1_line; + +% calculate E_field independent properties of the atom +% to be used as sub matrix templates for Liouville operator matrix +[L0m, polarizability_m]=L0_and_polarization_submatrices( ... +Nlevels, ... +H0, g_decay, g_dephasing, dipole_elements ... +); + +global atom_properties; +atom_properties.L0m=L0m; +atom_properties.polarizability_m=polarizability_m; +atom_properties.dipole_elements=dipole_elements; + +%light_positive_freq = [wp]; +E_field_drive = [0 ]; +E_field_probe = [Ep ]; +E_field_zero = [0 ]; +E_field_lab_pos_freq.linear = E_field_zero + (1.00000+0.00000i)*E_field_probe + (1.00000+0.00000i)*E_field_drive; + +% we define light as linearly polarized +% where phi is angle between light polarization and axis x +% only sign of modulation frequency is important now +% we define actual frequency later on +[E_field_lab_pos_freq.x, E_field_lab_pos_freq.y] = rotXpolarization(phi, E_field_lab_pos_freq.linear); +% we add required ellipticity +E_field_lab_pos_freq.x*=exp(I*psi_el); +E_field_lab_pos_freq.y*=exp(-I*psi_el); +E_field_lab_pos_freq.z=E_field_zero; + +E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq); + + +wp0=w_pf1-w_sf2; %Fg=2 -> Fe=1 +wp=wp0+detuning_freq; +light_positive_freq=[ wp]; +% we calculate dc and negative frequencies as well as amplitudes +[modulation_freq, E_field] = ... + light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(... + light_positive_freq, E_field_pos_freq); +freq_index=freq2index(wp,modulation_freq); + +atom_field_problem.E_field = E_field; +atom_field_problem.modulation_freq = modulation_freq; +atom_field_problem.freq_index = freq_index; + +problems_cell_array=atom_field_problem; + +[xi_linear, xi_left, xi_right]=susceptibility_steady_state_at_freq( problems_cell_array); + + +elapsed_time = etime (clock (), t0) +return + +endfunction + +% vim: ts=2:sw=2:fdm=indent diff --git a/faraday/run_me.m b/faraday/run_me.m index 8a29973..a8984e1 100644 --- a/faraday/run_me.m +++ b/faraday/run_me.m @@ -11,7 +11,7 @@ B_fields=linspace(-zeeman_splitting/gmg, zeeman_splitting/gmg, Nsteps); % phi is angle between linear polarization and axis x %phi=pi/4; -phi=0; +phi=pi/2; % theta is angle between lab z axis (light propagation direction) and magnetic field axis (z') theta=0; % psi_el is the ellipticity parameter (phase difference between left and right polarization) diff --git a/faraday/useful_constants.m b/faraday/useful_constants.m deleted file mode 100644 index 82d1d1a..0000000 --- a/faraday/useful_constants.m +++ /dev/null @@ -1,6 +0,0 @@ -1; - -% some physical constants -hbar=1; -im_one=0+1i; -el_charge=1; diff --git a/faraday/useful_constants.m b/faraday/useful_constants.m new file mode 120000 index 0000000..fc5b92d --- /dev/null +++ b/faraday/useful_constants.m @@ -0,0 +1 @@ +../useful_constants.m \ No newline at end of file diff --git a/faraday/useful_functions.m b/faraday/useful_functions.m deleted file mode 100644 index eb3e880..0000000 --- a/faraday/useful_functions.m +++ /dev/null @@ -1,337 +0,0 @@ -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 diff --git a/faraday/useful_functions.m b/faraday/useful_functions.m new file mode 120000 index 0000000..a2d3237 --- /dev/null +++ b/faraday/useful_functions.m @@ -0,0 +1 @@ +../useful_functions.m \ No newline at end of file -- cgit v1.2.3