-
Notifications
You must be signed in to change notification settings - Fork 2
/
calculate_r_peak.py
executable file
·152 lines (131 loc) · 5.56 KB
/
calculate_r_peak.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
#! /usr/bin/env python3
"""
Script for calculating radius of pressure maximum given inner and outer edges
of equilibrium torus.
Usage:
calculate_r_peak.py <fm or c> <spin> <r_in> <r_out>
This covers both the Fishbone-Moncrief (FM: 1976 ApJ 207 962) and Chakrabarti
(C: 1985 ApJ 288 1) tori. The latter is generally thinner for the same inputs.
Formulas also reference Abramowicz, Jaroszynski, & Sikora (AJS: 1978 AA 63 21)
and Penna, Kulkarni, & Narayan (PKN: 2013 AA 559 A116).
"""
# Python standard modules
import argparse
import warnings
# Numerical modules
import numpy as np
from scipy.optimize import brentq
# Main function
def main(**kwargs):
# Numerical parameter
r_in_factor = 1.01
# Calculate peak radius for Fishbone-Moncrief torus
if kwargs['torus_type'] == 'fm':
def res(r):
f_in = fm_f(kwargs['spin'], kwargs['r_in'], np.pi / 2.0, r)
f_out = fm_f(kwargs['spin'], kwargs['r_out'], np.pi / 2.0, r)
return f_in - f_out
r_peak = brentq(res, kwargs['r_in'], kwargs['r_out'])
# Calculate peak radius for Chakrabarti torus
if kwargs['torus_type'] == 'c':
def res(r):
c, n = c_cn(kwargs['spin'], kwargs['r_in'], r)
l_in = c_l(kwargs['spin'], kwargs['r_in'], np.pi / 2.0, c, n)
l_out = c_l(kwargs['spin'], kwargs['r_out'], np.pi / 2.0, c, n)
h_in = c_h(kwargs['spin'], kwargs['r_in'], np.pi / 2.0, l_in, c, n, kwargs['r_in'])
h_out = c_h(kwargs['spin'], kwargs['r_out'], np.pi / 2.0, l_out, c, n, kwargs['r_in'])
return h_out - h_in
r_peak = brentq(res, kwargs['r_in'] * r_in_factor, kwargs['r_out'])
# Report results
print('r_peak: {0:24.16e}'.format(r_peak))
# Geometric factors
def geometry(spin, r, theta):
s = np.sin(theta)
c = np.cos(theta)
delta = r ** 2 - 2.0 * r + spin ** 2
sigma = r ** 2 + spin ** 2 * c ** 2
aa = (r ** 2 + spin ** 2) ** 2 - spin ** 2 * delta * s ** 2
return s, delta, sigma, aa
# Metric factors
def metric(spin, r, theta):
s, _, sigma, _ = geometry(spin, r, theta)
g_tt = -1.0 + 2.0 * r / sigma
g_tphi = -2.0 * spin * r / sigma * s ** 2
g_phiphi = (r ** 2 + spin ** 2 + 2.0 * spin ** 2 * r / sigma * s ** 2) * s ** 2
return g_tt, g_tphi, g_phiphi
# Fishbone-Moncrief l_* (FM 3.8)
def fm_ls(spin, r_peak):
numerator = r_peak ** 4 + spin ** 2 * r_peak ** 2 - 2.0 * spin ** 2 * r_peak - spin * r_peak ** 0.5 * (r_peak ** 2 - spin ** 2)
denominator = r_peak ** 2 - 3.0 * r_peak + 2.0 * spin * r_peak ** 0.5
return r_peak ** -1.5 * numerator / denominator
# Fishbone-Moncrief f (cf. FM 3.6)
def fm_f(spin, r, theta, r_peak):
s, delta, sigma, aa = geometry(spin, r, theta)
ls = fm_ls(spin, r_peak)
with warnings.catch_warnings():
warnings.filterwarnings('ignore', 'invalid value encountered in double_scalars', RuntimeWarning)
warnings.filterwarnings('ignore', 'invalid value encountered in log', RuntimeWarning)
warnings.filterwarnings('ignore', 'divide by zero encountered in double_scalars', RuntimeWarning)
term_1 = 0.5 * np.log(aa / (delta * sigma) + ((aa / (delta * sigma)) ** 2 + 4.0 * ls ** 2 / (delta * s ** 2)) ** 0.5)
term_2 = -0.5 * (1.0 + 4.0 * ls ** 2 * delta * sigma ** 2 / (aa ** 2 * s ** 2)) ** 0.5
term_3 = -2.0 * spin * r * ls / aa
f = term_1 + term_2 + term_3
return f
# Abramowicz-Jaroszynski-Sikora u_t (AJS 5)
def ajs_u_t(spin, r, theta, l):
g_tt, g_tphi, g_phiphi = metric(spin, r, theta)
with warnings.catch_warnings():
warnings.filterwarnings('ignore', 'invalid value encountered in sqrt', RuntimeWarning)
u_t = -((g_tphi ** 2 - g_tt * g_phiphi) / (g_phiphi + 2.0 * l * g_tphi + l ** 2 * g_tt)) ** 0.5
return u_t
# Keplerian l (PKN 4)
def l_k(spin, r):
numerator = 1.0 - 2.0 * spin / r ** 1.5 + spin ** 2 / r ** 2
denominator = 1.0 - 2.0 / r + spin / r ** 1.5
return r ** 0.5 * numerator / denominator
# Von Zeipel parameter (C 3.7, 3.8)
def vz(spin, r, l):
numerator = r ** 3 + spin ** 2 * r + 2.0 * spin * (spin - l)
denominator = r + 2.0 * spin / l - 2.0
return (numerator / denominator) ** 0.5
# Chakrabarti c and n (C 2.14a)
def c_cn(spin, r_in, r_peak):
l_in = l_k(spin, r_in)
l_peak = l_k(spin, r_peak)
lambda_in = vz(spin, r_in, l_in)
lambda_peak = vz(spin, r_peak, l_peak)
n = np.log(l_peak / l_in) / np.log(lambda_peak / lambda_in)
c = l_in / lambda_in ** n
return c, n
# Chakrabarti l (C 2.5)
def c_l(spin, r, theta, c, n):
variable_l_min = 1.0
variable_l_max = 100.0
g_tt, g_tphi, g_phiphi = metric(spin, r, theta)
def res_c_l(l):
return (l / c) ** (2.0 / n) + (l * g_phiphi + l ** 2 * g_tphi) / (g_tphi + l * g_tt)
try:
l_val = brentq(res_c_l, variable_l_min, variable_l_max)
except ValueError:
l_val = np.nan
return l_val
# Chakrabarti h (C 2.16)
def c_h(spin, r, theta, l, c, n, r_in):
l_in = c_l(spin, r_in, np.pi / 2.0, c, n)
u_t = ajs_u_t(spin, r, theta, l)
u_t_in = ajs_u_t(spin, r_in, np.pi / 2.0, l_in)
h = u_t_in / u_t
if n == 1.0:
h *= (l_in / l) ** (c ** 2 / (c ** 2 - 1.0))
else:
h *= abs(1.0 - c ** (2.0 / n) * l ** (2.0 - 2.0 / n)) ** (n / (2.0 - 2.0 * n)) * abs(1.0 - c ** (2.0 / n) * l_in ** (2.0 - 2.0 / n)) ** (n / (2.0 * n - 2.0))
return h
# Parse inputs and execute main function
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('torus_type', choices=('fm','c'), help='"fm" for Fishbone-Moncrief or "c" for Chakrabarti')
parser.add_argument('spin', type=float, help='dimensionless spin')
parser.add_argument('r_in', type=float, help='inner edge in gravitational radii')
parser.add_argument('r_out', type=float, help='outer edge in gravitational radii')
args = parser.parse_args()
main(**vars(args))