summaryrefslogtreecommitdiff
path: root/faraday
diff options
context:
space:
mode:
Diffstat (limited to 'faraday')
l---------[-rw-r--r--]faraday/basis_transformation.m58
l---------[-rw-r--r--]faraday/dipole_elementRb87D1line.m178
-rw-r--r--faraday/ouput_psr_vs_detuning_combo.m36
-rw-r--r--faraday/output_results.m31
-rw-r--r--faraday/output_xi_results.m46
-rw-r--r--faraday/propagation_problem.m79
-rw-r--r--faraday/run_me.m2
l---------[-rw-r--r--]faraday/useful_constants.m7
l---------[-rw-r--r--]faraday/useful_functions.m338
9 files changed, 84 insertions, 691 deletions
diff --git a/faraday/basis_transformation.m b/faraday/basis_transformation.m
index 8b8eb04..28d8beb 100644..120000
--- a/faraday/basis_transformation.m
+++ b/faraday/basis_transformation.m
@@ -1,57 +1 @@
-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
+../basis_transformation.m \ No newline at end of file
diff --git a/faraday/dipole_elementRb87D1line.m b/faraday/dipole_elementRb87D1line.m
index 9235274..6a11745 100644..120000
--- a/faraday/dipole_elementRb87D1line.m
+++ b/faraday/dipole_elementRb87D1line.m
@@ -1,177 +1 @@
-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
+../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
index 82d1d1a..fc5b92d 100644..120000
--- a/faraday/useful_constants.m
+++ b/faraday/useful_constants.m
@@ -1,6 +1 @@
-1;
-
-% some physical constants
-hbar=1;
-im_one=0+1i;
-el_charge=1;
+../useful_constants.m \ No newline at end of file
diff --git a/faraday/useful_functions.m b/faraday/useful_functions.m
index eb3e880..a2d3237 100644..120000
--- a/faraday/useful_functions.m
+++ b/faraday/useful_functions.m
@@ -1,337 +1 @@
-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
+../useful_functions.m \ No newline at end of file
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 22:38:21 +0000