-
Notifications
You must be signed in to change notification settings - Fork 95
/
IsothermalDiskPotential.py
56 lines (39 loc) · 1.78 KB
/
IsothermalDiskPotential.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
###############################################################################
# IsothermalDiskPotential.py: class that implements the one-dimensional
# self-gravitating isothermal disk
###############################################################################
import numpy
from .linearPotential import linearPotential, _APY_LOADED
if _APY_LOADED:
from astropy import units
class IsothermalDiskPotential(linearPotential):
"""Class representing the one-dimensional self-gravitating isothermal disk
.. math::
\\rho(x) = \\mathrm{amp}\\,\\mathrm{sech}^2\\left(\\frac{x}{2H}\\right)
where the scale height :math:`H^2 = \\sigma^2/[8\\pi G \\,\\mathrm{amp}]`. The parameter to setup the disk is the velocity dispersion :math:`\\sigma`.
"""
def __init__(self,amp=1.,sigma=0.1,ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
Initialize an IsothermalDiskPotential
INPUT:
amp - an overall amplitude
sigma - velocity dispersion (can be a Quantity)
OUTPUT:
instance
HISTORY:
2018-04-11 - Written - Bovy (UofT)
"""
linearPotential.__init__(self,amp=amp,ro=ro,vo=vo)
if _APY_LOADED and isinstance(sigma,units.Quantity):
sigma= sigma.to(units.km/units.s).value/self._vo
self._sigma2= sigma**2.
self._H= sigma/numpy.sqrt(8.*numpy.pi*self._amp)
self._amp= 1. # Need to manually set to 1, because amp is now contained in the combination of H and sigma^2
self.hasC= True
def _evaluate(self,x,t=0.):
return 2.*self._sigma2*numpy.log(numpy.cosh(0.5*x/self._H))
def _force(self,x,t=0.):
return -self._sigma2*numpy.tanh(0.5*x/self._H)/self._H