# -*- coding: utf-8 -*-
r"""Code for extracting properties of two-body systems.
Phase Shifts
============
We consider here short-range potentials such that for $V(r>R) = 0$.
In this case, the potential may be described in terms of the phase
shift of the scattering states.
"""
from copy import copy
from functools import wraps
import numpy
import numpy as np
import scipy as sp
from numpy import pi, exp, sqrt, sinh, cosh, tanh, tan
from scipy.integrate import ode
from scipy.special import gamma
from mmf.math.integrate.odepack import odeint
import scipy.optimize
import mmf.math.interp
from mmf.math.differentiate import differentiate
# Don't set too large or you get a crash! Try
# SoonYongPotential.are_(-2.9) a couple of times with these set to 1e-14.
_rel_tol = 1e-12
_abs_tol = 1e-12
[docs]class Potential(object):
r"""Instances are callable and return compact potentials V(r).
The class allows one to specify parameters. The important
property is that when parameters are modified, the class ensures
that the radius returned by R = get_R() is such that V(R) <
abs_tol.
The parameter M should be the reduced mass of the system.
"""
[docs] def __init__(self,abs_tol=_abs_tol,rel_tol=_rel_tol,Rmin=None):
r"""Construct the potential object."""
self._modified = True
self.rel_tol = rel_tol
self.abs_tol = abs_tol
self.Rmin = Rmin
self.tabulate_options = mmf.math.interp.Options()
self.tabulate_options.plot_spline = False
self.tabulate_options.plot = True
self.verbose=False
def __call__(self,r):
r"""Return V(r)."""
raise Exception("Potential not defined!")
def __repr__(self):
return self.__class__.__name__+":"+str(self.get_ar())
[docs] def dVdV(self,r):
r"""Return V'(r)/V(r). Should be overloaded for speed and/or
accuracy.
"""
h0 = self.R/100.0
dVr = differentiate(self,r,h0=h0)
Vr = self(r)
return dVr/Vr
[docs] def compute_R(self):
r"""Return `R >= self.Rmin` such that `abs(self(R)) <
self.abs_tol`.
"""
if self.Rmin is not None:
R = self.Rmin
else:
R = 0.1
while self.abs_tol <= abs(self(R)):
R *= 2
return R
def _set_R(self, R):
self._R = R
def _get_R(self):
r"""Return a radius for which self(R) <= abs_tol."""
if (not '_R' in self.__dict__ or
self._R < self.Rmin or
abs(self(self._R)) > self.abs_tol):
self._R = self.compute_R()
return self._R
R = property(_get_R, _set_R)
[docs] def get_deq(self,f,jac,U0,method='odeint'):
r"""Return a "deq" object capable of integrating the system
defined by `U' = f(r,U)` and `diff(U',U) = jac(r,U)` with initial
condition `U(0.0) = U0`.
"""
if 'ode' == method:
deq = ode(f,jac)
deq.set_integrator('vode',
rtol=self.rel_tol,
atol=self.abs_tol,
nsteps=1000)
deq.set_initial_value(U0)
elif 'odeint' == method:
atol = self.abs_tol
rtol = self.rel_tol
class DEQ:
def __init__(self,f=f,jac=jac):
def f0(y,t):
return f(t,y)
def df0(y,t):
return jac(t,y)
self.f = f0
self.Dfun = df0
self.y0 = U0
def integrate(self,R,atol=atol,rtol=rtol,
verbose=self.verbose):
infodict = {'istate': -1}
y0 = self.y0
r0 = 0.0
while infodict['istate'] == -1:
y,infodict = odeint(self.f,
y0,
[r0,R],
(),
self.Dfun,
rtol=rtol,
atol=atol,
mxstep=500,
full_output=True,
printmessg=verbose)
y0 = y[1,:]
r0 = infodict['tcur'][0]
if verbose and r0 < R:
print "Integrating interval " + \
`[r0,R]`
return y[1,:]
def set_initial_value(self,y0):
self.y0 = copy(y0)
deq = DEQ()
return deq
[docs] def get_deq_ar(self,method='odeint'):
r"""Return the differential equation (u,du,uu,duu) for finding
the s-wave scattering length and the s-wave effective range.
Uses the following system:
u'(r) = du(r)
du'(r) = 2*m/hbar**2*V(r)*u(r)
uu'(r) = duu(r)
duu'(r) = 2*m/hbar**2*(V(r)*uu(r)-u(r))
with initial conditions
u(0) = 0.0
du(0) = 1.0
uu(0) = 0.0
duu(0) = 0.0
Call (u,du,uu,duu) = deq.integrate(R) to integrate to R. Note
that deq maintains state, so subsequent calls must have
increasing R.
"""
const = 2.0*self.M/self.hbar/self.hbar
def f(r,U,V=self,const=const):
r"""Return vector U' where U = [u,du,uu,duu]:
u' = du
du' = const*V*u
uu' = duu
duu' = const*(V*uu - u)
"""
u,du,uu,duu = U
Vr = V(r)
dU = numpy.array([du,const*Vr*u,duu,const*Vr*uu-u])
return dU
def jac(r,U,V=self,const=const):
r"""Return jacobian for problem."""
u,du,uu,duu = U
Vr = V(r)
return numpy.array([[0.0,1.0,0.0,0.0],
[const*Vr,0.0,0.0,0.0],
[0.0,0.0,0.0,1.0],
[-1.0,0.0,const*Vr,0.0]])
U0 = numpy.array([0.0,1.0,0.0,0.0])
return self.get_deq(f,jac,U0,method=method)
[docs] def get_deq_E(self,E,method='odeint'):
r"""Return solution to energy dependent differential equation
for (u,du) for finding the s-wave phase shift.
Uses the following system:
du'(r) = 2*m/hbar**2*(V(r)-E)*u(r)
u'(r) = du(r)
with initial conditions
u(0) = 0.0
du(0) = 1.0
Call (u,du) = deq.integrate(R) to integrate to R. Note
that deq maintains state, so subsequent calls must have
increasing R.
This is to be used by calling the method integrate(R).
"""
const = 2.0*self.M/self.hbar/self.hbar
def f(r,U,E=E,V=self,const=const):
r"""Return vector U' where U = [u,u']:
du = u'
du' = const*(V(r) - E)*u
"""
dU = numpy.array([U[1],const*(V(r) - E)*U[0]])
return dU
def jac(r,U,E=E,V=self,const=const):
r"""Return jacobian for problem."""
return numpy.array([[0.0,1.0],[const*(V(r) - E),0.0]])
U0 = numpy.array([0.0,1.0])
return self.get_deq(f,jac,U0,method=method)
[docs] def get_deq_gE(self, E, method='odeint'):
r"""Return solution to energy dependent differential equation
for `(Ap, Am)` for finding the s-wave phase shift which is
.. math::
g(E) = \kappa \frac{A_+ - A_-}{A_+ + A_-}
where $E = -\kappa^2$ and the s-wave radial wavefunction has the
asymptotic form
.. math::
u(r>R) = A_+ e^{\kappa r} + A_-e^{-\kappa r}.
We solve the radial Schrödinger equation
.. math::
u''(r) = [\tilde{V}(r) - \tilde{E}]u(r)
where $\tilde{V} = 2mV(r)/\hbar^2$ and $\tilde{E} = 2mE(r)/\hbar^2$
with initial conditions $u(0)=0$, and $u'(0)=1$ using the following
transformation
.. math::
A_{\pm}(r) = \frac{1}{2}\left(u(r) \pm \frac{u'(r)}{\kappa}\right)
e^{\mp\kappa r}
These should be constant at the asymptotic values once $r>R$ is larger
than the support of the potential. These functions satisfy the
following set of coupled first-order equations:
.. math::
\diff{}{r}
\begin{pmatrix}
A_+(r)\\
A_-(r)
\end{pmatrix} = \frac{\tilde{V}(r)}{2\kappa}\begin{pmatrix}
1 & e^{-2\kappa r}\\
-e^{2\kappa r} & -1
\end{pmatrix}
\begin{pmatrix}
A_+(r)\\
A_-(r)
\end{pmatrix}
with initial conditions $A_{\pm}(0) = \pm\frac{1}{2\kappa}$.
"""
const = 2.0*self.M/self.hbar/self.hbar
kappa2 = -const*E
kappa = sqrt(kappa2)
def f(r, U, E=E, V=self,
kappa2=kappa2, kappa=kappa, const=const):
r"""Return vector `U'` where `U = [Ap, Am]`:
Ap' = ((2*kappa**2 + VV)*fp + VV*fm)/2/kappa
Am' = -((2*kappa**2 + VV)*fm + VV*fp)/2/kappa
fp(0) = 1/2/kappa
fm(0) = -1/2/kappa
"""
Ap, Am = U
Vr = V(r)
VV = const*Vr
k2V = 2*kappa2 + VV
exp2kr = np.exp(2*kappa*r)
dU = VV*numpy.array([Ap + Am/exp2kr,
-Ap*exp2kr - Am])/2/kappa
return dU
def jac(r,U,E=E,V=self,
kappa2=kappa2,kappa=kappa,const=const):
r"""Return jacobian for problem."""
Ap, Am = U
Vr = V(r)
VV = const*Vr
k2V = 2*kappa2 + VV
exp2kr = np.exp(2*kappa*r)
return VV*numpy.array([[1, 1.0/exp2kr],
[-exp2kr, -1.0]])/2/kappa
U0 = numpy.array([0.5/kappa, -0.5/kappa])
return self.get_deq(f,jac,U0,method=method)
[docs] def get_deq_gE_old(self, E, method='odeint',
version=1):
r"""Old version in terms of fp and fm.
`version=0` ::
fp' = dfp'
dfp' = (kappa**2 + VV)*fp - V'/V*(kappa*fp - dfp)
fm' = dfm
dfm' = (kappa**2 + VV)*fm + V'/V*(kappa*fm + dfm)
fp(0) = 1/2/kappa
dfp(0) = 1/2
fm(0) = -1/2/kappa
dfm(0) = 1/2
`version=1` ::
fp' = ((2*kappa**2 + VV)*fp + VV*fm)/2/kappa
fm' = -((2*kappa**2 + VV)*fm + VV*fp)/2/kappa
fp(0) = 1/2/kappa
fm(0) = -1/2/kappa
Call (fp,dfp,fm,dfm) = deq.integrate(R) to integrate to R. Note
that deq maintains state, so subsequent calls must have
increasing R.
This is to be used by calling the method integrate(R).
Examples
--------
>>> V = DoubleGaussian(f_v=0.0)
>>> V._gE_neg(E=-1.0, version=0)
-1.874095805...
>>> V._gE_neg(E=-1.0, version=1)
-1.874095805...
>>> V._gE_neg(E=-1.0, version=2)
-1.874095805...
"""
const = 2.0*self.M/self.hbar/self.hbar
kappa2 = -const*E
kappa = sqrt(kappa2)
if 0 == version:
def f(r,U,E=E,V=self,dVdV=self.dVdV,
kappa2=kappa2,kappa=kappa,const=const):
r"""Return vector U' where U = [fp,dfp,fm,dfm]:
fp' = dfp
dfp' = (kappa**2 + const*V)*fp - V'/V*(kappa*fp - dfp)
fm' = dfm
dfm' = (kappa**2 + const*V)*fm + V'/V*(kappa*fm + dfm)
"""
fp,dfp,fm,dfm = U
Vr = V(r)
dVdVr = dVdV(r)
k2V = kappa2+const*Vr
dU = numpy.array([dfp,
k2V*fp - dVdVr*(kappa*fp - dfp),
dfm,
k2V*fm + dVdVr*(kappa*fm + dfm)])
return dU
def jac(r,U,E=E,V=self,dVdV=self.dVdV,
kappa2=kappa2,kappa=kappa,const=const):
r"""Return jacobian for problem."""
fp,dfp,fm,dfm = U
Vr = V(r)
dVdVr = dVdV(r)
k2V = kappa2+const*Vr
return numpy.array([[0.0,1.0,0.0,0.0],
[k2V-dVdVr*kappa,dVdVr,0.0,0.0],
[0.0,0.0,0.0,1.0],
[0.0,0.0,k2V+dVdVr*kappa,dVdVr]])
U0 = numpy.array([0.5/kappa,0.5,-0.5/kappa,0.5])
else:
def f(r, U, E=E, V=self,
kappa2=kappa2, kappa=kappa,const=const):
r"""Return vector `U'` where `U = [fp, fm]`:
fp' = ((2*kappa**2 + VV)*fp + VV*fm)/2/kappa
fm' = -((2*kappa**2 + VV)*fm + VV*fp)/2/kappa
fp(0) = 1/2/kappa
fm(0) = -1/2/kappa
"""
fp, fm = U
Vr = V(r)
VV = const*Vr
k2V = 2*kappa2 + VV
dU = numpy.array([(k2V*fp + VV*fm)/2.0/kappa,
-(k2V*fm + VV*fp)/2.0/kappa])
return dU
def jac(r,U,E=E,V=self,
kappa2=kappa2,kappa=kappa,const=const):
r"""Return jacobian for problem."""
fp, fm = U
Vr = V(r)
VV = const*Vr
k2V = 2*kappa2 + VV
return numpy.array([[k2V/2.0/kappa, VV/2.0/kappa],
[-VV/2.0/kappa, -k2V/2.0/kappa]])
U0 = numpy.array([0.5/kappa, -0.5/kappa])
return self.get_deq(f,jac,U0,method=method)
[docs] def get_ar(self,R=None,method='odeint'):
r"""Return (a,re) where a is the s-wave scattering length and
re is the s-wave effective range. Integrates the wavefunction
to radius R. Answer is independent of R if V(r>R) = 0.
"""
if R is None:
R = self.R
deq = self.get_deq_ar(method)
u,du,uu,duu = deq.integrate(R)
R2 = R*R
R3 = R*R2
#A = u - du*R
#AA = du*R3/6.0 - R2*u/2.0 - R*duu + uu
a = R - u/du
a2 = a*a
re = 2.0*R\
+2.0/3.0*R3/a2\
-2.0*R2/a\
-2.0*uu/du/a2\
+2.0*(R/a - 1.0)*duu/du/a
return (a,re)
def _gE(self,E,R=None,method='odeint'):
r"""Return k*cot(delta) as a function of energy E for the specified
scattering problem by matching at a specified radius R."""
if R is None:
R = self.R
deq = self.get_deq_E(E)
u,du = deq.integrate(R)
kappa = sqrt(2.0*self.M*abs(E))/self.hbar
if 0.0 == kappa:
return du/(u-du*R)
elif E >= 0:
k = kappa
x = du/u
t = tan(k*R)
return (x + k*t)/(1.0-x*t/k)
else:
return (du/u - kappa*tanh(kappa*R))/(1.0-du/u*tanh(kappa*R)/kappa)
def _gE_neg(self, E, R=None, method='odeint', version=2):
r"""Return a k*cot(delta) as a function of energy E for the specified
scattering problem by matching at a specified radius R. Use this with
`E < 0`.
"""
if R is None:
R = self.R
assert E < 0.0
kappa = sqrt(2.0*self.M*abs(E))/self.hbar
if version < 2:
deq = self.get_deq_gE_old(E, version=version)
if 0 == version:
fp, dfp, fm, dfm = deq.integrate(R)
else:
fp, fm = deq.integrate(R)
Ap = fp/exp(kappa*R)
Am = fm/exp(-kappa*R)
else:
deq = self.get_deq_gE(E)
Ap, Am = deq.integrate(R)
return kappa*(Ap-Am)/(Ap+Am)
[docs] def gE(self,E,R=None,method='odeint'):
r"""Return k*cot(delta) as a function of energy E for the specified
scattering problem by matching at a specified radius R."""
if R is None:
R = self.R
if E < 0.0:
return self._gE_neg(E,R=R,method=method)
else:
return self._gE(E,R=R,method=method)
[docs] def plot(self,Rmax=None):
r"""Plot potential."""
if Rmax is None:
Rmax = self.R
Vr = mmf.math.interp.tabulate(self,[0,Rmax],self.tabulate_options)
from matplotlib import pyplot as plt
plt.plot(Vr.x,numpy.array(Vr.y)[0,:])
[docs] def plot_gE(self,Erange,Rmax=None):
r"""Plot potential."""
if Rmax is None:
Rmax = self.R
@np.vectorize
def f(E):
return self.gE(E)
gE_interp = mmf.math.interp.tabulate(f,Erange,self.tabulate_options)
a,re = self.get_ar()
from matplotlib import pyplot as plt
E = gE_interp.x
gE = np.array(gE_interp.y)[0,:]
gE = np.where(gE > 1000,1000,gE)
k2 = E/self.hbar/self.hbar*2.0*self.M
y = -1./a + re*k2/2.0
y_bound = np.where(k2<0,-sqrt(abs(k2)),0.0)
plt.plot(E,gE,'-+')
plt.plot(E,y,':r',E,y_bound,':y')
[docs] def plot_u(self,E,Rmax=None):
r"""Plot radial wavefunction."""
if Rmax is None:
Rmax = self.R
deq = self.get_deq_E(E)
@np.vectorize
def f(r,E=E,deq=deq,U0=numpy.array([0.0,1.0])):
deq.set_initial_value(U0)
u,du = deq.integrate(r)
return u
u_interp = mmf.math.interp.tabulate(f,[0,Rmax],self.tabulate_options)
import pylab
r = u_interp.x
u = numpy.array(u_interp.y)[0,:]
pylab.plot(r,u,'-+')
[docs]class PotentialV0R0(Potential):
r"""Mixing for two-parameter potentials with a parameter $v_0$ governing the
overall strength and a parameter $r_0$ governing the effective range. The
potential can have additional parameters, but these two will be adjusted to
set the scattering-length and range."""
[docs] def set_ainv_re(self, ainv, re):
r"""Uses a non-linear root-finder to set the parameters `v0` and `r0`
so that the scattering length is `1.0/ainv` and the effective range is
`re`.
>>> V = PoschlTeller()
>>> V.set_ainv_re(ainv=0.0,re=1.0)
>>> np.allclose(V.r0, 1.0)
True
>>> np.allclose(V.v0, -1.0)
True
"""
def f(x,self=self,ainv=ainv,re=re):
r"""Function to zero: [ainv - 1/a(V),re-re(V)].
x = [v0, r0]"""
self.v0 = x[0]
self.r0 = x[1]
(a_, re_) = self.get_ar()
return [ainv - 1.0/a_,re-re_]
(x,infodict,ier,mesg) = \
scipy.optimize.fsolve(func=f,
x0=[self.v0, self.r0],
xtol=1.49012e-08,
maxfev=10000,
full_output=1)
if (ier != 1):
print "Error!!!"
print mesg
print infodict
x = [-1.0,0.0]
self.v0 = x[0]
self.r0 = x[1]
[docs]class PoschlTeller(PotentialV0R0):
r"""Extended Poschl-Teller (Rosen-Morse) potential used by Soon
Yong, Alex, Joe, and Stefano:
.. math::
V(r) &= v_0 \frac{\hbar^2}{M}\frac{\mu_0^2}{\cosh^2(\mu_0 r)}\\
V(r) &= v_0 \frac{\hbar^2}{M}\frac{4}{r_0^2}\sech^2\frac{2r}{r_0}\\
Units are set by $k_F = (3\pi^2N_+)^{1/3}/L$ for a box of length $L$ and
$\mu = A k_F$. As a test, at unitarity, $v_0 = 1$ and $r_e = 2/\mu_0$ for
two-body scattering where $m$ here is the individual mass of the particles
(`M=0.5` by default).
>>> V = PoschlTeller(v0=-1.0)
>>> (a,re) = V.get_ar()
>>> mu = 2.0/V.r0
>>> abs(a) > 6e8
True
>>> abs(re - 2.0/mu) < 1e-9
True
It is resonant $(a=\infty)$ for $v_0 = -n(2n - 1)$ for $n \in [1,2,3,...]$.
>>> def infty(n):
... return abs(PoschlTeller.a_(-n*(2.0*n-1.0)))
>>> infty(1) > 1e8
True
>>> infty(2) > 1e8
True
At these resonances one has $r_e/r_0 = H(2n-1)$ where $H(n)$ is the
Harmonic Number $H(n) = \sum_{j=1}^{n}j^{-1}$.
>>> def H(n):
... return (1./numpy.arange(1,n+1)).sum()
>>> def zero(n):
... return abs(PoschlTeller.re_(-n*(2.0*n-1.0))-H(2*n-1))
>>> zero(1) < 1e-10
True
>>> zero(2) < 1e-10
True
The scattering length is zero at the following points where the
effective range diverges. I am not sure yet what these correspond
to::
v0 = [0.0,-4.4438633257625231,-12.913326147200248,
-25.421577034814653,-41.955714059176017,
-62.508749555456298,...]
# kF_SY = (3.*pi)**(2./3.)
# EFG_SY = kF_SY**2/2.0*3.0/5.0
# EFG = EFG_SY/4.0/pi**2
# re = 2.0*L/(6.2527*(3.0*pi)**(2./3.))
#
# V = -2*hbar^2*mu^2*v0/m/cosh(mu*r)**2
# v0 = 1.0 => ainv = 0, re = 2.0/mu
# v0 = 1.30119 => kfa_SY = 1.0
"""
[docs] def __init__(self,
v0=-1.0, r0=2.0/6.2527,
hbar=1.0, M=0.5,
abs_tol=_abs_tol,
rel_tol=_rel_tol,
Rmin=0.1):
Potential.__init__(self,
abs_tol=abs_tol,
rel_tol=rel_tol,
Rmin=Rmin)
self.v0 = v0
self.r0 = r0
self.hbar = hbar
self.M = M
def __repr__(self):
return self.__class__.__name__+":"+\
str(([self.v0,self.r0],self.get_ar()))
def __call__(self,r):
#kF = (3*pi*pi*N)**(1./3.)/L
#mu0 = self.A0*kF
mu0 = 2.0/self.r0
hbar = self.hbar
M = self.M
m = self.M*2.0
tmp = hbar*hbar/self.M
v0 = self.v0*mu0*mu0*tmp
return v0/(cosh(mu0*r))**2
@property
[docs] def mu0(self):
r"""Stefano, Alex, etc. use this parameter instead, so we provide an
alias."""
return 2.0/self.r0
[docs] def dVdV(self,r):
r"""Return $V'(r)/V(r)$
.. code-block:: maple
V:=v0/cosh(mu0*r)^2+V1/cosh(mu1*r)^2;
with(codegen);
C(simplify(diff(V,r)/V),optimized);
"""
mu0 = 2.0/self.r0
hbar = self.hbar
M = self.M
m = self.M*2.0
tmp = hbar*hbar/M
v0 = self.v0*mu0*mu0*tmp
return -2.0*mu0*tanh(mu0*r)
[docs] def set_ainv(self,ainv=0.0):
r"""Set the coefficient `v0` so that the inverse scattering length is
`ainv`.
>>> V = PoschlTeller()
>>> V.v0 = -0.5
>>> V.set_ainv(0.0)
>>> np.allclose(V.v0, -1.0)
True
"""
def f(v0,self=self,ainv=ainv):
r"""Function to zero: ainv - 1/a(V)."""
self.v0 = v0
(a,re) = self.get_ar()
return ainv - 1.0/a
(v0_,infodict,ier,mesg) = \
scipy.optimize.fsolve(func=f,
x0 = self.v0,
xtol=1.49012e-08,
maxfev=10000,
full_output=1)
if (ier != 1):
print "Error!!!"
print mesg
print infodict
v0_ = -1.0
self.v0 = v0_
[docs] def set_ainv_re(self,ainv=0.0,re=0.07):
r"""Set the coefficients `v0` and `r0` so that the inverse scattering
length is `ainv` and the effective range is `re`.
>>> V = PoschlTeller(r0=1.0)
>>> re = 1.0
>>> V.set_ainv_re(ainv=0.0,re=re)
>>> np.allclose(V.r0, re)
True
>>> np.allclose(V.v0, -1.0)
True
"""
def f(x,self=self,ainv=ainv,re=re):
r"""Function to zero: [ainv - 1/a(V),re-re(V)].
x = [v0,r0]"""
self.v0 = x[0]
self.r0 = x[1]
(a_,re_) = self.get_ar()
return [ainv - 1.0/a_,re-re_]
(x,infodict,ier,mesg) = \
scipy.optimize.fsolve(func=f,
x0=[self.v0,self.r0],
xtol=1.49012e-08,
maxfev=10000,
full_output=1)
if (ier != 1):
print "Error!!!"
print mesg
print infodict
x = [-1.0,0.0]
self.v0 = x[0]
self.r0 = x[1]
@staticmethod
@numpy.vectorize
[docs] def are_(v0,Vtmp=[None]):
r"""Return the dimensionless function `a(v0,r0)/r0` which is
characteristic of this potential.
>>> abs(PoschlTeller.a_(-1.0)) > 1e6
True
>>>
"""
if Vtmp[0] is None:
Vtmp[0] = PoschlTeller(r0=1.0)
Vtmp[0].v0 = v0
return Vtmp[0].get_ar()
@staticmethod
@numpy.vectorize
[docs] def a_(v0):
r"""Return the dimensionless function `a(v0,r0)/r0` which is
characteristic of this potential.
"""
return PoschlTeller.are_(v0)[0]
@staticmethod
@numpy.vectorize
[docs] def re_(v0):
r"""Return the dimensionless function `re(v0,r0)/r0` which is
characteristic of this potential.
"""
return PoschlTeller.are_(v0)[1]
[docs] def k_cot_delta_analytic(self, E, l=None):
r"""This potential has some analytic properties (see Alex's thesis and
references therein):
.. math::
k\cot\delta_k &= -k\frac{\Im{Z}}{\Re{Z}},\\
Z &= \frac{2^{-i\tilde{k}}\Gamma(i\tilde{k})}
{\Gamma\left(\frac{1+\lambda + i\tilde{k}}{2}\right)
\Gamma\left(\frac{2-\lambda + i\tilde{k}}{2}\right)},\\
\tilde{k} &= k/\mu, \qquad
l = \frac{1 + \sqrt{1 - 8v_0}}{2}}
Examples
--------
>>> V = PoschlTeller(v0=-1.0,r0=2.0/6.2527)
>>> V.gE(E=1.0)
0.1599309...
>>> V.k_cot_delta_analytic(E=1.0)
0.1599309...
"""
k2 = abs(E*2*self.M)
if l is None:
l = (1 + np.sqrt(1 - 4*2*self.v0))/2
mu = 2/self.r0
k = sp.sqrt(k2)
k_mu = k/mu
ik_mu = 1j*k_mu*np.sign(E)
z = gamma(ik_mu)*2**(-ik_mu)/gamma((1+l+ik_mu)/2)/gamma((2-l+ik_mu)/2)
cot_delta = -z.imag/z.real
return k*cot_delta
SoonYongPotential = PoschlTeller
[docs]class SoonYongPotential2(Potential):
r"""Extended Potential used by Soon Yong.
V(r) = hbar**2/M*(v0*mu0**2/(cosh(mu0*r))**2 \
+ V1*mu1**2/(cosh(mu1*r))**2)
His units are set by
kF = (3*pi*pi*N)**(1./3.) / L for a box of length L
mu = A*kF
As a test, at unitarity, v0 = 1.0 and re = 2/mu for two-body
scattering where m here is the individual mass of the particles.
>>> V = SoonYongPotential2(v0=-1.0,V1=0.0)
>>> (a,re) = V.get_ar()
>>> mu = V.A0*V.kF
>>> abs(a) > 6e8
True
>>> abs(re - 2.0/mu) < 1e-9
True
# kF_SY = (3.*pi)**(2./3.)
# EFG_SY = kF_SY**2/2.0*3.0/5.0
# EFG = EFG_SY/4.0/pi**2
# re = 2.0*L/(6.2527*(3.0*pi)**(2./3.))
#
# V = -2*hbar^2*mu^2*v0/m/cosh(mu*r)**2
# v0 = 1.0 => ainv = 0, re = 2.0/mu
# v0 = 1.30119 => kfa_SY = 1.0
"""
[docs] def __init__(self,
v0=-1.0,A0=6.2527,
V1=0.0,A1=3.0,
hbar=1./2./pi,M=0.5,kF=1.0,
abs_tol=_abs_tol,
rel_tol=_rel_tol,
Rmin=0.1):
Potential.__init__(self,
abs_tol=abs_tol,
rel_tol=rel_tol,
Rmin=Rmin)
self.v0 = v0
self.V1 = V1
self.A0 = A0
self.A1 = A1
self.hbar = hbar
self.M = M
self.kF = kF
def __call__(self,r):
#kF = (3*pi*pi*N)**(1./3.)/L
#mu0 = self.A0*kF
#mu1 = self.A1*kF
mu0 = self.A0*self.kF
mu1 = self.A1*self.kF
hbar = self.hbar
M = self.M
m = self.M*2.0
tmp = hbar*hbar/self.M
v0 = self.v0*mu0*mu0*tmp
V1 = self.V1*mu1*mu1*tmp
return v0/(cosh(mu0*r))**2\
+ V1/(cosh(mu1*r))**2
[docs] def dVdV(self,r):
r"""Return V'(r)/V(r).
V:=v0/cosh(mu0*r)^2+V1/cosh(mu1*r)^2;
with(codegen);
C(simplify(diff(V,r)/V),optimized);
"""
mu0 = self.A0*self.kF
mu1 = self.A1*self.kF
hbar = self.hbar
M = self.M
m = self.M*2.0
tmp = hbar*hbar*2.0/m
v0 = self.v0*mu0*mu0*tmp
V1 = self.V1*mu1*mu1*tmp
t1 = mu0*r;
t2 = sinh(t1);
t4 = mu1*r;
t5 = cosh(t4);
t6 = t5*t5;
t10 = sinh(t4);
t12 = cosh(t1);
t13 = t12*t12;
t28 = -2.0*(v0*t2*mu0*t6*t5+V1*t10*mu1*t13*t12)/t12/t5/(v0*t6+V1*t13);
return t28
[docs] def set_ainv(self,ainv=0.0):
r"""Set the coefficient v0 so that the scattering length is
a.
>>> V = SoonYongPotential2()
>>> V.v0 = -0.5
>>> V.set_ainv(0.0)
>>> np.allclose(V.v0, -1.0)
True
"""
def f(v0,self=self,ainv=ainv):
r"""Function to zero: ainv - 1/a(V)."""
self.v0 = v0
(a,re) = self.get_ar()
return ainv - 1.0/a
(v0_,infodict,ier,mesg) = \
scipy.optimize.fsolve(func=f,
x0 = self.v0,
xtol=1.49012e-08,
maxfev=10000,
full_output=1)
if (ier != 1):
print "Error!!!"
print mesg
print infodict
v0_ = -1.0
self.v0 = v0_
[docs] def set_ainv_re(self,ainv=0.0,re=0.07):
r"""Set the coefficients v0 and A0 so that the scattering length is
1.0/ainv and the effective range is re.
>>> V = SoonYongPotential2(kF=1.0)
>>> V.V1 = 0.0
>>> V.set_ainv_re(ainv=0.0,re=1.0)
>>> np.allclose(V.A0, 2./V.kF)
True
>>> np.allclose(V.v0, -1.0)
True
"""
def f(x,self=self,ainv=ainv,re=re):
r"""Function to zero: [ainv - 1/a(V),re-re(V)].
x = [v0,A0]"""
self.v0 = x[0]
self.A0 = x[1]
(a_,re_) = self.get_ar()
return [ainv - 1.0/a_,re-re_]
(x,infodict,ier,mesg) = \
scipy.optimize.fsolve(func=f,
x0=[self.v0,self.A0],
xtol=1.49012e-08,
maxfev=10000,
full_output=1)
if (ier != 1):
print "Error!!!"
print mesg
print infodict
x = [-1.0,0.0]
self.v0 = x[0]
self.A0 = x[1]
[docs] def set_ainv_re_VV(self,ainv=0.0,re=0.07):
r"""Set the coefficients v0 and V1 so that the scattering length is
1.0/ainv and the effective range is re.
>>> V = SoonYongPotential2()
>>> V.A0 = 6.2527
>>> V.A1 = 3.0
>>> V.v0 = -1.5
>>> V.V1 = 0.6
>>> V.set_ainv_re_VV(ainv=0.0,re=0.0)
>>> np.allclose(V.v0, -1.4789900537007772)
True
>>> np.allclose(V.V1, 0.61932931745837538)
True
"""
def f(x,self=self,ainv=ainv,re=re):
r"""Function to zero: ainv - 1/a(V). x = [v0,V1]"""
self.v0 = x[0]
self.V1 = x[1]
(a_,re_) = self.get_ar()
return [ainv - 1.0/a_,re-re_]
(x,infodict,ier,mesg) = \
scipy.optimize.fsolve(func=f,
x0=[self.v0,self.V1],
xtol=1.49012e-08,
maxfev=10000,
full_output=1)
if (ier != 1):
print "Error!!!"
print mesg
print infodict
x = [-1.0,0.0]
self.v0 = x[0]
self.V1 = x[1]
[docs]class DoubleGaussian(PotentialV0R0):
r"""Double Gaussian potential:
.. math::
V(r) &= -v_0 \frac{\hbar^2}{M} \mu_0^2 \left(
\exp[-(\mu_0 r/2)^{2}] + f_v \exp[-(\mu_0 r f_r/2)^2]
\right)\\
&= -v_0 \frac{\hbar^2}{M} \mu_0^2 \left(
\exp[-(\mu_0 r/2)^{2}] + f_v \exp[-(\mu_0 r f_r/2)^2]
\right)\\
Note that $\mu_0 = 4/r_0$ here.
>>> v = DoubleGaussian(mu0=12.0)
>>> v.get_ar() # One of Stefano's values to check.
(-2887043864..., 0.3293...)
>>> v.set_ainv_re(ainv=0.0, re=1.0)
>>> v.v0, v.r0
(0.786097607..., 1.012180386...)
"""
[docs] def __init__(self, v0=0.7860976072874742,
r0=4.0/6.0,
mu0=None,
f_v=-4.0, f_r=2.0,
hbar=1.0, M=0.5,
abs_tol=_abs_tol,
rel_tol=_rel_tol,
Rmin=0.1):
Potential.__init__(self,
abs_tol=abs_tol,
rel_tol=rel_tol,
Rmin=Rmin)
self.hbar = hbar
self.M = M
self.f_v = f_v
self.f_r = f_r
self.v0 = v0
if mu0 is None:
self.r0 = r0
else:
self.mu0 = mu0
[docs] def set_mu0(self, mu0):
self.r0 = 4./mu0
[docs] def get_mu0(self):
return 4./self.r0
mu0 = property(get_mu0, set_mu0)
def __call__(self,r):
#kF = (3*pi*pi*N)**(1./3.)/L
#mu0 = self.A0*kF
mu0 = self.mu0
hbar = self.hbar
M = self.M
m = self.M*2.0
tmp = hbar*hbar/self.M
v0 = -self.v0*mu0*mu0*tmp
return v0*(np.exp(-(mu0*r/2)**2)
+ self.f_v*np.exp(-(mu0*r*self.f_r/2)**2))
[docs]class DoubleExponential(PotentialV0R0):
r"""Double Exponential potential:
.. math::
V(r) = -v_0 \frac{\hbar^2}{M} \mu_0^2 \left(
\exp[-\mu_0 r] + f_v \exp[-\mu_0 r f_r]
\right)
Note that $\mu_0 = 4/r_0$ here.
>>> v = DoubleExponential(mu0=14.0)
>>> v.get_ar() # One of Stefano's values to check.
(5852295157..., 0.3099...)
>>> v.set_ainv_re(ainv=0.0, re=1.0)
>>> v.v0, v.r0
(1.190915117966..., 0.9218838916...)
"""
[docs] def __init__(self, v0=1.1909151179666314,
r0=4./6.0,
mu0=None,
f_v=-2.0, f_r=2.0,
hbar=1.0, M=0.5,
abs_tol=_abs_tol,
rel_tol=_rel_tol,
Rmin=0.1):
Potential.__init__(self,
abs_tol=abs_tol,
rel_tol=rel_tol,
Rmin=Rmin)
self.hbar = hbar
self.M = M
self.f_v = f_v
self.f_r = f_r
self.v0 = v0
if mu0 is None:
self.r0 = r0
else:
self.mu0 = mu0
[docs] def set_mu0(self, mu0):
self.r0 = 4./mu0
[docs] def get_mu0(self):
return 4./self.r0
mu0 = property(get_mu0, set_mu0)
def __repr__(self):
return self.__class__.__name__+":"+\
str(([self.v0,self.r0],self.get_ar()))
def __call__(self,r):
#kF = (3*pi*pi*N)**(1./3.)/L
#mu0 = self.A0*kF
mu0 = self.mu0
hbar = self.hbar
M = self.M
m = self.M*2.0
tmp = hbar*hbar/self.M
v0 = -self.v0*mu0*mu0*tmp
return v0*(np.exp(-mu0*r) + self.f_v*np.exp(-mu0*r*self.f_r))
[docs]class BargmannPotential(Potential):
r"""Bargmann Potential where the effective range expansion is exact
with only a and r_e.
>> V = BargmannPotential(a=1.0,re=0.01)
>> (a,re) = V.get_ar()
>> (a,re)
(1.0, 0.1)
"""
[docs] def __init__(self,
a=1.0,
re=0.4,
hbar=1./2./pi,M=0.5,
**kwarg):
Potential.__init__(self,**kwarg)
self.a = a
self.re = re
self.hbar = hbar
self.M = M
def __call__(self,r):
re = self.re
a = self.a
hbar = self.hbar
M = self.M
m = 2.0*M
tmp = sqrt(1.0 - 2.0*re/a)
Ap = tmp + 1.0
Am = tmp - 1.0
b = Ap/re
beta = (Ap-Am)/(Ap+Am)
v0 = -8.0*hbar*hbar/m*b*b
num = beta*exp(-2.0*b*r)
den = 1 + num
den = den*den
return v0*num/den
[docs] def dVdV(self,r):
r"""Return V'(r)/V(r).
V:=v0*beta*exp(-2*b*r)/(1+beta*exp(-2*b*r))^2;
with(codegen);
C(simplify(diff(V,r)/V),optimized);
"""
re = self.re
a = self.a
hbar = self.hbar
M = self.M
m = 2.0*M
tmp = sqrt(1.0 - 2.0*re/a)
Ap = tmp + 1.0
Am = tmp - 1.0
b = Ap/re
beta = (Ap-Am)/(Ap+Am)
t3 = exp(-2.0*b*r);
t4 = beta*t3;
t10 = 2.0*b*(-1.0+t4)/(1.0+t4);
return t10
r"""
# First modify
from numpy import *
from pylab import *
ion()
import two_body as tb;reload(tb);V=tb.SoonYongPotential2()
deq = V.get_ar()
r = linspace(0.01,1.0,100)
U = []
for r0 in r:
U.append(deq.integrate(r0))
U = numpy.array(U)
plot(r,U[:,2])
P0 = polyfit(r[50:],U[50:,0],1)
P2 = polyfit(r[50:],U[50:,2],3)
A = P0[-1]
a = -A/P0[-2]
AA = P2[-1]
re = (P2[-2] + AA/a)*2.0/A
print P2[-3],-A/2
print P2[-4],A/6.0/a
"""