# Source code for galpy.potential.SpiralArmsPotential

###############################################################################
#  SpiralArmsPotential.py: class that implements the spiral arms potential
#                           from Cox and Gomez (2002)
#
#  https://arxiv.org/abs/astro-ph/0207635
#
#  Phi(r, phi, z) = -4*pi*G*H*rho0*exp(-(r-r0)/Rs)*sum(Cn/(Kn*Dn)*cos(n*gamma)*sech(Kn*z/Bn)^Bn)
###############################################################################
from __future__ import division
from .Potential import Potential
from ..util import conversion
import numpy
[docs]class SpiralArmsPotential(Potential):
"""Class that implements the spiral arms potential from (Cox and Gomez 2002 <https://arxiv.org/abs/astro-ph/0207635>__). Should be used to modulate an existing potential (density is positive in the arms, negative outside; note that because of this, a contour plot of this potential will appear to have twice as many arms, where half are the underdense regions).

.. math::

\\Phi(R, \\phi, z) = -4 \\pi GH \\,\\rho_0 exp \\left( -\\frac{R-r_{ref}}{R_s} \\right) \\sum{\\frac{C_n}{K_n D_n} \\,\cos(n \\gamma) \\,\\mathrm{sech}^{B_n} \\left( \\frac{K_n z}{B_n} \\right)}

where

.. math::
K_n &= \\frac{n N}{R \sin(\\alpha)} \\\\
B_n &= K_n H (1 + 0.4 K_n H) \\\\
D_n &= \\frac{1 + K_n H + 0.3 (K_n H)^2}{1 + 0.3 K_n H} \\\\

and

.. math::
\\gamma = N \\left[\\phi - \\phi_{ref} - \\frac{\\ln(R/r_{ref})}{\\tan(\\alpha)} \\right]

The default of :math:C_n=[1] gives a sinusoidal profile for the potential. An alternative from Cox and Gomez (2002) <https://arxiv.org/abs/astro-ph/0207635>__  creates a density that behaves approximately as a cosine squared in the arms but is separated by a flat interarm region by setting

.. math::
C_n = \\left[\\frac{8}{3 \\pi}\,,\\frac{1}{2} \\,, \\frac{8}{15 \\pi}\\right]

"""
normalize= property() # turn off normalize
[docs]    def __init__(self, amp=1, ro=None, vo=None, amp_units='density',
N=2, alpha=0.2, r_ref=1, phi_ref=0, Rs=0.3, H=0.125, omega=0, Cs=[1]):

"""
NAME:
__init__
PURPOSE:
initialize a spiral arms potential
INPUT:
:amp: amplitude to be applied to the potential (default: 1);
can be a Quantity with units of density. (:math:amp = 4 \\pi G \\rho_0)
:ro: distance scales for translation into internal units (default from configuration file)
:vo: velocity scales for translation into internal units (default from configuration file)
:N: number of spiral arms
:alpha: pitch angle of the logarithmic spiral arms in radians (can be Quantity)
:r_ref: fiducial radius where :math:\\rho = \\rho_0 (:math:r_0 in the paper by Cox and Gomez) (can be Quantity)
:phi_ref: reference angle (:math:\\phi_p(r_0) in the paper by Cox and Gomez) (can be Quantity)
:Rs: radial scale length of the drop-off in density amplitude of the arms (can be Quantity)
:H: scale height of the stellar arm perturbation (can be Quantity)
:Cs: list of constants multiplying the :math:\cos(n \\gamma) terms
:omega: rotational pattern speed of the spiral arms (can be Quantity)
OUTPUT:
(none)
HISTORY:
Started - 2017-05-12  Jack Hong (UBC)

Completed - 2017-07-04 Jack Hong (UBC)
"""

Potential.__init__(self, amp=amp, ro=ro, vo=vo, amp_units=amp_units)
alpha= conversion.parse_angle(alpha)
r_ref= conversion.parse_length(r_ref,ro=self._ro)
phi_ref= conversion.parse_angle(phi_ref)
Rs= conversion.parse_length(Rs,ro=self._ro)
H= conversion.parse_length(H,ro=self._ro)
omega= conversion.parse_frequency(omega,ro=self._ro,vo=self._vo)
self._N = -N  # trick to flip to left handed coordinate system; flips sign for phi and phi_ref, but also alpha.
self._alpha = -alpha  # we don't want sign for alpha to change, so flip alpha. (see eqn. 3 in the paper)
self._sin_alpha = numpy.sin(-alpha)
self._tan_alpha = numpy.tan(-alpha)
self._r_ref = r_ref
self._phi_ref = phi_ref
self._Rs = Rs
self._H = H
self._Cs = self._Cs0 = numpy.array(Cs)
self._ns = self._ns0 = numpy.arange(1, len(Cs) + 1)
self._omega = omega
self._rho0 = 1 / (4 * numpy.pi)
self._HNn = self._HNn0 = self._H * self._N * self._ns0

self.isNonAxi = True   # Potential is not axisymmetric
self.hasC = True       # Potential has C implementation to speed up orbit integrations
self.hasC_dxdv = True  # Potential has C implementation of second derivatives

def _evaluate(self, R, z, phi=0, t=0):
"""
NAME:
_evaluate
PURPOSE:
Evaluate the potential at the given coordinates. (without the amp factor; handled by super class)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: Phi(R, z, phi, t)
HISTORY:
2017-05-12  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

return -self._H * numpy.exp(-(R-self._r_ref) / self._Rs) \
* numpy.sum(self._Cs / Ks / Ds * numpy.cos(self._ns * self._gamma(R, phi - self._omega * t)) / numpy.cosh(Ks * z / Bs) ** Bs,axis=0)

def _Rforce(self, R, z, phi=0, t=0):
"""
NAME:
_Rforce
PURPOSE:
Evaluate the radial force for this potential at the given coordinates. (-dPhi/dR)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the radial force
HISTORY:
2017-05-12  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

He = self._H * numpy.exp(-(R-self._r_ref)/self._Rs)

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

dKs_dR = self._dK_dR(R)
dBs_dR = self._dB_dR(R)
dDs_dR = self._dD_dR(R)

g = self._gamma(R, phi - self._omega * t)
dg_dR = self._dgamma_dR(R)

cos_ng = numpy.cos(self._ns * g)
sin_ng = numpy.sin(self._ns * g)

zKB = z * Ks / Bs
sechzKB = 1 / numpy.cosh(zKB)

return -He * numpy.sum(self._Cs * sechzKB**Bs / Ds * ((self._ns * dg_dR / Ks * sin_ng
+ cos_ng * (z * numpy.tanh(zKB) * (dKs_dR/Ks - dBs_dR/Bs)
- dBs_dR / Ks * numpy.log(sechzKB)
+ dKs_dR / Ks**2
+ dDs_dR / Ds / Ks))
+ cos_ng / Ks / self._Rs),axis=0)

def _zforce(self, R, z, phi=0, t=0):
"""
NAME:
_zforce
PURPOSE:
Evaluate the vertical force for this potential at the given coordinates. (-dPhi/dz)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the vertical force
HISTORY:
2017-05-25  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)
zK_B = z * Ks / Bs

return -self._H * numpy.exp(-(R-self._r_ref) / self._Rs) \
* numpy.sum(self._Cs / Ds * numpy.cos(self._ns * self._gamma(R, phi - self._omega * t))
* numpy.tanh(zK_B) / numpy.cosh(zK_B)**Bs,axis=0)

def _phiforce(self, R, z, phi=0, t=0):
"""
NAME:
_phiforce
PURPOSE:
Evaluate the azimuthal force in cylindrical coordinates. (-dPhi/dphi)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the azimuthal force
HISTORY:
2017-05-25  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

g = self._gamma(R, phi - self._omega * t)
Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

return -self._H * numpy.exp(-(R-self._r_ref) / self._Rs) \
* numpy.sum(self._N * self._ns * self._Cs / Ds / Ks / numpy.cosh(z * Ks / Bs)**Bs * numpy.sin(self._ns * g),axis=0)

def _R2deriv(self, R, z, phi=0, t=0):
"""
NAME:
_R2deriv
PURPOSE:
Evaluate the second (cylindrical) radial derivative of the potential.
(d^2 potential / d R^2)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the second radial derivative
HISTORY:
2017-05-31  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

Rs = self._Rs
He = self._H * numpy.exp(-(R-self._r_ref)/self._Rs)

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

dKs_dR = self._dK_dR(R)
dBs_dR = self._dB_dR(R)
dDs_dR = self._dD_dR(R)

R_sina = R * self._sin_alpha
HNn_R_sina = self._HNn / R_sina
HNn_R_sina_2 = HNn_R_sina**2
x = R * (0.3 * HNn_R_sina + 1) * self._sin_alpha

d2Ks_dR2 = 2 * self._N * self._ns / R**3 / self._sin_alpha
d2Bs_dR2 = HNn_R_sina / R**2 * (2.4 * HNn_R_sina + 2)
d2Ds_dR2 = self._sin_alpha / R / x * (self._HNn* (0.18 * self._HNn * (HNn_R_sina + 0.3 * HNn_R_sina_2 + 1) / x**2
+ 2 / R_sina
- 0.6 * HNn_R_sina * (1 + 0.6 * HNn_R_sina) / x
- 0.6 * (HNn_R_sina + 0.3 * HNn_R_sina_2 + 1) / x
+ 1.8 * self._HNn / R_sina**2))

g = self._gamma(R, phi - self._omega * t)
dg_dR = self._dgamma_dR(R)
d2g_dR2 = self._N / R**2 / self._tan_alpha

sin_ng = numpy.sin(self._ns * g)
cos_ng = numpy.cos(self._ns * g)

zKB = z * Ks / Bs
sechzKB = 1 / numpy.cosh(zKB)
sechzKB_Bs = sechzKB**Bs
log_sechzKB = numpy.log(sechzKB)
tanhzKB = numpy.tanh(zKB)
ztanhzKB = z * tanhzKB

return -He / Rs * (numpy.sum(self._Cs * sechzKB_Bs / Ds
* ((self._ns * dg_dR / Ks * sin_ng
+ cos_ng * (ztanhzKB * (dKs_dR/Ks - dBs_dR/Bs)
- dBs_dR / Ks * log_sechzKB
+ dKs_dR / Ks**2
+ dDs_dR / Ds / Ks))
- (Rs * (1 / Ks * ((ztanhzKB * (dBs_dR / Bs * Ks - dKs_dR)
+ log_sechzKB * dBs_dR)
- dDs_dR / Ds) * (self._ns * dg_dR * sin_ng
+ cos_ng * (ztanhzKB * Ks * (dKs_dR/Ks - dBs_dR/Bs)
- dBs_dR * log_sechzKB
+ dKs_dR / Ks
+ dDs_dR / Ds))
+ (self._ns * (sin_ng * (d2g_dR2 / Ks - dg_dR / Ks**2 * dKs_dR)
+ dg_dR**2 / Ks * cos_ng * self._ns)
+ z * (-sin_ng * self._ns * dg_dR * tanhzKB * (dKs_dR/Ks - dBs_dR/Bs)
+ cos_ng * (z * (dKs_dR/Bs - dBs_dR/Bs**2 * Ks) * (1-tanhzKB**2) * (dKs_dR/Ks - dBs_dR/Bs)
+ tanhzKB * (d2Ks_dR2/Ks-(dKs_dR/Ks)**2 - d2Bs_dR2/Bs + (dBs_dR/Bs)**2)))
+ (cos_ng * (dBs_dR/Ks * ztanhzKB * (dKs_dR/Bs - dBs_dR/Bs**2*Ks)
-(d2Bs_dR2/Ks-dBs_dR*dKs_dR/Ks**2) * log_sechzKB)
+ dBs_dR/Ks * log_sechzKB * sin_ng * self._ns * dg_dR)
+ ((cos_ng * (d2Ks_dR2 / Ks**2 - 2 * dKs_dR**2 / Ks**3)
- dKs_dR / Ks**2 * sin_ng * self._ns * dg_dR)
+ (cos_ng * (d2Ds_dR2 / Ds / Ks
- (dDs_dR/Ds)**2 / Ks
- dDs_dR / Ds / Ks**2 * dKs_dR)
- sin_ng * self._ns * dg_dR * dDs_dR / Ds / Ks))))
- 1 / Ks * (cos_ng / Rs
+ (cos_ng * ((dDs_dR * Ks + Ds * dKs_dR) / (Ds * Ks)
-  (ztanhzKB * (dBs_dR / Bs * Ks - dKs_dR)
+ log_sechzKB * dBs_dR))
+ sin_ng * self._ns * dg_dR)))),axis=0))

def _z2deriv(self, R, z, phi=0, t=0):
"""
NAME:
_z2deriv
PURPOSE:
Evaluate the second (cylindrical) vertical derivative of the potential.
(d^2 potential / d z^2)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the second vertical derivative
HISTORY:
2017-05-26  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

g = self._gamma(R, phi - self._omega * t)
Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)
zKB = z * Ks / Bs
tanh2_zKB = numpy.tanh(zKB)**2

return -self._H * numpy.exp(-(R-self._r_ref)/self._Rs) \
* numpy.sum(self._Cs * Ks / Ds * ((tanh2_zKB - 1) / Bs + tanh2_zKB) * numpy.cos(self._ns * g) / numpy.cosh(zKB)**Bs,axis=0)

def _phi2deriv(self, R, z, phi=0, t=0):
"""
NAME:
_phi2deriv
PURPOSE:
Evaluate the second azimuthal derivative of the potential in cylindrical coordinates.
(d^2 potential / d phi^2)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: d^2 potential / d phi^2
HISTORY:
2017-05-29 Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

g = self._gamma(R, phi - self._omega * t)
Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

return self._H * numpy.exp(-(R-self._r_ref) / self._Rs) \
* numpy.sum(self._Cs * self._N**2. * self._ns**2. / Ds / Ks / numpy.cosh(z*Ks/Bs)**Bs * numpy.cos(self._ns*g),axis=0)

def _Rzderiv(self, R, z, phi=0., t=0.):
"""
NAME:
_Rzderiv
PURPOSE:
Evaluate the mixed (cylindrical) radial and vertical derivative of the potential
(d^2 potential / dR dz).
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: d^2 potential / dR dz
HISTORY:
2017-05-12  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

Rs = self._Rs
He = self._H * numpy.exp(-(R-self._r_ref)/self._Rs)

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

dKs_dR = self._dK_dR(R)
dBs_dR = self._dB_dR(R)
dDs_dR = self._dD_dR(R)

g = self._gamma(R, phi - self._omega * t)
dg_dR = self._dgamma_dR(R)

cos_ng = numpy.cos(self._ns * g)
sin_ng = numpy.sin(self._ns * g)

zKB = z * Ks / Bs
sechzKB = 1 / numpy.cosh(zKB)
sechzKB_Bs = sechzKB**Bs
log_sechzKB = numpy.log(sechzKB)
tanhzKB = numpy.tanh(zKB)

return - He * numpy.sum(sechzKB_Bs * self._Cs / Ds * (Ks * tanhzKB * (self._ns * dg_dR / Ks * sin_ng
+ cos_ng * (z * tanhzKB * (dKs_dR/Ks - dBs_dR/Bs)
- dBs_dR / Ks * log_sechzKB
+ dKs_dR / Ks**2
+ dDs_dR / Ds / Ks))
- cos_ng * ((zKB * (dKs_dR/Ks - dBs_dR/Bs) * (1 - tanhzKB**2)
+ tanhzKB * (dKs_dR/Ks - dBs_dR/Bs)
+ dBs_dR / Bs * tanhzKB)
- tanhzKB / Rs)),axis=0)

def _Rphideriv(self, R, z, phi=0,t=0):
"""
NAME:
_Rphideriv
PURPOSE:
Return the mixed radial and azimuthal derivative of the potential in cylindrical coordinates
(d^2 potential / dR dphi)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the mixed radial and azimuthal derivative
HISTORY:
2017-06-09  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

He = self._H * numpy.exp(-(R - self._r_ref) / self._Rs)

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

dKs_dR = self._dK_dR(R)
dBs_dR = self._dB_dR(R)
dDs_dR = self._dD_dR(R)

g = self._gamma(R, phi - self._omega * t)
dg_dR = self._dgamma_dR(R)

cos_ng = numpy.cos(self._ns * g)
sin_ng = numpy.sin(self._ns * g)
zKB = z * Ks / Bs
sechzKB = 1 / numpy.cosh(zKB)
sechzKB_Bs = sechzKB ** Bs

return - He * numpy.sum(self._Cs * sechzKB_Bs / Ds * self._ns * self._N
* (- self._ns * dg_dR / Ks * cos_ng
+ sin_ng * (z * numpy.tanh(zKB) * (dKs_dR / Ks - dBs_dR / Bs)
+ 1/Ks * (-dBs_dR * numpy.log(sechzKB)
+ dKs_dR / Ks
+ dDs_dR / Ds
+ 1 / self._Rs))),axis=0)

def _phizderiv(self, R, z, phi=0, t=0):
"""
NAME:
_phizderiv
PURPOSE:
Evaluate the mixed azimuthal, vertical derivative for this potential at the given coordinates. (-dPhi/dz)
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: mixed azimuthal, vertical derivative
HISTORY:
2021-04-30 - Jo Bovy (UofT)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)
zK_B = z * Ks / Bs

return -self._H * numpy.exp(-(R-self._r_ref) / self._Rs) \
* numpy.sum(self._Cs / Ds * self._ns * self._N * numpy.sin(self._ns * self._gamma(R, phi - self._omega * t))
* numpy.tanh(zK_B) / numpy.cosh(zK_B)**Bs,axis=0)

def _dens(self, R, z, phi=0, t=0):
"""
NAME:
_dens
PURPOSE:
Evaluate the density. If not given, the density is computed using the Poisson equation
from the first and second derivatives of the potential (if all are implemented).
INPUT:
:param R: galactocentric cylindrical radius
:param z: vertical height
:param phi: azimuth
:param t: time
OUTPUT:
:return: the density
HISTORY:
2017-05-12  Jack Hong (UBC)
"""
if isinstance(R,numpy.ndarray) or isinstance(z,numpy.ndarray):
nR= len(R) if isinstance(R,numpy.ndarray) else len(z)
self._Cs=numpy.transpose(numpy.array([self._Cs0,]*nR))
self._ns=numpy.transpose(numpy.array([self._ns0,]*nR))
self._HNn=numpy.transpose(numpy.array([self._HNn0,]*nR))
else:
self._Cs=self._Cs0
self._ns=self._ns0
self._HNn=self._HNn0

g = self._gamma(R, phi - self._omega * t)

Ks = self._K(R)
Bs = self._B(R)
Ds = self._D(R)

ng = self._ns * g
zKB = z * Ks / Bs
sech_zKB = 1 / numpy.cosh(zKB)
tanh_zKB = numpy.tanh(zKB)
log_sech_zKB = numpy.log(sech_zKB)

# numpy of E as defined in the appendix of the paper.
E = 1 + Ks * self._H / Ds * (1 - 0.3 / (1 + 0.3 * Ks * self._H) ** 2) - R / self._Rs \
- (Ks * self._H) * (1 + 0.8 * Ks * self._H) * log_sech_zKB \
- 0.4 * (Ks * self._H) ** 2 * zKB * tanh_zKB

# numpy array of rE' as define in the appendix of the paper.
rE = -Ks * self._H / Ds * (1 - 0.3 * (1 - 0.3 * Ks * self._H) / (1 + 0.3 * Ks * self._H) ** 3) \
+ (Ks * self._H / Ds * (1 - 0.3 / (1 + 0.3 * Ks * self._H) ** 2)) - R / self._Rs \
+ Ks * self._H * (1 + 1.6 * Ks * self._H) * log_sech_zKB \
- (0.4 * (Ks * self._H) ** 2 * zKB * sech_zKB) ** 2 / Bs \
+ 1.2 * (Ks * self._H) ** 2 * zKB * tanh_zKB

return numpy.sum(self._Cs * self._rho0 * (self._H / (Ds * R)) * numpy.exp(-(R - self._r_ref) / self._Rs)
* sech_zKB**Bs * (numpy.cos(ng) * (Ks * R * (Bs + 1) / Bs * sech_zKB**2
- 1 / Ks / R * (E**2 + rE))
- 2 * numpy.sin(ng)* E * numpy.cos(self._alpha)),axis=0)

def OmegaP(self):
"""
NAME:
OmegaP
PURPOSE:
Return the pattern speed. (used to compute the Jacobi integral for orbits).
INPUT:
:param self
OUTPUT:
:return: the pattern speed
HISTORY:
2017-06-09  Jack Hong (UBC)
"""
return self._omega

def _gamma(self, R, phi):
"""Return gamma. (eqn 3 in the paper)"""
return self._N * (phi - self._phi_ref - numpy.log(R / self._r_ref) / self._tan_alpha)

def _dgamma_dR(self, R):
"""Return the first derivative of gamma wrt R."""
return -self._N / R / self._tan_alpha

def _K(self, R):
"""Return numpy array from K1 up to and including Kn. (eqn. 5)"""
return self._ns * self._N / R / self._sin_alpha

def _dK_dR(self, R):
"""Return numpy array of dK/dR from K1 up to and including Kn."""
return -self._ns * self._N / R**2 / self._sin_alpha

def _B(self, R):
"""Return numpy array from B1 up to and including Bn. (eqn. 6)"""
HNn_R = self._HNn / R

return HNn_R / self._sin_alpha * (0.4 * HNn_R / self._sin_alpha + 1)

def _dB_dR(self, R):
"""Return numpy array of dB/dR from B1 up to and including Bn."""
return -self._HNn / R**3 / self._sin_alpha**2 * (0.8 * self._HNn + R * self._sin_alpha)

def _D(self, R):
"""Return numpy array from D1 up to and including Dn. (eqn. 7)"""
return (0.3 * self._HNn**2 / self._sin_alpha / R
+ self._HNn + R * self._sin_alpha) / (0.3 * self._HNn + R * self._sin_alpha)

def _dD_dR(self, R):
"""Return numpy array of dD/dR from D1 up to and including Dn."""
HNn_R_sina = self._HNn / R / self._sin_alpha

return HNn_R_sina * (0.3 * (HNn_R_sina + 0.3 * HNn_R_sina**2. + 1) / R / (0.3 * HNn_R_sina + 1)**2
- (1/R * (1 + 0.6 * HNn_R_sina) / (0.3 * HNn_R_sina + 1)))