Calculations for the Sodium FADOF from FITS files
The notebook calculates the transmission spectra and other details of the sodium FADOF. Since the original calculations take quite some time, it makes use of precalculated FITS files.

Copyright (C) 2011, 2014 Joanna A. Zielinska, Morgan W. Mitchell, Wilhelm Kiefer, Ilja Gerhardt
This program is free software : you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.You should have received a copy of the GNU General Public License along with this program.If not, see <http://www.gnu.org/licenses/> .

We would appreciate all publications using the output of this program to contain a brief acknowledgement “The authors acknowledge the courtesy of Joanna A. Zielinska, Morgan W. Mitchell (ICFO, Barcelona) and Wilhelm Kiefer, Ilja Gerhardt (MPI, Stuttgart) for permission to use this program”. Thank you.

This notebook can be compiled by pressing ”Ctrl - A” and then ”Shift - Enter”

SetOptions[{ListPlot},{Joined->True,Axes->False,Frame->True,FrameStyle->Thick,PlotStyle->{Thick,Black},PlotRange->All,Background->None,ImageSize->{Large},BaseStyle->{FontFamily->"Arial",FontSize->16},GridLines->Automatic,GridLinesStyle->Directive[{Red,Dashed}],Filling->0}];

First, some definitions :

line="d1"; (* which line should be plotted *)

temp=273+50; (* only values from 100degC to 250degC are allowed *)

length=0.1; (* the length of the cell in meters *)

size=""; (* ideally you take the small ('_small') data first *)

importurl="http://localhost/cesium/"; (* I keep the localhost in here, ppl using this program should know what to do *)

Make everything pretty

si1=8;si2=16;si3=18;fo="Arial"; (* Different font sizes *)

SetOptions[ListDensityPlot,{ColorFunction->"AvocadoColors",InterpolationOrder->0,FrameStyle->Thick,FrameTicks->{{yti,Automatic},{ftix,None}},PlotRange->{0,1},ColorFunctionScaling->False,InterpolationOrder->0,FrameStyle->Thick,ImageSize->600,BaseStyle->{FontFamily->fo,FontSize->si2},ImageSize->400,FrameLabel->{"Laser Detuning [GHz]","Magnetic Field [mT]"},FrameLabel->{"Laser Detuning [GHz]","Magnetic Field [mT]"}}];

Function definitions :

c=3 (10^8); (* speed of light *)

energyd1=335116048.807 ;(* cesium D1 transition in units of MHz, Steck, Table 4, Line 1 *)

energyd2=351725718.50 ;(* cesium D2 transition in units of MHz, Steck, Table 3, Line 1 *)

EIn={1/Sqrt[2],1/Sqrt[2]}; (* vectors in and outgoing light *)
EDet={1/Sqrt[2],-1/Sqrt[2]};

GenMatrix[data_]:= (* let's generate all matrices *)
Table[
thisf=freqtab[[ff]]+ToExpression["energy"<>line]; (* get the base frequency of the present line *)
nR=10^6Sqrt[1+10^(-6)(data[[1,1,bf,ff]]+I data[[1,3,bf,ff]])];
nL=10^6Sqrt[1+10^(-6)(data[[1,2,bf,ff]]+I data[[1,4,bf,ff]])];
M={{Exp[I nR length 2π *thisf/c],0},{0,Exp[I nL length 2π *thisf/c]} };
M
,{ff,1,Dimensions[data][[-1]]} (* all frequencies *)
,{bf,1,Dimensions[data][[-2]]}]; (* all bfields *)

calvec[vec_]:=(Abs[vec[[1]]]^2+Abs[vec[[2]]]^2) {Abs[vec[[1]]]^2-Abs[vec[[2]]]^2,2 Re[Conjugate[vec[[1]]] vec[[2]]],2 Im[Conjugate[vec[[1]]] vec[[2]]]} (* calculate from jones -> stokes *)

calvecr[vec_]:={Abs[vec[[1]]]^2-Abs[vec[[2]]]^2,2 Re[Conjugate[vec[[1]]] vec[[2]]],2 Im[Conjugate[vec[[1]]] vec[[2]]]} (* same, but w/o length *)

lhp={1,0};lvp={0,1};lplus=1/Sqrt[2] {1,1};lminus=1/Sqrt[2] {1,-1};rcp=1/Sqrt[2] {1,I};lcp=1/Sqrt[2] {1,-I}; (* a variety of input states *)

co:=Map[calvec,vout];

winki2:=Map[VectorAngle[calvec[lminus],#]&,co];

cs:=Graphics3D[{Lighter[Gray],Line[Table[{Sin[xx],Cos[xx],0},{xx,-.1,2π,0.01}]],Line[Table[{0,Sin[xx],Cos[xx]},{xx,-.1,2π,0.01}]],Line[Table[{Sin[xx],0,Cos[xx]},{xx,-0.1,2π,0.01}]]}];

ar[vec_,col_: Red]:=Graphics3D[{col,Thick,Arrow[{{0,0,0},calvec[vec]}]}]; (* draw an arrow to the resulting polarization *)

coords:=Graphics3D[{Line[{{-1.2,0,0},{1.2,0,0}}],Line[{{0,-1.2,0},{0,1.2,0}}],Line[{{0,0,-1.2},{0,0,1.2}}]}]; (* put some lines in there, xyz coordinates *)

t:=Graphics3D[{Blue,Thick,Line[threed]}];

threed:=Map[Normalize,Map[calvecr,vout]]; (* we want to have the normalized vector here *)

Import the data :

data=Import[importurl<>line<>"_data/"<>line<>"_calc_"<>ToString[temp]<>size<>".fits","RawData"]; (* let's fetch the data from the website *)

freqtab=Table[xx,{xx,-20000,20000,(40000/(Dimensions[data][[-1]]-1))}]; (* calculated from -20 to +20 GHz *)

Calculate one spectrum

g=GenMatrix[data];

Assume one BField only

bfield=40; (* this is the line which shows 40mT *)

sf=0.2; (* 5 lines = 1 mT *)

bfpos=Round[1/sf bfield];

ListPlot[
{
Table[{freqtab[[xx]]/1000,Abs[EDet.g[[xx,bfpos]].EIn]^2},{xx,1,Length[g]}], (* FADOF Transmission *)
Table[{freqtab[[xx]]/1000,Abs[EIn.g[[xx,bfpos]].EIn]^2},{xx,1,Length[g]}], (* X FADOF Transmission *)
Table[{freqtab[[xx]]/1000,Abs[EDet.g[[xx,bfpos]].EIn]^2+Abs[EIn.g[[xx,bfpos]].EIn]^2},{xx,1,Length[g]}], (* Doppler Transmission *)
Table[{freqtab[[xx]]/1000,((Abs[g[[xx,bfpos,1]].{1/Sqrt[2],1/Sqrt[2]}])^2)},{xx,1,Length[g]}],(* transmission for sigma (-) *)
Table[{freqtab[[xx]]/1000,((Abs[g[[xx,bfpos,2]].{1/Sqrt[2],1/Sqrt[2]}])^2)},{xx,1,Length[g]}](* transmission for sigma (+) *)
},PlotStyle->{Thick,Thick,Thick,Thin,Thin},Joined-> True,FrameLabel-> {"Frequncy Detuning [GHz]","Transmission"}]

vout:=g[[All,bfpos]].lplus; (* let the matrix act on a certain polarization (lplus, assumed) *)

ListPlot[Transpose[{freqtab/1000,(Cos[winki2]+1)/2}],Joined-> True,FrameLabel-> {"Frequncy Detuning [GHz]","Transmission"}] (* this is the transmission according to the rotation, but no (Doppler-) abs assumed *)

ups=Show[{cs,coords,t,ar[lplus],ar[lhp]},Boxed-> False,ViewPoint->{0.9,-1.1,0.6}] (* how does the vector change, when I scan the magnetic field, plot as stokes vector *)

Calculate one spectral line

detuning=4200;

pos=Round[(detuning+20000)/(40000/Dimensions[data][[-1]]),1]+1;(*this is the line we need for this particular detuning*)

ListPlot[
{
Table[{sf xx,Abs[EDet.g[[pos,xx]].EIn]^2},{xx,1,Length[g[[pos]]]}], (* FADOF Transmission *)
Table[{sf xx,Abs[EIn.g[[pos,xx]].EIn]^2},{xx,1,Length[g[[pos]]]}], (* X FADOF Transmission *)
Table[{sf xx,Abs[EDet.g[[pos,xx]].EIn]^2+Abs[EIn.g[[pos,xx]].EIn]^2},{xx,1,Length[g[[pos]]]}], (* Doppler Transmission *)
Table[{sf xx,((Abs[g[[pos,xx,1]].{1/Sqrt[2],1/Sqrt[2]}])^2)},{xx,1,Length[g[[pos]]]}],(* transmission für sigma (-) *)
Table[{sf xx,((Abs[g[[pos,xx,2]].{1/Sqrt[2],1/Sqrt[2]}])^2)},{xx,1,Length[g[[pos]]]}](* transmission für sigma (+) *)
},PlotStyle->{Thick,Thick,Thick,Thin,Thin},Joined-> True,FrameLabel-> {"Magnetic Field [mT]","Transmission"}]

vout=g[[pos,All]].lplus;

ListPlot[(Cos[winki2]+1)/2,Joined-> True,FrameLabel-> {"Magnetic Field [mT]","Transmission"}]

ups=Show[{cs,coords,t,ar[lplus],ar[lhp]},Boxed-> False,ViewPoint->{0.9,-1.1,0.6}] (* how does the vector change, when I scan the magnetic field, plot as stokes vector *)

Make 2 D Plot

ftix=Table[{N[(xx-Fit[freqtab,{1,x},x][[1]])/Fit[freqtab,{1,x},x][[2,1]]],StringReplace[ToString[N[xx/1000]],RegularExpression["\\.$"]->""]},{xx,-20000,20000,5000}];

xpo=1900;ypo=190; (* position of the legend *)

rang=10;absmax=100;absmin=0; (* the nice legend of the transmission plots *)l1=Show[Graphics[{White,Rectangle[{-0.1,-0.1},{1.1,rang+1.1}],Flatten[Table[{ColorData["AvocadoColors"][x],Rectangle[{0,0+rang x},{1,1+rang x}],ColorData["AvocadoColors"][1-x],Text[StringReplace[ToString[Round[(absmax-absmin)/rang IntegerPart[rang x],0.01]],"."->""],{0.5,rang x+0.5},{0,0},BaseStyle->{FontFamily->fo,FontSize->si1}]},{x,0,1,1/rang}]]}]];

fadofout=Table[Abs[EDet.g[[x,y]].EIn]^2,{y,1,Length[g[[1]]]},{x,1,Length[g]}]; (* calculate the FADOF spectrum for all values (2D Plot) *)

xfadofout=Table[Abs[EIn.g[[x,y]].EIn]^2,{y,1,Length[g[[1]]]},{x,1,Length[g]}];(* calculate the Cross FADOF spectrum for all values (2D Plot) *)

yti=Table[{xx/sf,xx},{xx,0,100,20}]; (* the magnetic field scaling *)

aplot=ListDensityPlot[fadofout,PlotLabel->"FADOF Transmission"];
aaplot=Show[{aplot,Graphics[{Inset[l1,{xpo,ypo}]}]}];

bplot=ListDensityPlot[xfadofout,PlotLabel->"Cross FADOF Transmission"];
bbplot=Show[{bplot,Graphics[{Inset[l1,{xpo,ypo}]}]}];

cplot=ListDensityPlot[fadofout+xfadofout,PlotLabel->"Doppler Transmission"]; (* the Doppler transmission comes out, when you have both output ports *)
ccplot=Show[{cplot,Graphics[{Inset[l1,{xpo,ypo}]}]}];

opl=Show[GraphicsGrid[{{aaplot,bbplot,ccplot}}],PlotLabel->"Sodium "<>ToUpperCase[line]<>" transition at "<>StringReplace[ToString[temp-273],"."->""]<>"°C, "<>StringReplace[ToString[1000 length],"."->""]<>" mm length",BaseStyle->{FontFamily->fo,FontSize->si3}]