M87* Calculations

import math
import matplotlib.pyplot as plt
from mpmath import mp

mp.dps = 50


#https://astronomy.stackexchange.com/questions/8413/what-is-the-temperature-of-an-accretion-disc-surrounding-a-supermassive-black-ho

G = 6.67408e-11 #m³ / kg s²
c = 299792458.0 #m / s
σ = 5.670367e-8 #W / m² K⁴
π = math.pi
k_B = 1.38064852e-23 #J / K
h = 6.62607015e-34 #J s

def r_s(M):
    return 2.0*G*M/(c*c)

def T(r,M,Mdot):
    term1 = ( 3.0 * G * M * Mdot ) / ( 8.0 * π * σ * r*r*r )
    term2 = 1.0 - math.sqrt( r_s(M) / r )
    return (term1*term2)**0.25

def B(λ,T):
    return 2.0*h*c*c / ( λ*λ*λ*λ*λ * ( mp.exp((h*c)/(λ*k_B*T)) - 1.0 ) )



#M87*

M = 7.22e9 * 1.98847e30 #7.22 billion solar masses

Mdot = 1.98847e30 / (10.0*86400.0*365.2524) #1 solar mass every ten years



def plot_T():
    rf0 = 1.0
    rf1 = 4.0

    rfs = []
    dr = (rf1-rf0)*0.0001
    rf = rf0
    while rf <= rf1:
        rfs.append(rf)
        rf += dr

    Ts = [T(rf*r_s(M),M,Mdot) for rf in rfs]

    plt.plot(rfs,Ts)

    plt.xlabel("Schwarzschild Radii")
    plt.ylabel("Temperature (K)")

    plt.title("Accretion Disk Temperature for M87*")

    plt.show()

def plot_bb():
    T = 3.0e9

    λs = [ λfm for λfm in range(1,5000,1) ]
    Bs = [ B(λfm*1.0e-15,T) for λfm in λs ]

    plt.plot(λs,Bs)

    plt.xlabel("λ (fm)")
    plt.ylabel("Spectral Radiance")

    plt.title("Blackbody Radiation")

    plt.show()

plot_T()
#plot_bb()