-
Notifications
You must be signed in to change notification settings - Fork 95
/
FlattenedPowerPotential.py
154 lines (136 loc) · 5.48 KB
/
FlattenedPowerPotential.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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
###############################################################################
# FlattenedPowerPotential.py: Power-law potential that is flattened in the
# potential (NOT the density)
#
# amp
# phi(R,z)= --------- ; m^2 = R^2 + z^2/q^2
# m^\alpha
###############################################################################
import numpy
from ..util import conversion
from .Potential import Potential
_CORE = 10**-8
class FlattenedPowerPotential(Potential):
"""Class that implements a power-law potential that is flattened in the potential (NOT the density)
.. math::
\\Phi(R,z) = -\\frac{\\mathrm{amp}\\,r_1^\\alpha}{\\alpha\\,\\left(R^2+(z/q)^2+\\mathrm{core}^2\\right)^{\\alpha/2}}
and the same as LogarithmicHaloPotential for :math:`\\alpha=0`
See Figure 1 in `Evans (1994) <http://adsabs.harvard.edu/abs/1994MNRAS.267..333E>`_ for combinations of alpha and q that correspond to positive densities
"""
def __init__(
self,
amp=1.0,
alpha=0.5,
q=0.9,
core=_CORE,
normalize=False,
r1=1.0,
ro=None,
vo=None,
):
"""
Initialize a flattened power-law potential.
Parameters
----------
amp : float or Quantity, optional
Amplitude to be applied to the potential. Can be a Quantity with units of velocity squared.
alpha : float, optional
Power-law exponent.
q : float, optional
Flattening parameter.
core : float or Quantity, optional
Core radius.
normalize : bool or float, optional
If True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.
r1 : float or Quantity, optional
Reference radius for amplitude.
ro : float or Quantity, optional
Distance scale for translation into internal units (default from configuration file).
vo : float or Quantity, optional
Velocity scale for translation into internal units (default from configuration file).
Notes
-----
- 2013-01-09 - Written - Bovy (IAS)
"""
Potential.__init__(self, amp=amp, ro=ro, vo=vo, amp_units="velocity2")
core = conversion.parse_length(core, ro=self._ro)
r1 = conversion.parse_length(r1, ro=self._ro)
self.alpha = alpha
self.q2 = q**2.0
self.core2 = core**2.0
# Back to old definition
self._amp *= r1**self.alpha
if normalize or (
isinstance(normalize, (int, float)) and not isinstance(normalize, bool)
): # pragma: no cover
self.normalize(normalize)
self.hasC = True
self.hasC_dxdv = True
self.hasC_dens = True
def _evaluate(self, R, z, phi=0.0, t=0.0):
if self.alpha == 0.0:
return 1.0 / 2.0 * numpy.log(R**2.0 + z**2.0 / self.q2 + self.core2)
else:
m2 = self.core2 + R**2.0 + z**2.0 / self.q2
return -(m2 ** (-self.alpha / 2.0)) / self.alpha
def _Rforce(self, R, z, phi=0.0, t=0.0):
if self.alpha == 0.0:
return -R / (R**2.0 + z**2.0 / self.q2 + self.core2)
else:
m2 = self.core2 + R**2.0 + z**2.0 / self.q2
return -(m2 ** (-self.alpha / 2.0 - 1.0)) * R
def _zforce(self, R, z, phi=0.0, t=0.0):
if self.alpha == 0.0:
return -z / self.q2 / (R**2.0 + z**2.0 / self.q2 + self.core2)
else:
m2 = self.core2 + R**2.0 + z**2.0 / self.q2
return -(m2 ** (-self.alpha / 2.0 - 1.0)) * z / self.q2
def _R2deriv(self, R, z, phi=0.0, t=0.0):
if self.alpha == 0.0:
denom = 1.0 / (R**2.0 + z**2.0 / self.q2 + self.core2)
return denom - 2.0 * R**2.0 * denom**2.0
else:
m2 = self.core2 + R**2.0 + z**2.0 / self.q2
return -(m2 ** (-self.alpha / 2.0 - 1.0)) * (
(self.alpha + 2) * R**2.0 / m2 - 1.0
)
def _z2deriv(self, R, z, phi=0.0, t=0.0):
if self.alpha == 0.0:
denom = 1.0 / (R**2.0 + z**2.0 / self.q2 + self.core2)
return denom / self.q2 - 2.0 * z**2.0 * denom**2.0 / self.q2**2.0
else:
m2 = self.core2 + R**2.0 + z**2.0 / self.q2
return (
-1.0
/ self.q2
* m2 ** (-self.alpha / 2.0 - 1.0)
* ((self.alpha + 2) * z**2.0 / m2 / self.q2 - 1.0)
)
def _dens(self, R, z, phi=0.0, t=0.0):
if self.alpha == 0.0:
return (
1.0
/ 4.0
/ numpy.pi
/ self.q2
* (
(2.0 * self.q2 + 1.0) * self.core2
+ R**2.0
+ (2.0 - 1.0 / self.q2) * z**2.0
)
/ (R**2.0 + z**2.0 / self.q2 + self.core2) ** 2.0
)
else:
m2 = self.core2 + R**2.0 + z**2.0 / self.q2
return (
1.0
/ self.q2
* (
self.core2 * (1.0 + 2.0 * self.q2)
+ R**2.0 * (1.0 - self.alpha * self.q2)
+ z**2.0 * (2.0 - (1.0 + self.alpha) / self.q2)
)
* m2 ** (-self.alpha / 2.0 - 2.0)
/ 4.0
/ numpy.pi
)