How can I numerically (and efficiently) integrate and plot my function?

Question

I want to do a plot of this equation below:

Problem 1: You see... since my function is a function of ν I have to calculate my integral to each ν in my domain. My question is: what is the best way to do that?

I thought about using scipy to do the integral and a for loop to calculate it several times to each ν, but it seems a very inelegant way to solve my problem. Does someone know a better alternative? Does someone have a different idea?

Problem 2: When I write my code I get some errors, mainly because I think that the exponential has a very small expoent. Do you have any ideas of how should I change it so I can do this plot using Python?

Oh, if you try with a different method, it is supposed to look like this

Here is the code I was working on. I'm coming back to Python now, so maybe there is some errors. The plot I'm getting is very different from the one that this is supposed to look.

from scipy.integrate import quad
from scipy.constants import c, Planck, k, pi
import numpy as np
import matplotlib.pyplot as plt

def luminosity_integral(r, x):
    T_est = 4000 
    R_est = 2.5 * (696.34*1e6)

    Temp = ((2/(3*pi))**(1/4)) * T_est * ((R_est/r)**(3/4))
    termo1 = ((4 * (pi**2) * Planck * (x**4) ) / (c**2))
    termo2 = ((Planck * x) / (k*Temp))
   

    return ((termo1 * r ) / (np.exp(termo2) - 1))


freqs = np.linspace(1e10, 1e16)
y = np.array([])
for i in freqs:
I = quad(luminosity_integral, (6 * 2.5 * (696.34*1e6)), (7e4 *  2.5 * (696.34*1e6)), args = (i))
temp = np.array([I[0]])
y = np.concatenate((y, temp))


plt.loglog(freqs, y)
plt.show()

user_12750353 · Answer

Reuse the term R_est instead instead of writing its expression 3 times (better if you want to change that parameter).
you used a pi**2 in the constant multiplying the integral (don't affect the shape)
The shape resembles what you put as reference, but not in the suggested range.
You are using the value of T as T_*, are you sure about that?

Try this version of the code

from scipy.integrate import quad
from scipy.constants import c, Planck, k, pi
import numpy as np
import matplotlib.pyplot as plt

R_est = 2.5 * (696.34e6)
def luminosity_integral(r, x):
    T_est = 4000 

    termo1 = ((4 * pi * Planck * (x**4) ) / (c**2))
    termo2 = ((Planck * x) / (k*T_est)) * (3*pi/2 * (r/R_est)**3)**0.25
    termo3 = np.exp(-termo2)
    return ((termo1 * r ) * termo3 / (1 - termo3))


freqs = np.logspace(6, 16)
y = np.zeros_like(freqs)

for i, nu in enumerate(freqs):
    y[i] = quad(luminosity_integral, (6* R_est), (7e4 *  R_est), args = (nu))[0]
    

plt.loglog(freqs, y)
plt.ylim([1e6, 1e25])
plt.show()

How can I numerically (and efficiently) integrate and plot my function?

1 Answers