summaryrefslogtreecommitdiff
path: root/faraday
diff options
context:
space:
mode:
Diffstat (limited to 'faraday')
-rw-r--r--faraday/basis_transformation.m57
-rw-r--r--faraday/data_file_name.m5
-rw-r--r--faraday/dipole_elementRb87D1line.m177
-rw-r--r--faraday/drawing.gpbin0 -> 132234 bytes
-rw-r--r--faraday/field_description.m14
-rw-r--r--faraday/make_representative_psr_vs_detuning_for_given_B_and_psi_el.m35
-rw-r--r--faraday/ouput_psr_vs_detuning_combo.m36
-rw-r--r--faraday/output_psr_results_vs_detuning.m59
-rw-r--r--faraday/output_psr_results_vs_power.m70
-rw-r--r--faraday/output_results.m31
-rw-r--r--faraday/output_xi_results.m46
-rw-r--r--faraday/psr_vs_detuning.m161
-rw-r--r--faraday/psr_vs_detuning.ps1711
-rw-r--r--faraday/psr_vs_detuning_combo.m149
-rw-r--r--faraday/psr_vs_power.m147
-rw-r--r--faraday/rb87_D1_line.m148
-rw-r--r--faraday/useful_constants.m6
-rw-r--r--faraday/useful_functions.m337
18 files changed, 3189 insertions, 0 deletions
diff --git a/faraday/basis_transformation.m b/faraday/basis_transformation.m
new file mode 100644
index 0000000..8b8eb04
--- /dev/null
+++ b/faraday/basis_transformation.m
@@ -0,0 +1,57 @@
+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/data_file_name.m b/faraday/data_file_name.m
new file mode 100644
index 0000000..5d6e106
--- /dev/null
+++ b/faraday/data_file_name.m
@@ -0,0 +1,5 @@
+function fname=data_file_name(prefix, title, suffix, B_field, theta,phi,psi_el)
+ str_title=sprintf("%s B field=%.5f Gauss, ellipticity=%.2f rad, theta=%.2f, phi=%.2f", title, B_field, psi_el, theta, phi);
+ fname=strcat(prefix, str_title, '.', suffix);
+endfunction
+
diff --git a/faraday/dipole_elementRb87D1line.m b/faraday/dipole_elementRb87D1line.m
new file mode 100644
index 0000000..9235274
--- /dev/null
+++ b/faraday/dipole_elementRb87D1line.m
@@ -0,0 +1,177 @@
+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/drawing.gp b/faraday/drawing.gp
new file mode 100644
index 0000000..3bcddd3
--- /dev/null
+++ b/faraday/drawing.gp
Binary files differ
diff --git a/faraday/field_description.m b/faraday/field_description.m
new file mode 100644
index 0000000..71a1cdb
--- /dev/null
+++ b/faraday/field_description.m
@@ -0,0 +1,14 @@
+1;
+
+%EM field definition
+%Ed=0.2; %drive
+%Edc=conj(Ed);
+Ep=0.2; %probe
+Epc=conj(Ep);
+Em=-Ep; % opposite sideband (resulting from EOM modulation of drive)
+Emc=conj(Em);
+%wd=w13;
+%wp=w12;
+%wm=wd-(wp-wd);
+%light_positive_freq = [wp, wd, wp-wd];
+%E_field = [Ep, Ed, 0 ];
diff --git a/faraday/make_representative_psr_vs_detuning_for_given_B_and_psi_el.m b/faraday/make_representative_psr_vs_detuning_for_given_B_and_psi_el.m
new file mode 100644
index 0000000..46d8262
--- /dev/null
+++ b/faraday/make_representative_psr_vs_detuning_for_given_B_and_psi_el.m
@@ -0,0 +1,35 @@
+function ...
+ [ ...
+ 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 ...
+ ] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el)
+ %%%%%%%%%%%%%%%%%%%%%
+ %printf("%f\n", B_field);
+ Ep=sqrt(0.01);
+ [psr_rad_tnEp_pos_el]=psr_vs_detuning(detuning_freq, Ep, psi_el, B_field, theta, phi) ;
+ [psr_rad_tnEp_neg_el]=psr_vs_detuning(detuning_freq, Ep,-psi_el, B_field, theta, phi) ;
+
+ Ep=sqrt(0.1);
+ [psr_rad_smEp_pos_el]=psr_vs_detuning(detuning_freq, Ep, psi_el, B_field, theta, phi) ;
+ [psr_rad_smEp_neg_el]=psr_vs_detuning(detuning_freq, Ep,-psi_el, B_field, theta, phi) ;
+
+ Ep=sqrt(1.0);
+ [psr_rad_lgEp_pos_el]=psr_vs_detuning(detuning_freq, Ep, psi_el, B_field, theta, phi) ;
+ [psr_rad_lgEp_neg_el]=psr_vs_detuning(detuning_freq, Ep,-psi_el, B_field, theta, phi) ;
+
+ Ep=sqrt(10.0);
+ [psr_rad_grEp_pos_el]=psr_vs_detuning(detuning_freq, Ep, psi_el, B_field, theta, phi) ;
+ [psr_rad_grEp_neg_el]=psr_vs_detuning(detuning_freq, Ep,-psi_el, B_field, theta, phi) ;
+
+
+ fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el)
+ save(fname, ...
+ 'detuning_freq', 'B_field', 'theta', 'phi', 'psi_el', ...
+ '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' );
+
+endfunction
diff --git a/faraday/ouput_psr_vs_detuning_combo.m b/faraday/ouput_psr_vs_detuning_combo.m
new file mode 100644
index 0000000..154effc
--- /dev/null
+++ b/faraday/ouput_psr_vs_detuning_combo.m
@@ -0,0 +1,36 @@
+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_psr_results_vs_detuning.m b/faraday/output_psr_results_vs_detuning.m
new file mode 100644
index 0000000..93912ae
--- /dev/null
+++ b/faraday/output_psr_results_vs_detuning.m
@@ -0,0 +1,59 @@
+function psr_rad=output_psr_results_vs_detuning()
+load '/tmp/xi_vs_detuning.mat' ;
+
+Er=(1+I*xi_right)*E_field_pos_freq.right;
+El=(1+I*xi_left) *E_field_pos_freq.left;
+
+Ex=(Er+El)/sqrt(2);
+Ey=I*(Er-El)/sqrt(2);
+
+%extra rotation to compensate rotation due to ellipticity
+% actually no need for it since x-polarization shifts by positive phase
+% and y-pol by negative phase
+%el_rot=0*psi_el;
+%Ex=cos(el_rot)*Ex-sin(el_rot)*Ey;
+%Ey=sin(el_rot)*Ex+cos(el_rot)*Ey;
+
+Ipos=(abs(Ey).^2)/2;
+Ineg=(abs(Ex).^2)/2;
+
+figure(1);
+hold off;
+plot(detuning_freq, real(xi_left-xi_right), '-');
+title("differential real xi");
+xlabel("two photon detuning (MHz)");
+
+figure(2);
+hold off;
+plot(detuning_freq, imag(xi_left-xi_right), '-');
+title("differential imag xi");
+xlabel("two photon detuning (MHz)");
+
+figure(3);
+hold off;
+plot(detuning_freq, real(xi_left), '-', detuning_freq, real(xi_right), '-');
+title("real xi");
+xlabel("two photon detuning (MHz)");
+
+figure(4);
+hold off;
+plot(detuning_freq, imag(xi_left), '-', detuning_freq, imag(xi_right), '-');
+title("imag xi");
+xlabel("two photon detuning (MHz)");
+
+figure(5);
+hold off;
+I_probe=E_field_probe^2;
+psr_rad=(Ipos-Ineg)/(2*I_probe);
+plot(detuning_freq, psr_rad, '-');
+subt_str=sprintf("Laser Rabi freq normalized to upper state decay %.3f, ellipticity %.1f degree, \n B field ground level splitting %.3f Gauss", I_probe, psi_el*180/pi, B_field);
+title(cstrcat("BPD normilized PSR signal at F_g=2 to F_e=1,2.\n ",subt_str) );
+xlabel("two photon detuning (MHz)");
+ylabel("PSR (radians)");
+
+print("psr_vs_detuning.ps");
+
+fname= sprintf("psr_vs_detuning_Fg=2toFe=1,2_Ip=%.3f_el_%.1f_B=%.3fG.mat", I_probe, psi_el*180/pi,B_field);
+save(fname,'detuning_freq', 'psr_rad');
+
+return;
diff --git a/faraday/output_psr_results_vs_power.m b/faraday/output_psr_results_vs_power.m
new file mode 100644
index 0000000..2d28565
--- /dev/null
+++ b/faraday/output_psr_results_vs_power.m
@@ -0,0 +1,70 @@
+1;
+
+
+load '/tmp/xi_vs_power.mat' ;
+
+Er=(1+I*xi_right)*E_field_pos_freq.right;
+El=(1+I*xi_left) *E_field_pos_freq.left;
+
+Ex=(Er+El)/sqrt(2);
+Ey=I*(Er-El)/sqrt(2);
+
+%extra rotation to compensate rotation due to ellipticity
+% actually no need for it since x-polarization shifts by positive phase
+% and y-pol by negative phase
+%el_rot=0*psi_el;
+%Ex=cos(el_rot)*Ex-sin(el_rot)*Ey;
+%Ey=sin(el_rot)*Ex+cos(el_rot)*Ey;
+
+Ipos=(abs(Ey).^2)/2;
+Ineg=(abs(Ex).^2)/2;
+
+figure(1);
+hold off;
+plot(Ep.^2, real(xi_left-xi_right), '-');
+title("differential real xi");
+xlabel("two photon detuning");
+
+figure(2);
+hold off;
+plot(Ep.^2, imag(xi_left-xi_right), '-');
+title("differential imag xi");
+xlabel("two photon detuning");
+
+figure(3);
+hold off;
+plot(Ep.^2, imag(xi_left), '-', Ep, imag(xi_right), '-');
+title("imag xi");
+xlabel("two photon detuning");
+
+figure(4);
+hold off;
+plot(Ep.^2, real(xi_left), '-', Ep.^2, real(xi_right), '-');
+title("real xi");
+xlabel("two photon detuning");
+
+figure(5);
+hold off;
+%plot(Ep.^2, (Ipos-Ineg), '-');
+semilogx(Ep.^2, (Ipos-Ineg), '-');
+%semilogx(Ep.^2, (Ipos-Ineg)./(Ep.^2), '-');
+title("BPD signal 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/output_results.m b/faraday/output_results.m
new file mode 100644
index 0000000..eca6085
--- /dev/null
+++ b/faraday/output_results.m
@@ -0,0 +1,31 @@
+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
new file mode 100644
index 0000000..7dc3855
--- /dev/null
+++ b/faraday/output_xi_results.m
@@ -0,0 +1,46 @@
+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/psr_vs_detuning.m b/faraday/psr_vs_detuning.m
new file mode 100644
index 0000000..acf7d43
--- /dev/null
+++ b/faraday/psr_vs_detuning.m
@@ -0,0 +1,161 @@
+function [psr_rad]=psr_vs_detuning(detuning_freq, Ep, psi_el, B_field, theta, phi)
+% calculates psr in rad vs detunings of the probe field
+% for given laser probe and B field.
+% 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;
+%four_levels_with_polarization;
+%four_levels;
+%three_levels;
+%two_levels;
+
+% load EM field description
+%field_description;
+
+%Nfreq=length(modulation_freq);
+
+
+
+%tune probe frequency
+detuning_p=0;
+%detuning_p_min=-B_field*gmg*4; % span +/-4 Zeeman splitting
+%detuning_freq=zeros(1,N_detun_steps+1);
+N_detun_steps=length(detuning_freq);
+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;
+
+fprintf (stderr, "calculating atom properties\n");
+fflush (stderr);
+pfile='atomic_B_field.mat'; % the parent file where B_field is stored. This is the parameter for calculated L0_and_polarization_submatrices
+cfile='L0m_and_polarizability_calculated.mat'; % the child file to which calculated matrices are written
+need_update=false;
+[s, err_p, msg] = stat (pfile);
+if(err_p)
+ %file does not exist
+ need_update=true;
+ else
+ B_field_cur=B_field;
+ load (pfile); % loading old B_field value
+ if (B_field ~= B_field_cur)
+ % old and current B field are different
+ B_field=B_field_cur;
+ need_update=true;
+ else
+ need_update=false;
+ endif
+endif
+
+[s, err, msg] = stat (cfile);
+if(err)
+ %file does not exist
+ need_update=true;
+endif;
+if ( !need_update)
+ % matrices already calculated and up to date, all we need to load them
+ load(cfile);
+ else
+ % 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 ...
+ );
+ save(pfile, 'B_field');
+ save(cfile, 'L0m', 'polarizability_m');
+ endif
+elapsed_time = etime (clock (), t0);
+fprintf (stderr, "elapsed time so far is %.3f sec\n",elapsed_time);
+fflush (stderr);
+
+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;
+%E_field_lab_pos_freq.right = E_field_zero + (0.00000+0.00000i)*E_field_probe + (0.00000+0.00000i)*E_field_drive;
+%E_field_lab_pos_freq.left = E_field_zero + (0.00000+0.00000i)*E_field_probe + (0.00000+0.00000i)*E_field_drive;
+
+% phi is angle between linear polarization and axis x
+%phi=pi*2/8;
+% 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)
+%psi_el=-30/180*pi;
+
+% 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);
+
+
+fprintf (stderr, "tuning laser in forloop to set conditions vs detuning\n");
+fflush (stderr);
+for detuning_p_cntr=1:length(detuning_freq);
+ wp0=w_pf1-w_sf2; %Fg=2 -> Fe=1
+ %wd=w_pf1-w_hpf_ground;
+ %detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1);
+ detuning_p=detuning_freq(detuning_p_cntr);
+ wp=wp0+detuning_p;
+ light_positive_freq=[ wp];
+ % we calculate dc and negative frequiencies 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{detuning_p_cntr}=atom_field_problem;
+
+
+ %kappa_p(detuning_p_cntr)=susceptibility_steady_state_at_freq( atom_field_problem);
+ %detuning_freq(detuning_p_cntr)=detuning_p;
+endfor
+
+save '/tmp/problem_definition.mat' problems_cell_array atom_properties detuning_freq ;
+fprintf (stderr, "now really hard calculations begin\n");
+fflush (stderr);
+% once we define all problems the main job is done here
+[xi_linear, xi_left, xi_right]=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
+%[xi_linear, xi_left, xi_right]=cellfun( @susceptibility_steady_state_at_freq, problems_cell_array);
+
+%save '/tmp/relative_transmission_vs_detuning.mat' detuning_freq relative_transmission_vs_detuning;
+save '/tmp/xi_vs_detuning.mat' detuning_freq xi_linear xi_left xi_right E_field_pos_freq E_field_probe B_field psi_el;
+
+psr_rad=output_psr_results_vs_detuning;
+
+elapsed_time = etime (clock (), t0)
+return
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/faraday/psr_vs_detuning.ps b/faraday/psr_vs_detuning.ps
new file mode 100644
index 0000000..a5695fa
--- /dev/null
+++ b/faraday/psr_vs_detuning.ps
@@ -0,0 +1,1711 @@
+%!PS-Adobe-2.0
+%%Title: psr_vs_detuning.ps
+%%Creator: gnuplot 4.2 patchlevel 2
+%%CreationDate: Sat Jul 10 11:37:59 2010
+%%DocumentFonts: (atend)
+%%BoundingBox: 50 50 612 792
+%%Orientation: Portrait
+%%Pages: (atend)
+%%EndComments
+%%BeginProlog
+/gnudict 256 dict def
+gnudict begin
+%
+% The following 6 true/false flags may be edited by hand if required
+% The unit line width may also be changed
+%
+/Color false def
+/Blacktext false def
+/Solid false def
+/Dashlength 1 def
+/Landscape false def
+/Level1 false def
+/Rounded false def
+/TransparentPatterns false def
+/gnulinewidth 5.000 def
+/userlinewidth gnulinewidth def
+%
+/vshift -46 def
+/dl1 {
+ 10.0 Dashlength mul mul
+ Rounded { currentlinewidth 0.75 mul sub dup 0 le { pop 0.01 } if } if
+} def
+/dl2 {
+ 10.0 Dashlength mul mul
+ Rounded { currentlinewidth 0.75 mul add } if
+} def
+/hpt_ 31.5 def
+/vpt_ 31.5 def
+/hpt hpt_ def
+/vpt vpt_ def
+Level1 {} {
+/SDict 10 dict def
+systemdict /pdfmark known not {
+ userdict /pdfmark systemdict /cleartomark get put
+} if
+SDict begin [
+ /Title (psr_vs_detuning.ps)
+ /Subject (gnuplot plot)
+ /Creator (gnuplot 4.2 patchlevel 2 )
+ /Author (Eugeniy Mikhailov)
+% /Producer (gnuplot)
+% /Keywords ()
+ /CreationDate (Sat Jul 10 11:37:59 2010)
+ /DOCINFO pdfmark
+end
+} ifelse
+%
+% Gnuplot Prolog Version 4.2 (August 2006)
+%
+/M {moveto} bind def
+/L {lineto} bind def
+/R {rmoveto} bind def
+/V {rlineto} bind def
+/N {newpath moveto} bind def
+/Z {closepath} bind def
+/C {setrgbcolor} bind def
+/f {rlineto fill} bind def
+/vpt2 vpt 2 mul def
+/hpt2 hpt 2 mul def
+/Lshow {currentpoint stroke M 0 vshift R
+ Blacktext {gsave 0 setgray show grestore} {show} ifelse} def
+/Rshow {currentpoint stroke M dup stringwidth pop neg vshift R
+ Blacktext {gsave 0 setgray show grestore} {show} ifelse} def
+/Cshow {currentpoint stroke M dup stringwidth pop -2 div vshift R
+ Blacktext {gsave 0 setgray show grestore} {show} ifelse} def
+/UP {dup vpt_ mul /vpt exch def hpt_ mul /hpt exch def
+ /hpt2 hpt 2 mul def /vpt2 vpt 2 mul def} def
+/DL {Color {setrgbcolor Solid {pop []} if 0 setdash}
+ {pop pop pop 0 setgray Solid {pop []} if 0 setdash} ifelse} def
+/BL {stroke userlinewidth 2 mul setlinewidth
+ Rounded {1 setlinejoin 1 setlinecap} if} def
+/AL {stroke userlinewidth 2 div setlinewidth
+ Rounded {1 setlinejoin 1 setlinecap} if} def
+/UL {dup gnulinewidth mul /userlinewidth exch def
+ dup 1 lt {pop 1} if 10 mul /udl exch def} def
+/PL {stroke userlinewidth setlinewidth
+ Rounded {1 setlinejoin 1 setlinecap} if} def
+% Default Line colors
+/LCw {1 1 1} def
+/LCb {0 0 0} def
+/LCa {0 0 0} def
+/LC0 {1 0 0} def
+/LC1 {0 1 0} def
+/LC2 {0 0 1} def
+/LC3 {1 0 1} def
+/LC4 {0 1 1} def
+/LC5 {1 1 0} def
+/LC6 {0 0 0} def
+/LC7 {1 0.3 0} def
+/LC8 {0.5 0.5 0.5} def
+% Default Line Types
+/LTw {PL [] 1 setgray} def
+/LTb {BL [] LCb DL} def
+/LTa {AL [1 udl mul 2 udl mul] 0 setdash LCa setrgbcolor} def
+/LT0 {PL [] LC0 DL} def
+/LT1 {PL [4 dl1 2 dl2] LC1 DL} def
+/LT2 {PL [2 dl1 3 dl2] LC2 DL} def
+/LT3 {PL [1 dl1 1.5 dl2] LC3 DL} def
+/LT4 {PL [6 dl1 2 dl2 1 dl1 2 dl2] LC4 DL} def
+/LT5 {PL [3 dl1 3 dl2 1 dl1 3 dl2] LC5 DL} def
+/LT6 {PL [2 dl1 2 dl2 2 dl1 6 dl2] LC6 DL} def
+/LT7 {PL [1 dl1 2 dl2 6 dl1 2 dl2 1 dl1 2 dl2] LC7 DL} def
+/LT8 {PL [2 dl1 2 dl2 2 dl1 2 dl2 2 dl1 2 dl2 2 dl1 4 dl2] LC8 DL} def
+/Pnt {stroke [] 0 setdash gsave 1 setlinecap M 0 0 V stroke grestore} def
+/Dia {stroke [] 0 setdash 2 copy vpt add M
+ hpt neg vpt neg V hpt vpt neg V
+ hpt vpt V hpt neg vpt V closepath stroke
+ Pnt} def
+/Pls {stroke [] 0 setdash vpt sub M 0 vpt2 V
+ currentpoint stroke M
+ hpt neg vpt neg R hpt2 0 V stroke
+ } def
+/Box {stroke [] 0 setdash 2 copy exch hpt sub exch vpt add M
+ 0 vpt2 neg V hpt2 0 V 0 vpt2 V
+ hpt2 neg 0 V closepath stroke
+ Pnt} def
+/Crs {stroke [] 0 setdash exch hpt sub exch vpt add M
+ hpt2 vpt2 neg V currentpoint stroke M
+ hpt2 neg 0 R hpt2 vpt2 V stroke} def
+/TriU {stroke [] 0 setdash 2 copy vpt 1.12 mul add M
+ hpt neg vpt -1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt 1.62 mul V closepath stroke
+ Pnt} def
+/Star {2 copy Pls Crs} def
+/BoxF {stroke [] 0 setdash exch hpt sub exch vpt add M
+ 0 vpt2 neg V hpt2 0 V 0 vpt2 V
+ hpt2 neg 0 V closepath fill} def
+/TriUF {stroke [] 0 setdash vpt 1.12 mul add M
+ hpt neg vpt -1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt 1.62 mul V closepath fill} def
+/TriD {stroke [] 0 setdash 2 copy vpt 1.12 mul sub M
+ hpt neg vpt 1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt -1.62 mul V closepath stroke
+ Pnt} def
+/TriDF {stroke [] 0 setdash vpt 1.12 mul sub M
+ hpt neg vpt 1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt -1.62 mul V closepath fill} def
+/DiaF {stroke [] 0 setdash vpt add M
+ hpt neg vpt neg V hpt vpt neg V
+ hpt vpt V hpt neg vpt V closepath fill} def
+/Pent {stroke [] 0 setdash 2 copy gsave
+ translate 0 hpt M 4 {72 rotate 0 hpt L} repeat
+ closepath stroke grestore Pnt} def
+/PentF {stroke [] 0 setdash gsave
+ translate 0 hpt M 4 {72 rotate 0 hpt L} repeat
+ closepath fill grestore} def
+/Circle {stroke [] 0 setdash 2 copy
+ hpt 0 360 arc stroke Pnt} def
+/CircleF {stroke [] 0 setdash hpt 0 360 arc fill} def
+/C0 {BL [] 0 setdash 2 copy moveto vpt 90 450 arc} bind def
+/C1 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 0 90 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C2 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 90 180 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C3 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 0 180 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C4 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 180 270 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C5 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 0 90 arc
+ 2 copy moveto
+ 2 copy vpt 180 270 arc closepath fill
+ vpt 0 360 arc} bind def
+/C6 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 90 270 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C7 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 0 270 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C8 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 270 360 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C9 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 270 450 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C10 {BL [] 0 setdash 2 copy 2 copy moveto vpt 270 360 arc closepath fill
+ 2 copy moveto
+ 2 copy vpt 90 180 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C11 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 0 180 arc closepath fill
+ 2 copy moveto
+ 2 copy vpt 270 360 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C12 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 180 360 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C13 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 0 90 arc closepath fill
+ 2 copy moveto
+ 2 copy vpt 180 360 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/C14 {BL [] 0 setdash 2 copy moveto
+ 2 copy vpt 90 360 arc closepath fill
+ vpt 0 360 arc} bind def
+/C15 {BL [] 0 setdash 2 copy vpt 0 360 arc closepath fill
+ vpt 0 360 arc closepath} bind def
+/Rec {newpath 4 2 roll moveto 1 index 0 rlineto 0 exch rlineto
+ neg 0 rlineto closepath} bind def
+/Square {dup Rec} bind def
+/Bsquare {vpt sub exch vpt sub exch vpt2 Square} bind def
+/S0 {BL [] 0 setdash 2 copy moveto 0 vpt rlineto BL Bsquare} bind def
+/S1 {BL [] 0 setdash 2 copy vpt Square fill Bsquare} bind def
+/S2 {BL [] 0 setdash 2 copy exch vpt sub exch vpt Square fill Bsquare} bind def
+/S3 {BL [] 0 setdash 2 copy exch vpt sub exch vpt2 vpt Rec fill Bsquare} bind def
+/S4 {BL [] 0 setdash 2 copy exch vpt sub exch vpt sub vpt Square fill Bsquare} bind def
+/S5 {BL [] 0 setdash 2 copy 2 copy vpt Square fill
+ exch vpt sub exch vpt sub vpt Square fill Bsquare} bind def
+/S6 {BL [] 0 setdash 2 copy exch vpt sub exch vpt sub vpt vpt2 Rec fill Bsquare} bind def
+/S7 {BL [] 0 setdash 2 copy exch vpt sub exch vpt sub vpt vpt2 Rec fill
+ 2 copy vpt Square fill Bsquare} bind def
+/S8 {BL [] 0 setdash 2 copy vpt sub vpt Square fill Bsquare} bind def
+/S9 {BL [] 0 setdash 2 copy vpt sub vpt vpt2 Rec fill Bsquare} bind def
+/S10 {BL [] 0 setdash 2 copy vpt sub vpt Square fill 2 copy exch vpt sub exch vpt Square fill
+ Bsquare} bind def
+/S11 {BL [] 0 setdash 2 copy vpt sub vpt Square fill 2 copy exch vpt sub exch vpt2 vpt Rec fill
+ Bsquare} bind def
+/S12 {BL [] 0 setdash 2 copy exch vpt sub exch vpt sub vpt2 vpt Rec fill Bsquare} bind def
+/S13 {BL [] 0 setdash 2 copy exch vpt sub exch vpt sub vpt2 vpt Rec fill
+ 2 copy vpt Square fill Bsquare} bind def
+/S14 {BL [] 0 setdash 2 copy exch vpt sub exch vpt sub vpt2 vpt Rec fill
+ 2 copy exch vpt sub exch vpt Square fill Bsquare} bind def
+/S15 {BL [] 0 setdash 2 copy Bsquare fill Bsquare} bind def
+/D0 {gsave translate 45 rotate 0 0 S0 stroke grestore} bind def
+/D1 {gsave translate 45 rotate 0 0 S1 stroke grestore} bind def
+/D2 {gsave translate 45 rotate 0 0 S2 stroke grestore} bind def
+/D3 {gsave translate 45 rotate 0 0 S3 stroke grestore} bind def
+/D4 {gsave translate 45 rotate 0 0 S4 stroke grestore} bind def
+/D5 {gsave translate 45 rotate 0 0 S5 stroke grestore} bind def
+/D6 {gsave translate 45 rotate 0 0 S6 stroke grestore} bind def
+/D7 {gsave translate 45 rotate 0 0 S7 stroke grestore} bind def
+/D8 {gsave translate 45 rotate 0 0 S8 stroke grestore} bind def
+/D9 {gsave translate 45 rotate 0 0 S9 stroke grestore} bind def
+/D10 {gsave translate 45 rotate 0 0 S10 stroke grestore} bind def
+/D11 {gsave translate 45 rotate 0 0 S11 stroke grestore} bind def
+/D12 {gsave translate 45 rotate 0 0 S12 stroke grestore} bind def
+/D13 {gsave translate 45 rotate 0 0 S13 stroke grestore} bind def
+/D14 {gsave translate 45 rotate 0 0 S14 stroke grestore} bind def
+/D15 {gsave translate 45 rotate 0 0 S15 stroke grestore} bind def
+/DiaE {stroke [] 0 setdash vpt add M
+ hpt neg vpt neg V hpt vpt neg V
+ hpt vpt V hpt neg vpt V closepath stroke} def
+/BoxE {stroke [] 0 setdash exch hpt sub exch vpt add M
+ 0 vpt2 neg V hpt2 0 V 0 vpt2 V
+ hpt2 neg 0 V closepath stroke} def
+/TriUE {stroke [] 0 setdash vpt 1.12 mul add M
+ hpt neg vpt -1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt 1.62 mul V closepath stroke} def
+/TriDE {stroke [] 0 setdash vpt 1.12 mul sub M
+ hpt neg vpt 1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt -1.62 mul V closepath stroke} def
+/PentE {stroke [] 0 setdash gsave
+ translate 0 hpt M 4 {72 rotate 0 hpt L} repeat
+ closepath stroke grestore} def
+/CircE {stroke [] 0 setdash
+ hpt 0 360 arc stroke} def
+/Opaque {gsave closepath 1 setgray fill grestore 0 setgray closepath} def
+/DiaW {stroke [] 0 setdash vpt add M
+ hpt neg vpt neg V hpt vpt neg V
+ hpt vpt V hpt neg vpt V Opaque stroke} def
+/BoxW {stroke [] 0 setdash exch hpt sub exch vpt add M
+ 0 vpt2 neg V hpt2 0 V 0 vpt2 V
+ hpt2 neg 0 V Opaque stroke} def
+/TriUW {stroke [] 0 setdash vpt 1.12 mul add M
+ hpt neg vpt -1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt 1.62 mul V Opaque stroke} def
+/TriDW {stroke [] 0 setdash vpt 1.12 mul sub M
+ hpt neg vpt 1.62 mul V
+ hpt 2 mul 0 V
+ hpt neg vpt -1.62 mul V Opaque stroke} def
+/PentW {stroke [] 0 setdash gsave
+ translate 0 hpt M 4 {72 rotate 0 hpt L} repeat
+ Opaque stroke grestore} def
+/CircW {stroke [] 0 setdash
+ hpt 0 360 arc Opaque stroke} def
+/BoxFill {gsave Rec 1 setgray fill grestore} def
+/Density {
+ /Fillden exch def
+ currentrgbcolor
+ /ColB exch def /ColG exch def /ColR exch def
+ /ColR ColR Fillden mul Fillden sub 1 add def
+ /ColG ColG Fillden mul Fillden sub 1 add def
+ /ColB ColB Fillden mul Fillden sub 1 add def
+ ColR ColG ColB setrgbcolor} def
+/BoxColFill {gsave Rec PolyFill} def
+/PolyFill {gsave Density fill grestore grestore} def
+/h {rlineto rlineto rlineto gsave fill grestore} bind def
+%
+% PostScript Level 1 Pattern Fill routine for rectangles
+% Usage: x y w h s a XX PatternFill
+% x,y = lower left corner of box to be filled
+% w,h = width and height of box
+% a = angle in degrees between lines and x-axis
+% XX = 0/1 for no/yes cross-hatch
+%
+/PatternFill {gsave /PFa [ 9 2 roll ] def
+ PFa 0 get PFa 2 get 2 div add PFa 1 get PFa 3 get 2 div add translate
+ PFa 2 get -2 div PFa 3 get -2 div PFa 2 get PFa 3 get Rec
+ gsave 1 setgray fill grestore clip
+ currentlinewidth 0.5 mul setlinewidth
+ /PFs PFa 2 get dup mul PFa 3 get dup mul add sqrt def
+ 0 0 M PFa 5 get rotate PFs -2 div dup translate
+ 0 1 PFs PFa 4 get div 1 add floor cvi
+ {PFa 4 get mul 0 M 0 PFs V} for
+ 0 PFa 6 get ne {
+ 0 1 PFs PFa 4 get div 1 add floor cvi
+ {PFa 4 get mul 0 2 1 roll M PFs 0 V} for
+ } if
+ stroke grestore} def
+%
+/languagelevel where
+ {pop languagelevel} {1} ifelse
+ 2 lt
+ {/InterpretLevel1 true def}
+ {/InterpretLevel1 Level1 def}
+ ifelse
+%
+% PostScript level 2 pattern fill definitions
+%
+/Level2PatternFill {
+/Tile8x8 {/PaintType 2 /PatternType 1 /TilingType 1 /BBox [0 0 8 8] /XStep 8 /YStep 8}
+ bind def
+/KeepColor {currentrgbcolor [/Pattern /DeviceRGB] setcolorspace} bind def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop 0 0 M 8 8 L 0 8 M 8 0 L stroke}
+>> matrix makepattern
+/Pat1 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop 0 0 M 8 8 L 0 8 M 8 0 L stroke
+ 0 4 M 4 8 L 8 4 L 4 0 L 0 4 L stroke}
+>> matrix makepattern
+/Pat2 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop 0 0 M 0 8 L
+ 8 8 L 8 0 L 0 0 L fill}
+>> matrix makepattern
+/Pat3 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop -4 8 M 8 -4 L
+ 0 12 M 12 0 L stroke}
+>> matrix makepattern
+/Pat4 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop -4 0 M 8 12 L
+ 0 -4 M 12 8 L stroke}
+>> matrix makepattern
+/Pat5 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop -2 8 M 4 -4 L
+ 0 12 M 8 -4 L 4 12 M 10 0 L stroke}
+>> matrix makepattern
+/Pat6 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop -2 0 M 4 12 L
+ 0 -4 M 8 12 L 4 -4 M 10 8 L stroke}
+>> matrix makepattern
+/Pat7 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop 8 -2 M -4 4 L
+ 12 0 M -4 8 L 12 4 M 0 10 L stroke}
+>> matrix makepattern
+/Pat8 exch def
+<< Tile8x8
+ /PaintProc {0.5 setlinewidth pop 0 -2 M 12 4 L
+ -4 0 M 12 8 L -4 4 M 8 10 L stroke}
+>> matrix makepattern
+/Pat9 exch def
+/Pattern1 {PatternBgnd KeepColor Pat1 setpattern} bind def
+/Pattern2 {PatternBgnd KeepColor Pat2 setpattern} bind def
+/Pattern3 {PatternBgnd KeepColor Pat3 setpattern} bind def
+/Pattern4 {PatternBgnd KeepColor Landscape {Pat5} {Pat4} ifelse setpattern} bind def
+/Pattern5 {PatternBgnd KeepColor Landscape {Pat4} {Pat5} ifelse setpattern} bind def
+/Pattern6 {PatternBgnd KeepColor Landscape {Pat9} {Pat6} ifelse setpattern} bind def
+/Pattern7 {PatternBgnd KeepColor Landscape {Pat8} {Pat7} ifelse setpattern} bind def
+} def
+%
+%
+%End of PostScript Level 2 code
+%
+/PatternBgnd {
+ TransparentPatterns {} {gsave 1 setgray fill grestore} ifelse
+} def
+%
+% Substitute for Level 2 pattern fill codes with
+% grayscale if Level 2 support is not selected.
+%
+/Level1PatternFill {
+/Pattern1 {0.250 Density} bind def
+/Pattern2 {0.500 Density} bind def
+/Pattern3 {0.750 Density} bind def
+/Pattern4 {0.125 Density} bind def
+/Pattern5 {0.375 Density} bind def
+/Pattern6 {0.625 Density} bind def
+/Pattern7 {0.875 Density} bind def
+} def
+%
+% Now test for support of Level 2 code
+%
+Level1 {Level1PatternFill} {Level2PatternFill} ifelse
+%
+/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont
+dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall
+currentdict end definefont pop
+/MFshow {
+ { dup 5 get 3 ge
+ { 5 get 3 eq {gsave} {grestore} ifelse }
+ {dup dup 0 get findfont exch 1 get scalefont setfont
+ [ currentpoint ] exch dup 2 get 0 exch R dup 5 get 2 ne {dup dup 6
+ get exch 4 get {show} {stringwidth pop 0 R} ifelse }if dup 5 get 0 eq
+ {dup 3 get {2 get neg 0 exch R pop} {pop aload pop M} ifelse} {dup 5
+ get 1 eq {dup 2 get exch dup 3 get exch 6 get stringwidth pop -2 div
+ dup 0 R} {dup 6 get stringwidth pop -2 div 0 R 6 get
+ show 2 index {aload pop M neg 3 -1 roll neg R pop pop} {pop pop pop
+ pop aload pop M} ifelse }ifelse }ifelse }
+ ifelse }
+ forall} bind def
+/MFwidth {0 exch { dup 5 get 3 ge { 5 get 3 eq { 0 } { pop } ifelse }
+ {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont
+ 6 get stringwidth pop add} {pop} ifelse} ifelse} forall} bind def
+/MLshow { currentpoint stroke M
+ 0 exch R
+ Blacktext {gsave 0 setgray MFshow grestore} {MFshow} ifelse } bind def
+/MRshow { currentpoint stroke M
+ exch dup MFwidth neg 3 -1 roll R
+ Blacktext {gsave 0 setgray MFshow grestore} {MFshow} ifelse } bind def
+/MCshow { currentpoint stroke M
+ exch dup MFwidth -2 div 3 -1 roll R
+ Blacktext {gsave 0 setgray MFshow grestore} {MFshow} ifelse } bind def
+/XYsave { [( ) 1 2 true false 3 ()] } bind def
+/XYrestore { [( ) 1 2 true false 4 ()] } bind def
+end
+%%EndProlog
+%%Page: 1 1
+gnudict begin
+gsave
+50 50 translate
+0.100 0.100 scale
+0 setgray
+newpath
+(Helvetica) findfont 140 scalefont setfont
+gsave % colour palette begin
+/maxcolors 64 def
+/HSV2RGB { exch dup 0.0 eq {pop exch pop dup dup} % achromatic gray
+ { /HSVs exch def /HSVv exch def 6.0 mul dup floor dup 3 1 roll sub
+ /HSVf exch def /HSVi exch cvi def /HSVp HSVv 1.0 HSVs sub mul def
+ /HSVq HSVv 1.0 HSVs HSVf mul sub mul def
+ /HSVt HSVv 1.0 HSVs 1.0 HSVf sub mul sub mul def
+ /HSVi HSVi 6 mod def 0 HSVi eq {HSVv HSVt HSVp}
+ {1 HSVi eq {HSVq HSVv HSVp}{2 HSVi eq {HSVp HSVv HSVt}
+ {3 HSVi eq {HSVp HSVq HSVv}{4 HSVi eq {HSVt HSVp HSVv}
+ {HSVv HSVp HSVq} ifelse} ifelse} ifelse} ifelse} ifelse
+ } ifelse} def
+/Constrain {
+ dup 0 lt {0 exch pop}{dup 1 gt {1 exch pop} if} ifelse} def
+/YIQ2RGB {
+ 3 copy -1.702 mul exch -1.105 mul add add Constrain 4 1 roll
+ 3 copy -0.647 mul exch -0.272 mul add add Constrain 5 1 roll
+ 0.621 mul exch -0.956 mul add add Constrain 3 1 roll } def
+/CMY2RGB { 1 exch sub exch 1 exch sub 3 2 roll 1 exch sub 3 1 roll exch } def
+/XYZ2RGB { 3 copy -0.9017 mul exch -0.1187 mul add exch 0.0585 mul exch add
+ Constrain 4 1 roll 3 copy -0.0279 mul exch 1.999 mul add exch
+ -0.9844 mul add Constrain 5 1 roll -0.2891 mul exch -0.5338 mul add
+ exch 1.91 mul exch add Constrain 3 1 roll} def
+/SelectSpace {ColorSpace (HSV) eq {HSV2RGB}{ColorSpace (XYZ) eq {
+ XYZ2RGB}{ColorSpace (CMY) eq {CMY2RGB}{ColorSpace (YIQ) eq {YIQ2RGB}
+ if} ifelse} ifelse} ifelse} def
+/InterpolatedColor true def
+/grayindex {/gidx 0 def
+ {GrayA gidx get grayv ge {exit} if /gidx gidx 1 add def} loop} def
+/dgdx {grayv GrayA gidx get sub GrayA gidx 1 sub get
+ GrayA gidx get sub div} def
+/redvalue {RedA gidx get RedA gidx 1 sub get
+ RedA gidx get sub dgdxval mul add} def
+/greenvalue {GreenA gidx get GreenA gidx 1 sub get
+ GreenA gidx get sub dgdxval mul add} def
+/bluevalue {BlueA gidx get BlueA gidx 1 sub get
+ BlueA gidx get sub dgdxval mul add} def
+/interpolate {
+ grayindex grayv GrayA gidx get sub abs 1e-5 le
+ {RedA gidx get GreenA gidx get BlueA gidx get}
+ {/dgdxval dgdx def redvalue greenvalue bluevalue} ifelse} def
+/GrayA [0 .0159 .0317 .0476 .0635 .0794 .0952 .1111 .127 .1429 .1587 .1746
+ .1905 .2063 .2222 .2381 .254 .2698 .2857 .3016 .3175 .3333 .3492 .3651
+ .381 .3968 .4127 .4286 .4444 .4603 .4762 .4921 .5079 .5238 .5397 .5556
+ .5714 .5873 .6032 .619 .6349 .6508 .6667 .6825 .6984 .7143 .7302 .746
+ .7619 .7778 .7937 .8095 .8254 .8413 .8571 .873 .8889 .9048 .9206 .9365
+ .9524 .9683 .9841 1 ] def
+/RedA [.057 .0642 .0715 .0787 .086 .0932 .1004 .1077 .1187 .1559 .1932
+ .2305 .2677 .305 .3423 .3795 .4168 .4541 .4914 .5286 .5659 .6032 .6404
+ .6777 .7054 .7172 .7289 .7406 .7524 .7641 .7759 .7876 .7994 .8111 .8229
+ .8346 .8464 .8581 .8698 .8816 .8627 .8254 .7882 .7509 .7136 .6764 .6391
+ .6018 .5645 .5273 .49 .4527 .4155 .3782 .3409 .3037 .2824 .2634 .2444
+ .2254 .2065 .1875 .1685 .1495 ] def
+/GreenA [.057 .0642 .0715 .0787 .086 .0932 .1004 .1077 .1187 .1559 .1932
+ .2305 .2677 .305 .3423 .3795 .4168 .4541 .4914 .5286 .5659 .6032 .6404
+ .6777 .7054 .7172 .7289 .7406 .7524 .7641 .7759 .7876 .7994 .8111 .8229
+ .8346 .8464 .8581 .8698 .8816 .8627 .8254 .7882 .7509 .7136 .6764 .6391
+ .6018 .5645 .5273 .49 .4527 .4155 .3782 .3409 .3037 .2824 .2634 .2444
+ .2254 .2065 .1875 .1685 .1495 ] def
+/BlueA [.057 .0642 .0715 .0787 .086 .0932 .1004 .1077 .1187 .1559 .1932
+ .2305 .2677 .305 .3423 .3795 .4168 .4541 .4914 .5286 .5659 .6032 .6404
+ .6777 .7054 .7172 .7289 .7406 .7524 .7641 .7759 .7876 .7994 .8111 .8229
+ .8346 .8464 .8581 .8698 .8816 .8627 .8254 .7882 .7509 .7136 .6764 .6391
+ .6018 .5645 .5273 .49 .4527 .4155 .3782 .3409 .3037 .2824 .2634 .2444
+ .2254 .2065 .1875 .1685 .1495 ] def
+/pm3dround {maxcolors 0 gt {dup 1 ge
+ {pop 1} {maxcolors mul floor maxcolors 1 sub div} ifelse} if} def
+/pm3dGamma 1.0 1.5 div def
+/ColorSpace (RGB) def
+Color true and { % COLOUR vs. GRAY map
+ InterpolatedColor { %% Interpolation vs. RGB-Formula
+ /g {stroke pm3dround /grayv exch def interpolate
+ SelectSpace setrgbcolor} bind def
+ }{
+ /g {stroke pm3dround dup cF7 Constrain exch dup cF5 Constrain exch cF15 Constrain
+ SelectSpace setrgbcolor} bind def
+ } ifelse
+}{
+ /g {stroke pm3dround pm3dGamma exp setgray} bind def
+} ifelse
+0.500 UL
+LTb
+394 2016 M
+88 0 V
+4241 0 R
+-88 0 V
+stroke
+310 2016 M
+[ [(Helvetica) 120.0 0.0 true true 0 (-0.0004)]
+] -40.0 MRshow
+0.500 UL
+LTb
+394 2723 M
+88 0 V
+4241 0 R
+-88 0 V
+stroke
+310 2723 M
+[ [(Helvetica) 120.0 0.0 true true 0 (-0.0002)]
+] -40.0 MRshow
+0.500 UL
+LTb
+394 3429 M
+88 0 V
+4241 0 R
+-88 0 V
+stroke
+310 3429 M
+[ [(Helvetica) 120.0 0.0 true true 0 (0)]
+] -40.0 MRshow
+0.500 UL
+LTb
+394 4136 M
+88 0 V
+4241 0 R
+-88 0 V
+stroke
+310 4136 M
+[ [(Helvetica) 120.0 0.0 true true 0 (0.0002)]
+] -40.0 MRshow
+0.500 UL
+LTb
+394 4843 M
+88 0 V
+4241 0 R
+-88 0 V
+stroke
+310 4843 M
+[ [(Helvetica) 120.0 0.0 true true 0 (0.0004)]
+] -40.0 MRshow
+0.500 UL
+LTb
+394 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+394 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (-200)]
+] -40.0 MCshow
+0.500 UL
+LTb
+1115 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+1115 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (0)]
+] -40.0 MCshow
+0.500 UL
+LTb
+1837 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+1837 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (200)]
+] -40.0 MCshow
+0.500 UL
+LTb
+2558 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+2558 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (400)]
+] -40.0 MCshow
+0.500 UL
+LTb
+3280 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+3280 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (600)]
+] -40.0 MCshow
+0.500 UL
+LTb
+4001 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+4001 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (800)]
+] -40.0 MCshow
+0.500 UL
+LTb
+4723 1663 M
+0 88 V
+0 3445 R
+0 -88 V
+stroke
+4723 1523 M
+[ [(Helvetica) 120.0 0.0 true true 0 (1000)]
+] -40.0 MCshow
+0.500 UL
+LTb
+0.500 UL
+LTb
+394 5196 N
+0 -3533 V
+4329 0 V
+0 3533 V
+-4329 0 V
+Z stroke
+LCb setrgbcolor
+-432 3429 M
+currentpoint gsave translate 90 rotate 0 0 moveto
+[ [(Helvetica) 120.0 0.0 true true 0 (PSR \(radians\))]
+] -40.0 MCshow
+grestore
+LTb
+LCb setrgbcolor
+2558 1313 M
+[ [(Helvetica) 120.0 0.0 true true 0 (two photon detuning \(MHz\))]
+] -40.0 MCshow
+LTb
+2558 5686 M
+[ [(Helvetica) 120.0 0.0 true true 0 (BPD normilized PSR signal at F)]
+[(Helvetica) 96.0 -36.0 true true 0 (g)]
+[(Helvetica) 120.0 0.0 true true 0 (=2 to F)]
+[(Helvetica) 96.0 -36.0 true true 0 (e)]
+[(Helvetica) 120.0 0.0 true true 0 (=1,2.)]
+] -28.0 MCshow
+2558 5566 M
+[ [(Helvetica) 120.0 0.0 true true 0 ( Laser Rabi freq normalized to upper state decay 10.000, ellipticity -30.0 degree, )]
+] -40.0 MCshow
+2558 5446 M
+[ [(Helvetica) 120.0 0.0 true true 0 ( B field ground level splitting -1.000 Gauss)]
+] -40.0 MCshow
+1.000 UP
+0.500 UL
+LTb
+0.500 UL
+LT0
+394 3229 M
+4 -1 V
+5 -1 V
+4 -1 V
+4 -2 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -2 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -2 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -2 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 1 V
+4 0 V
+5 0 V
+4 1 V
+4 1 V
+5 0 V
+4 1 V
+4 1 V
+5 1 V
+849 3152 L
+4 1 V
+5 1 V
+4 1 V
+4 2 V
+5 1 V
+4 2 V
+4 1 V
+5 2 V
+4 2 V
+4 2 V
+5 2 V
+4 2 V
+4 2 V
+5 3 V
+4 2 V
+4 3 V
+5 3 V
+4 3 V
+4 3 V
+5 3 V
+4 3 V
+4 4 V
+5 3 V
+4 4 V
+4 4 V
+5 4 V
+4 4 V
+4 4 V
+5 5 V
+4 4 V
+4 5 V
+5 5 V
+4 5 V
+4 5 V
+5 5 V
+4 6 V
+4 5 V
+5 6 V
+4 6 V
+4 6 V
+5 6 V
+4 6 V
+4 6 V
+5 7 V
+4 6 V
+4 7 V
+5 6 V
+4 7 V
+4 7 V
+5 7 V
+4 7 V
+4 7 V
+5 7 V
+4 7 V
+4 8 V
+5 7 V
+4 7 V
+4 8 V
+5 7 V
+4 7 V
+4 7 V
+5 8 V
+4 7 V
+4 7 V
+5 7 V
+4 8 V
+4 7 V
+5 7 V
+4 7 V
+4 7 V
+5 7 V
+4 7 V
+4 7 V
+5 6 V
+4 7 V
+4 6 V
+5 7 V
+4 6 V
+4 6 V
+5 6 V
+4 6 V
+4 6 V
+5 6 V
+4 5 V
+4 6 V
+5 5 V
+4 5 V
+4 5 V
+5 5 V
+4 4 V
+4 5 V
+5 4 V
+4 4 V
+4 4 V
+5 4 V
+4 4 V
+4 4 V
+5 3 V
+4 4 V
+4 3 V
+5 3 V
+4 3 V
+4 3 V
+5 2 V
+1304 3666 L
+4 2 V
+5 3 V
+4 2 V
+4 2 V
+5 2 V
+4 2 V
+4 1 V
+5 2 V
+4 1 V
+4 2 V
+5 1 V
+4 1 V
+4 1 V
+5 1 V
+4 1 V
+4 1 V
+5 1 V
+4 0 V
+4 1 V
+5 0 V
+4 1 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -2 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -2 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -2 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -2 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -2 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -2 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -2 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+1759 3615 L
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -2 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -2 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+2214 3523 L
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+2669 3479 L
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+3124 3448 L
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+3579 3401 L
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 1 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 1 V
+4 0 V
+4 1 V
+5 0 V
+4 1 V
+4 0 V
+5 1 V
+4 0 V
+4 1 V
+5 1 V
+4 1 V
+4 0 V
+5 1 V
+4 1 V
+4 1 V
+5 1 V
+4 1 V
+4 2 V
+5 1 V
+4 1 V
+4 1 V
+5 2 V
+4 1 V
+4 2 V
+5 1 V
+4 2 V
+4 2 V
+5 1 V
+4 2 V
+4 2 V
+5 2 V
+4 2 V
+4 2 V
+5 2 V
+4 2 V
+4 2 V
+5 2 V
+4034 3421 L
+4 2 V
+5 3 V
+4 2 V
+4 2 V
+5 2 V
+4 3 V
+4 2 V
+5 2 V
+4 3 V
+4 2 V
+5 2 V
+4 3 V
+4 2 V
+5 2 V
+4 3 V
+4 2 V
+5 3 V
+4 2 V
+4 2 V
+5 3 V
+4 2 V
+4 2 V
+5 3 V
+4 2 V
+4 2 V
+5 2 V
+4 3 V
+4 2 V
+5 2 V
+4 2 V
+4 2 V
+5 2 V
+4 2 V
+4 2 V
+5 2 V
+4 2 V
+4 2 V
+5 2 V
+4 1 V
+4 2 V
+5 2 V
+4 1 V
+4 2 V
+5 2 V
+4 1 V
+4 2 V
+5 1 V
+4 1 V
+4 2 V
+5 1 V
+4 1 V
+4 1 V
+5 2 V
+4 1 V
+4 1 V
+5 1 V
+4 1 V
+4 1 V
+5 1 V
+4 0 V
+4 1 V
+5 1 V
+4 1 V
+4 1 V
+5 0 V
+4 1 V
+4 1 V
+5 0 V
+4 1 V
+4 0 V
+5 1 V
+4 0 V
+4 1 V
+5 0 V
+4 0 V
+4 1 V
+5 0 V
+4 0 V
+4 1 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 1 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4489 3545 L
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 0 V
+5 -1 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 0 V
+4 -1 V
+5 0 V
+4 -1 V
+4 0 V
+5 0 V
+4 -1 V
+1.000 UP
+stroke
+LTb
+stroke
+grestore
+end
+showpage
+%%Trailer
+%%DocumentFonts: Helvetica
+%%Pages: 1
diff --git a/faraday/psr_vs_detuning_combo.m b/faraday/psr_vs_detuning_combo.m
new file mode 100644
index 0000000..b2ea9d3
--- /dev/null
+++ b/faraday/psr_vs_detuning_combo.m
@@ -0,0 +1,149 @@
+1;
+
+data_dir='results/';
+output_dir='results/';
+N_detun_steps=1000;
+detuning_p_min=-200.0;
+%detuning_p_max=-detuning_p_min;
+detuning_p_max=1000;
+detuning_freq=linspace(detuning_p_min,detuning_p_max,N_detun_steps);
+
+gmg=.7; % gyro magnetic ration for ground level
+zeeman_splitting=+0.000;
+B_field=zeeman_splitting/gmg;
+
+%[psr_rad]=psr_vs_detuning(Ep, psi_el, B_field, theta, phi)
+
+% phi is angle between linear polarization and axis x
+phi=pi/4;
+% 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)
+psi_el=-30/180*pi;
+
+
+figure(6);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% zero magnetic field,, 30 degree ellipticity
+zeeman_splitting=+0.000;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, psi_el, theta, phi);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% 0.1 G magnetic field,, 30 degree ellipticity
+
+zeeman_splitting=+0.070;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% 0.0001 G magnetic field,, 30 degree ellipticity
+
+zeeman_splitting=+0.000070;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% 1.0 G magnetic field,, 30 degree ellipticity
+
+zeeman_splitting=+0.70;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% -0.1 G magnetic field,, 30 degree ellipticity
+
+zeeman_splitting=-0.070;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% -0.0001 G magnetic field,, 30 degree ellipticity
+
+zeeman_splitting=-0.000070;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+% -1.0 G magnetic field,, 30 degree ellipticity
+
+zeeman_splitting=-0.70;
+B_field=zeeman_splitting/gmg;
+psi_el=30/180*pi;
+
+%[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] =make_representative_psr_vs_detuning_for_given_B_and_psi_el(detuning_freq, B_field, theta, phi, psi_el);
+fname=data_file_name('results/', 'PSR.','mat', B_field, theta,phi,psi_el);
+load(fname)
+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 );
diff --git a/faraday/psr_vs_power.m b/faraday/psr_vs_power.m
new file mode 100644
index 0000000..97b9fab
--- /dev/null
+++ b/faraday/psr_vs_power.m
@@ -0,0 +1,147 @@
+1;
+clear all;
+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;
+%four_levels_with_polarization;
+%four_levels;
+%three_levels;
+%two_levels;
+
+% load EM field description
+field_description;
+
+%Nfreq=length(modulation_freq);
+
+
+
+%tune probe frequency
+detuning_p=0;
+N_detun_steps=100;
+%detuning_p_min=-B_field*gmg*4; % span +/-4 Zeeman splitting
+detuning_p_min=-200.0;
+detuning_p_max=-detuning_p_min;
+detuning_p_max=1000;
+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;
+
+fprintf (stderr, "calculating atom properties\n");
+fflush (stderr);
+pfile='rb87_D1_line.m'; % the parent file from which L0_and_polarization_submatrices calculated
+cfile='L0m_and_polarizability_calculated.mat'; % the child file to which calculated matrices writen
+[s, err, msg] = stat (pfile);
+if(err)
+ %file does not exist
+ disp('Big troubles are coming, no file to define Hamiltonian)');
+ msg=cstrcat('File: ', pfile, ' is missing...exiting');
+ disp(msg);
+ return;
+else
+ pfile_mtime=s.mtime;
+endif
+[s, err, msg] = stat (cfile);
+if(err)
+ %file does not exist
+ cfile_mtime=0;
+else
+ cfile_mtime=s.mtime;
+endif;
+if ( cfile_mtime >= pfile_mtime)
+ % matrices already calculated and up to date, all we need to load them
+ load(cfile);
+ else
+ % 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 ...
+ );
+ save(cfile, 'L0m', 'polarizability_m');
+ endif
+elapsed_time = etime (clock (), t0);
+fprintf (stderr, "elapsed time so far is %.3f sec\n",elapsed_time);
+fflush (stderr);
+
+global atom_properties;
+atom_properties.L0m=L0m;
+atom_properties.polarizability_m=polarizability_m;
+atom_properties.dipole_elements=dipole_elements;
+
+
+% phi is angle between linear polarization and axis x
+phi=pi*2/8;
+% 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)
+psi_el=-5/180*pi;
+
+
+
+fprintf (stderr, "tuning laser in forloop to set conditions vs detuning\n");
+fflush (stderr);
+wp=w_pf1-w_sf2 +80; %Fg=2 -> Fe=1 +80 MHz
+Ep=logspace(-2,1,100);
+for cntr=1:length(Ep);
+
+ %light_positive_freq = [wp];
+ E_field_drive = [0 ];
+ E_field_probe = [Ep(cntr) ];
+ 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);
+
+
+ light_positive_freq=[ wp];
+ % we calculate dc and negative frequiencies 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{cntr}=atom_field_problem;
+
+endfor
+
+save '/tmp/problem_definition.mat' problems_cell_array atom_properties Ep ;
+fprintf (stderr, "now really hard calculations begin\n");
+fflush (stderr);
+% once we define all problems the main job is done here
+[xi_linear, xi_left, xi_right]=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
+%[xi_linear, xi_left, xi_right]=cellfun( @susceptibility_steady_state_at_freq, problems_cell_array);
+
+%save '/tmp/relative_transmission_vs_detuning.mat' detuning_freq relative_transmission_vs_detuning;
+save '/tmp/xi_vs_power.mat' Ep xi_linear xi_left xi_right E_field_pos_freq wp;
+
+%output_psr_results_vs_detuning;
+output_psr_results_vs_power;
+
+elapsed_time = etime (clock (), t0)
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/faraday/rb87_D1_line.m b/faraday/rb87_D1_line.m
new file mode 100644
index 0000000..74a52f0
--- /dev/null
+++ b/faraday/rb87_D1_line.m
@@ -0,0 +1,148 @@
+1;
+useful_constants;
+
+% note all frequency values are in MHz
+
+% 87Rb D1 line
+%
+% m=-2 m=-1 m=0 m=1 m=2
+% ---- ---- ---- ---- ---- |P,F=2>
+%
+% ---- ---- ---- |P,F=1>
+% m=-1 m=0 m=1
+%
+%
+%
+%
+%
+% m=-2 m=-1 m=0 m=1 m=2
+% ---- ---- ---- ---- ---- |S,F=2>
+%
+% ---- ---- ---- |S,F=1>
+% m=-1 m=0 m=1
+
+
+w_hpf_ground=6834;
+w_hpf_exited=817;
+w_sf2 = w_hpf_ground; % Distance from |S,F=1> to |S,F=2>
+w_pf1 =1e9; % something big Distance from |S,F=1> to |P,F=1>
+w_pf2 = w_pf1+w_hpf_exited; %Distance from |S,F=1> to |P,F=2>
+gmg=.7; % gyro magnetic ration for ground level
+gme=.23; % gyro magnetic ration for exited level % CHECKME
+
+
+%bottom level |F=1>
+levels( 1)=struct( "ang_momentum", 0, "total_momentum", 1, "m", -1, "energy", 0, "gm", -gmg);
+levels( 2)=struct( "ang_momentum", 0, "total_momentum", 1, "m", 0, "energy", 0, "gm", -gmg);
+levels( 3)=struct( "ang_momentum", 0, "total_momentum", 1, "m", 1, "energy", 0, "gm", -gmg);
+
+%second bottom level |F=2>
+levels( 4)=struct( "ang_momentum", 0, "total_momentum", 2, "m", -2, "energy", w_sf2, "gm", gmg);
+levels( 5)=struct( "ang_momentum", 0, "total_momentum", 2, "m", -1, "energy", w_sf2, "gm", gmg);
+levels( 6)=struct( "ang_momentum", 0, "total_momentum", 2, "m", 0, "energy", w_sf2, "gm", gmg);
+levels( 7)=struct( "ang_momentum", 0, "total_momentum", 2, "m", 1, "energy", w_sf2, "gm", gmg);
+levels( 8)=struct( "ang_momentum", 0, "total_momentum", 2, "m", 2, "energy", w_sf2, "gm", gmg);
+
+% first exited level |F=1>
+levels( 9)=struct( "ang_momentum", 1, "total_momentum", 1, "m", -1, "energy", w_pf1, "gm", -gme);
+levels(10)=struct( "ang_momentum", 1, "total_momentum", 1, "m", 0, "energy", w_pf1, "gm", -gme);
+levels(11)=struct( "ang_momentum", 1, "total_momentum", 1, "m", 1, "energy", w_pf1, "gm", -gme);
+
+% second exited level |F=2>
+levels(12)=struct( "ang_momentum", 1, "total_momentum", 2, "m", -2, "energy", w_pf2, "gm", gme);
+levels(13)=struct( "ang_momentum", 1, "total_momentum", 2, "m", -1, "energy", w_pf2, "gm", gme);
+levels(14)=struct( "ang_momentum", 1, "total_momentum", 2, "m", 0, "energy", w_pf2, "gm", gme);
+levels(15)=struct( "ang_momentum", 1, "total_momentum", 2, "m", 1, "energy", w_pf2, "gm", gme);
+levels(16)=struct( "ang_momentum", 1, "total_momentum", 2, "m", 2, "energy", w_pf2, "gm", gme);
+
+Nlevels=size(levels)(2);
+
+
+H0=zeros(Nlevels);
+energy = [1:Nlevels].*0;
+ang_momentum = [1:Nlevels].*0;
+total_momentum = [1:Nlevels].*0;
+m = [1:Nlevels].*0;
+gm = [1:Nlevels].*0;
+
+for i=1:Nlevels
+ energy(i) = levels(i).energy;
+ ang_momentum(i) = levels(i).ang_momentum;
+ total_momentum(i) = levels(i).total_momentum;
+ m(i) = levels(i).m;
+ gm(i) = levels(i).gm;
+endfor
+H0=diag(energy)*hbar;
+
+
+dipole_elements.left = zeros(Nlevels);
+dipole_elements.right = zeros(Nlevels);
+dipole_elements.linear = zeros(Nlevels);
+% define dipole elements for transition from level j->k
+for j=1:Nlevels
+ for k=1:Nlevels
+ if ( abs(ang_momentum(j) - ang_momentum(k)) == 1)
+ %transition allowed for L =L' +/- 1
+ % but they correct within a factor, but not sign
+ if ( ( H0(j,j) < H0(k,k)) )
+ de=dipole_elementRb87D1line(...
+ total_momentum(j),m(j), ...
+ total_momentum(k),m(k) ...
+ );
+ dipole_elements.left(j,k) = de.left;
+ dipole_elements.right(j,k) = de.right;
+ dipole_elements.linear(j,k) = de.linear;
+ endif
+ endif
+ endfor
+endfor
+dipole_elements.linear+=dipole_elements.linear';
+dipole_elements.left+=dipole_elements.left';
+dipole_elements.right+=dipole_elements.right';
+
+maximum_dipole_elements=abs(dipole_elements.linear) + abs(dipole_elements.left) + abs(dipole_elements.right);
+
+%defasing matrix
+g_deph=0;
+g_dephasing=zeros(Nlevels);
+% dephasing only for non zero dipole elements (am I right?)
+g_dephasing=g_deph*(abs(maximum_dipole_elements) != 0);
+
+% decay matrix g(i,j) correspnds to decay from i-->j
+gamma=6;
+g_decay=zeros(Nlevels);
+for i=1:Nlevels
+ for j=1:Nlevels
+ if ( H0(i,i) > H0(j,j) )
+ % only upper levels decaying and decay is always positive
+ g_decay(i,j)=gamma * abs( maximum_dipole_elements(i,j) );
+ endif
+ endfor
+endfor
+
+% ground level mixing need to be artificial
+gamma_hpf=.0001;
+for i=1:Nlevels
+ for j=1:Nlevels
+ % it would be better to introduce a level corresponding to the bath
+ % but this slows calculations
+ % So
+ % ground hyperfine are equally mixed together
+ if ( (abs( H0(i,i) - H0(j,j)) == w_hpf_ground*hbar ) || ...
+ (abs( H0(i,i) - H0(j,j)) == 0 ) )
+ % i and j correspond to 2 ground levels
+ % we cannot decay to itself
+ if ( i != j )
+ % total eight ground levels F=2 and F=1 so we decay in 7 channels
+ g_decay(i,j)+=gamma_hpf/7;
+ endif
+ endif
+ endfor
+endfor
+
+% apply B field Zeeman splitting
+energy+=B_field * ( gm .* m);
+% convert frequency to energy units
+H0=diag(energy)*hbar;
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/faraday/useful_constants.m b/faraday/useful_constants.m
new file mode 100644
index 0000000..82d1d1a
--- /dev/null
+++ b/faraday/useful_constants.m
@@ -0,0 +1,6 @@
+1;
+
+% some physical constants
+hbar=1;
+im_one=0+1i;
+el_charge=1;
diff --git a/faraday/useful_functions.m b/faraday/useful_functions.m
new file mode 100644
index 0000000..eb3e880
--- /dev/null
+++ b/faraday/useful_functions.m
@@ -0,0 +1,337 @@
+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
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 11:26:36 +0000