import numpy as np 

print('Vector light shift induced by tight focusing, parameters from Lukin PRL')
lam_trap = 815e-9
print("FORT wavelength %i nm" % (lam_trap*1e9))
c = 3e8
d1 = c/lam_trap - c/794.98e-9 
d2 = c/lam_trap - c/780.24e-9 
kb=1.38e-23
h=6.626e-34
hbar = 1e-34
U_trapdepth = 0.82e-3 # in [K]
print("trap depth %2.2f mK" % (U_trapdepth*1e3))
U0 = U_trapdepth*kb/h # trap depth in lin frequency

alpha = 0.43
print("numerical aperture: %2.2f " % (alpha))
dCdy = 2.6 * alpha *np.sin(alpha) / lam_trap # = 0.1 [1/um]
print("-> dCdy =  %2.2f [1/um]. (in paper: 0.57 [1/um] )" % (dCdy*1e-6))

muB = 9.274e-24 #Bohr magneton 
Zeeman = muB/h *0.5 # [Hz/Tesla] = 0.7 [MHz/Gauss] Steck RB 87
dUdy = U0 * dCdy * (d2 - d1)/(2*d1+d2) * 0.5  / Zeeman # eq 1 of lukin paper; units [Tesla/m] 
print("-> dUdy =  %2.2f [Gauss/um]. (in paper: 1.4 [Gauss/um] )" % (dUdy*1e-6*1e4))

trap_freq = 2*np.pi*100e3 # radial trap freq
u = 1.66e-27
mass = 87*u
B_z = 10.5e-4 # = 10.5 Gauss
eta_pg = 3/2*muB * (dUdy)**2/(B_z *3*mass*trap_freq**2) #
print("-> eta_pg =  %2.4f [1/s]" % (eta_pg))
T=40e-6
T2 = 0.97 * 2*hbar / (kb*T*eta_pg)
print("-> T2 =  %2.5f [s]" % (T2))
print("-> 1 / T2 =  %2.3f [1/ms]. (in fig 1d of paper about 4 [1/ms]) \n" % (1/T2*1e-3))


print('Vector light shift induced by tight focusing, our parameters')
lam_trap = 850e-9
print("FORT wavelength %i nm" % (lam_trap*1e9))
c = 3e8
d1 = c/lam_trap - c/794.98e-9 
d2 = c/lam_trap - c/780.24e-9 
kb=1.38e-23
h=6.626e-34
hbar = 1e-34
U_trapdepth = 0.88e-3 # in [K]
print("trap depth %2.2f mK" % (U_trapdepth*1e3))
U0 = U_trapdepth*kb/h # trap depth in lin frequency

alpha = 0.20 # effective NA for trap beam
# NA = sin (divergence angle) -> NA = sin (lam/(np.pi*w0)) #wo = 1.4e-6 beam waist
# -> NA = alpha = 0.2
# alternative: using gaussian beam waist after lens in par-axial approx -> to get 1.4um waist one needs F/w1 = 5.208 -> NA = 0.2

# alternative scale NA from lukin paper: alpha = 0.43 / 1.5 = 0.28
print("numerical aperture: %2.2f " % (alpha))
dCdy = 2.6 * alpha *np.sin(alpha) / lam_trap # = 0.1 [1/um]
print("-> dCdy =  %2.2f [1/um]. " % (dCdy*1e-6))

muB = 9.274e-24 #Bohr magneton 
Zeeman = muB/h *0.5 # [Hz/Tesla] = 0.7 [MHz/Gauss] Steck RB 87
dUdy = U0 * dCdy * (d2 - d1)/(2*d1+d2) * 0.5  / Zeeman # eq 1 of lukin paper; units [Tesla/m] 
print("-> dUdy =  %2.2f [Gauss/um]. " % (dUdy*1e-6*1e4))

trap_freq = 2*np.pi*100e3 # radial trap freq
u = 1.66e-27
mass = 87*u
B_z = 14.4e-4 # = 14.4 Gauss
eta_pg = 3/2*muB * (dUdy)**2/(B_z *3*mass*trap_freq**2) #
print("-> eta_pg =  %2.5f [1/s]" % (eta_pg))
T=40e-6
T2 = 0.97 * 2*hbar / (kb*T*eta_pg)
print("-> T2 =  %2.5f [s]" % (T2))
print("-> 1 / T2 =  %2.3f [1/ms]. \n" % (1/T2*1e-3))


#waist = 2*lam/np.pi * F/w0
#NA = lam/(np.pi*w0) #wo = beam waist