function psr_rad=output_faraday_results_vs_B()
load '/tmp/xi_vs_B.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(B_fields, real(xi_left-xi_right), '-');  
title("differential real xi");
xlabel("two photon detuning (MHz)");

figure(2);
hold off;
plot(B_fields, imag(xi_left-xi_right), '-');  
title("differential imag xi");
xlabel("two photon detuning (MHz)");

figure(3);
hold off;
plot(B_fields, real(xi_left), '-', detuning_freq, real(xi_right), '-');  
title("real xi");
xlabel("two photon detuning (MHz)");

figure(4);
hold off;
plot(B_fields, 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(B_fields, 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, detuning_freq); 
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_B_Fg=2toFe=1,2_Ip=%.3f_el_%.1f_detun=%.3fMHz.mat", I_probe, psi_el*180/pi,detuning_freq); 
save(fname,'B_fields', 'psr_rad', 'xi_left', 'xi_right', 'xi_linear');

return;