from math import pi, sqrt
G, c, Ms = 6.67430e-11, 299792458.0, 1.98847e30
R = 12_000.0
def rs(M):
    return 2*G*M/c**2
def nominal_density(M, radius):
    return 3*M/(4*pi*radius**3)
def pressure_ratio(q):
    assert 0 <= q < 8/9
    return (1-sqrt(1-q))/(3*sqrt(1-q)-1)

# Table 1 and uniform-density pressure diagnostic
for solar_mass in (1.4, 2.0, 3.0, 4.0):
    M = solar_mass*Ms
    q = rs(M)/R
    ratio = pressure_ratio(q) if q < 8/9 else None
    print(solar_mass, rs(M)/1000, q, ratio)
print("Horizon equality mass", R*c**2/(2*G*Ms))
print("Buchdahl limiting mass", 4*R*c**2/(9*G*Ms))

# Tables 2 and 3
for solar_mass in (3.0, 10.0, 4e6, 1e9):
    M = solar_mass*Ms
    tau = 4*G*M/(3*c**3)
    print(solar_mass, rs(M)/1000,
          nominal_density(M, rs(M)), tau)

# Reference star and capture estimate
mstar, MBH = 1.4*Ms, 1e6*Ms
rho = nominal_density(mstar, R)
print("Reference density", rho)
print("Central pressure Pa", rho*c**2*pressure_ratio(rs(mstar)/R))
print("Dynamical time seconds", sqrt(R**3/(G*mstar)))
rt = R*(MBH/mstar)**(1/3)
print("Tidal radius km and ratio", rt/1000, rt/rs(MBH))