copy psr to faraday

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
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