<<PhysicalConstants`

m=SI[7 AMU]/Kilogram;
μB=SI[ BohrMagneton]/(Joule/Tesla);
μbar=((μB)/ℏ)/(Hz/Tesla); (*1.4  MHz/Gauss*)
c=SpeedOfLight/(Meter/Second);
kb=BoltzmannConstant/(Joule/Kelvin);
ℏ=PlanckConstantReduced /(Joule Second);
h=PlanckConstant/(Joule Second);
g=AccelerationDueToGravity/(Meter/(Second^2));
λ=670.9  (10^(-9)) ;
kL=((2 π)/λ);
w0=2π c/λ;
Γ=(2π)5.9 (10^6) ;
 i=(3/2);s=(1/2); l=0; j=(1/2);
gj=(3/2)+((s(s+1)-l(l+1))/(2j(j+1)));
gf[f_]:=gj ((f(f+1)-i(i+1)+j(j+1))/(2f(f+1)));
μ[f_,mf_]:=gf[f] mf μB;
vth[T_]:=(Sqrt[((2 kb T ((10^(-6))))/m)])(* [T]=μK*)

Isat=((ℏ (w0^3) )/(12 π (c^2)))Γ;(*[Is]= W/(m^2)*)
Isatm=((ℏ (w0^3) )/(12 π (c^2)))Γ (((10^3)/(10^4)));(*[Is]= mW/(cm^2)*)
σ0=(3/(2 π))(λ^2);
aB=5.3 (10^(-11));(*Bohr radius 5.3 Angst*)
as=-27.6aB; (*in (2,2) *)
σs=8 π (as^2); (*scattering cross section in (m^2)*)
σscm=σs (10^4);
λdB[T_]:=(Sqrt[((2 π (ℏ^2))/(m kb T))]); (*de Broglie wavelength in meters*)
C6=1.328 (10^(-76)); (*Van Der Waals for 7Li in (2,2) in units J/(m^6)*)
ac=(((((m C6)/(ℏ^2))))^(1/4)); (*range of Van der Waals potential m*)

Gf[x_,y_,z_,σx0_,σy0_,σz0_,T_,t_]:=(1/((((2 π))^(3/2))σx[σx0,T,t]σy[σy0,T,t] σz[σy0,T,t]))Exp[-(((((x (10^(-2))))^2)/(2 ((σx[σx0,T,t])^2))))-(((((y (10^(-2))))^2)/(2 ((σy[σy0,T,t])^2))))-(((((z (10^(-2))))^2)/(2 (((σz[σz0,T,t]/(Sqrt[2])))^2))))]; (*Gaussian cloud, x, y, z in cm *)
σx[T_,σx0_,t_]:=(Sqrt[(((σx0 (10^(-2))))^2)+(((vth[T](t (10^(-3)))))^2)]);(* t in ms*)
σy[T_,σy0_,t_]:=(Sqrt[(((σy0 (10^(-2))))^2)+(((vth[T](t (10^(-3)))))^2)]);(* t in ms*)
σz[T_,σz0_,t_]:=(Sqrt[(((σz0 (10^(-2))))^2)+(((vth[T](t (10^(-3)))))^2)]);(* t in ms*)
Manipulate[ContourPlot[{Gf[x,1,z,1,1,1,1000,t]},{z,-5,5},{x,-5,5},PlotRange->Full],{t,0,3}]

ClearAll[P,rc]
P=25; (*in mW per arm*)
rc=2.5;(*beam size in cm*)
Int=(P/(π (rc^2)));
σ0=1 (10^(-2)); (*1cm in m*)
α[δ_]:=2ℏ (kL^2)(((2 (Int/Isatm))(((2 δ)/Γ)))/(((1+(((((2 δ)/Γ)))^2)))^2));(* Molasses damping constant, δ in MHz*)
σ[δ_,T_,t__]:=(Sqrt[((σ0)^2)+2((kb T (10^(-3)))/(α[δ ]))t (10^(-3))])(* in m*)

GfM[x_,y_,z_,δ_,T_,t_]:=(1/((((2 π))^(3/2))σ[δ,T,t]σ[δ,T,t] σ[δ,T,t]))*Exp[-(((((x (10^(-2))))^2)/(2 ((σ[δ,T,t])^2))))-(((((y (10^(-2))))^2)/(2 ((σ[δ,T,t])^2))))-(((((z (10^(-2))))^2)/(2 (((σ[δ,T,t]/(Sqrt[2])))^2))))]; 
Animate[ContourPlot[{GfM[x,0,z,2 π 40,1,t]},{z,-5,5},{x,-5,5},PlotRange->Full],{t,0,3}]

GfM[1,1,1,2 π 40,1,0]
Gf[1,1,1,1,1,1,1,0]