r"""Various Optical Traps.
"""
from __future__ import division
import numpy as np
from numpy import pi
from mmf.objects import StateVars, Container, process_vars
from mmf.objects import Computed, Delegate
[docs]class Units(Container):
"""Unit conversions.
>>> u = Units()
>>> kB = 1.3806505e-23*u.m*u.m*u.kg/u.s/u.s/u.K
>>> kB
1.0
>>> hbar = 1.05457168e-34*u.m*u.m*u.kg/u.s
>>> hbar
1.0
>>> Li_m = 6.0151214*u.amu
>>> Li_m
1.0
>>> aosc = np.sqrt(hbar/Li_m/(2.*pi*325./u.s))
>>> aosc
1.0
"""
_state_vars = [
('s', 2042.0352248333654, "1 second in natural units"),
('m', 439785.04589706744, "1 meter in natural units"),
('kg', 1.0011662693416561e+26, "1 kg in natural units"),
('K', 64112752.248760417, "1 K in natural units"),
('amu', 1./6.0151214, "1 amu in natural units"),
('Hz', Computed),
('micron', Computed),
('mm', Computed),
('nK', Computed),
('J', Computed),
('hbar', Computed),
('Li6', Computed, "Mass of Lithium 6 species"),
('aosc', Computed, "Oscillator length for Li6"),
]
process_vars()
def __init__(self, *v, **kw):
# Some derived units.
self.Hz = 1./self.s
self.micron = 1e-6*self.m
self.mm = 1e-3*self.m
self.nK = 1e-9*self.K
self.J = self.kg*(self.m/self.s)**2
self.Li6 = Container(m=6.0151214*self.amu)
self.hbar = 1.05457168e-34*self.m**2*self.kg/self.s
self.kB = 1.3806505e-23*self.m**2*self.kg/self.m**2/self.K
self.aosc = np.sqrt(self.hbar/self.Li6.m/(2*pi*325/self.s))
# Here are the real units: use the unum.unit package to check dimensions.
UNITS = Units()
[docs]class Hulet(StateVars):
r"""Parameters for Rice trapping potential: cont-mat/0908455
>>> t = Hulet()
>>> t.U(0, 0)
0.0
"""
_state_vars = [
('units', Delegate(Units, ['Li6'])),
('omega_B', Computed),
('w_0', Computed),
('z_0', Computed),
('m=units.Li6.m'),
('U_0', Computed),
]
process_vars()
def __init__(self, *v, **kw):
self.omega_B = 2*pi*6.5*self.units.Hz
self.w_0 = 26*self.units.micron
self.z_0 = 1.7*self.units.mm
self.U_0 = 540*self.units.nK
def _w(self, z):
return self.w_0*(1.0-(z/self.z_0)**2)
def U(self, r, z):
r"""Return the trapping Potential.
Parameters
----------
r, z : float
Radius and height in trap.
"""
return self.m*(self.omega_B*z)**2+ \
self.U_0*(1.0 - (self.w_0/self._w(z))**2* \
np.exp(-2.0*(r/self._w(z))**2))
[docs]class Model:
"""This class represents a model for the thermodynamics.
The compute() method must be called if the parameters are changed.
"""
param = Container(xi=0.42,
y0=-0.5,
y1=-0.05,
units=UNITS)
def __init__(self):
"""Calls compute()"""
self.compute()
[docs] def compute(self):
"""Recomputes values for speed."""
self.beta =\
(2*self.param.units.m/self.param.units.hbar**2)**(3/2)/\
(6*pi**2)
self.h0 = 1
self.h1 = (1 + self.param.y1)/(2*self.param.xi)**(3/5)
[docs] def h(self,y):
"""The function h(y) must be provided.
Here we provide a simple version interpolating between y0 and
y1 linearly.
>>> m1 = Model()
>>> y = np.linspace(-1.0,1.0,100)
>>> h = np.vectorize(m1.h)(y)
>> p1 = plot(y,h)
>>> yc = m1.special_y()
>>> hc = np.vectorize(m1.h)(yc)
>> p2 = plot(yc,hc,'x')
"""
xi,y0,y1 = self.param.xi,self.param.y0,self.param.y1
if y <= y0:
return self.h0
elif y <= y1:
return self.h1*(y-y0)/(y1-y0)+self.h0*(y1-y)/(y1-y0)
else:
return (1.0+y)/(2.0*xi)**(3.0/5.0)
[docs] def dh(self, y):
"""The function dh(y) must also be provided. If h(y) has any
kinks, then one of the limiting values should be provided and
the kink locations returned by special_y."""
xi, y0, y1 = self.param.xi, self.param.y0, self.param.y1
if y <= self.param.y0:
return 0.0
elif y <= self.param.y1:
return (self.h1 - self.h0)/(y1 - y0)
else:
return 1.0/(2.0*self.param.xi)**(3.0/5.0)
[docs] def special_y(self):
"""Return an array of special y's where h(y) has a kink."""
return np.array([self.param.y0,self.param.y1])
[docs] def nab(self,mua,mub):
"""Return the density pair (na,nb) for given chemical potentials."""
if mua <= 0 and mub <= 0:
return (0.0,0.0)
elif mua < mub:
(nb,na) = self.nab(mub,mua)
return (na,nb)
else:
y = mub/mua
h = self.h(y)
dh = self.dh(y)
tmp = self.beta*(mua*h)**(3./2.)
return (tmp*(h-y*dh),tmp*dh)
[docs]def mupm2muab((mup, mum)):
"""Return (mua,mub) from (mup,mum)."""
return (mup + mum, mup - mum)
[docs]def muab2mupm((mua,mub)):
"""Return (mup,mum) from (mua,mub)."""
return ((mua + mub)/2.,(mua - mub)/2.)
[docs]def mup_eff(x, y, z, trap=None):
"""Return the effective chemical potential mu_+ - lambda_+ due to
trap at position (x,y,z) """
if trap is None:
trap = Hulet()
r = np.hypot(x, y)
return -trap.U(r, z)
[docs]def nab(x, y, z, mua, mub, trap=None, model=None):
"""Return the density at point (x,y,z) in specified trap with
chemical potentials (mua,mub).
>>> import matplotlib.pyplot as plt
>> plt.clf()
>>> mua = 12.0
>>> mub = 1.0
>>> z = np.linspace(0,500*UNITS.micron,100)
>>> (na,nb) = np.vectorize(nab)(0,0,z,mua,mub)
>> p=plot(z,na,z,nb)
>>> r = np.linspace(0,500*UNITS.micron,100)
>>> (na,nb) = np.vectorize(nab)(0,r,0,mua,mub)
>> p=plot(r,na,r,nb)
>> plt.show()
"""
if trap is None:
trap = Hulet()
if model is None:
model = Model()
mu_eff = mup_eff(x, y, z, trap)
(mua, mub) = (mua + mu_eff, mub + mu_eff)
return model.nab(mua, mub)