# Source code for contactpp.swell

"""
Square well pseudopotential.
"""

import numpy as np
from scipy.optimize import brentq

class SquareWellPotential(object):
[docs]    """
Object representing a Square well potential.

Attributes
----------

height : float
Potential height (negative for a well).

"""

self.height = height
self.height,0.)

@property
def V(self):
return self._V

class _SquareWellBaseGenerator(object):
"""
Base class for square well generators.
"""

def __init__(self,a,c):
self.a = a
self.c = c # value or None
self._check_inputs()

def make_pseudopotential(self):
if self.c is None:
R = self._calc_V_R()
else:
R = self.c
V = self._calc_V_from_R(R)
return SquareWellPotential(V,R)

class TopHatGenerator(_SquareWellBaseGenerator):

def _calc_gamma(self,R,V):
return R*np.sqrt(V)

def _check_inputs(self):
if self.c is not None:
if self.a >= self.c:
raise ValueError("Cannot generate a top hat potential for \
scattering lengths a > c, where c is the top hat radius.")
if self.a < 0.:
raise ValueError("Scattering length a <= 0.")

def _calc_a(self,R,V):
g = self._calc_gamma(R,V)
if g<1e-8:
# Use Taylor expansion
return R*((g**2)/3. - (g**4)*2./15. + (g**6)*17./315.)
else:
return R*(1.-np.tanh(g)/g)

def _calc_reff(self,R,V):
g = self._calc_gamma(R,V)
return R* ( 1 + (3*np.tanh(g)-g*(3+g**2) )/(3*g*(g - np.tanh(g))**2))

def _calc_V_from_R(self,R):
f = lambda V: self._calc_a(R,V)-self.a
xlow, xhigh = bracket_root(f, 0.)
return brentq(f,xlow,xhigh)

def _calc_V_R(self):
"""
Calculate V and R such that the potential has zero effective range.
"""
def objective_func(R):
V = self._calc_V_from_R(R)
reff = self._calc_reff(R,V)
return reff
Rlow, Rhigh = bracket_root(
objective_func,1.001*self.a,init_step=1e-3)
R = brentq(objective_func,Rlow,Rhigh)
return R

class SquareWellGenerator(_SquareWellBaseGenerator):

def _calc_gamma(self,R,V):
return R*np.sqrt(abs(V))

def _check_inputs(self):
if self.a > 0.:
raise ValueError("Scattering length a >= 0.")

def _calc_a(self,R,V):
g = self._calc_gamma(R,V)
if g < 1e-8:
# Use Taylor expansion
return -R*((g**2)/3. + (g**4)*2./15. +(g**6)*17./315.)
else:
return R*(1.-np.tan(g)/g)

def _calc_V_from_R(self,R):
f = lambda V: self._calc_a(R,V) - self.a
# bracket a root between V = 0. and V corresponds to first divergence of a.
xlow, xhigh = bracket_root(f,0.,direction=-1.,xmax=-np.pi**2/(4.*R**2)+1e-8)
return brentq(f,xlow,xhigh)

def bracket_root(f,xstart,direction=1,xmax=None,init_step=0.1,step_ratio=1.3):
"""
Bracket the first root of f beyond a certain point.

Returns a tuple (xlow, xhigh) of points around a root of f. The
algorithm will only search x >= xstart if direction == 1, and only
x <= xstart if direction == -1. Stop search if no root is found
between xstart and xmax.

Arguments
---------

f : scalar function f(x) -> y.
xstart : float
Will start looking for roots at xstart.
direction : +1 or -1, optional.
Direction in which to look for the root. If direction == +1,
search for x >= xstart, else search for x <= xstart.
default: +1
xmax : float, optional
Only search for roots between xstart and xmax. By default,
direction*infty
init_step : float, optional
Initial step size when searching. Reduce this if you think the
root is very close to xstart or if you think the algorithm
is skipping over roots.
step_ratio : float, optional
How much the step size increases at each iteration. Increasing this
increases the likelihood the algorithm will miss roots, at the cost
of computational expense. It's unlikely you would ever want to set
this to less than 1.
Default : 1.3
"""
if xmax is None:
xmax = direction*np.infty
else:
# Check that the sign of xmax is correct
if np.sign(xmax-xstart) != direction:
raise ValueError("Argument xstart has the wrong sign.")
step = init_step
xlow = xstart
xhigh = xlow + direction*step
flow = f(xlow)
fhigh = f(xhigh)
while flow*fhigh >= 0.:
if fhigh == 0.:
# Fallen exactly on a root.
xhigh = min(xmax,xhigh+init_step)
break
elif flow == 0.:
# Starts exactly on a root.
raise ValueError("xstart is a root.")
step *= step_ratio
xlow, xhigh = xhigh, xhigh + direction*step
if xhigh >= xmax:
xhigh = xmax
flow = f(xlow)
fhigh = f(xhigh)
if xhigh == xmax and flow*fhigh >= 0.:
raise RuntimeError(
"Failed to find a root between {} and {}.".format(
xlow,xhigh))
return xlow,xhigh

"""
Create a SquareWellPotential object. Can be either a top-hat
potential (for positive scattering lengths) or a well
(for negative scattering lengths.

Parameters
----------

scattering_length : float
Float indicating the scattering length.

Returns
-------

SquareWellPotential
"""