Построение кривой излучения черного тела с помощью matplotlib, но я не понимаю, какую ошибку я получаю?

Question

import math
import random
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

h=6.62607015E-34 # Planck constant, units of J*s
hbar=h/(2*math.pi) # Reduced Planck constant 
k=1.380649E-23 # Boltzmann constant, units of J/K
c=299792458.00 # speed of light M/s
sb=(2*(math.pi**5)*(k**4))/(15*(c**2)*(h**3)) # Stefan-Boltzmann constant, units of W/(m^2 K^4)

def one_star():
    mass=random.uniform(0.1,50) # random mass in relation to the sun's mass
    print("The mass of the star is "+str(mass)+" times the mass of the Sun")
    if mass <= 1.4: # given by the mass-radius relationship
        radius=mass**0.9
        print("The radius of the star is "+str(radius)+" times the radius of the Sun")
    else:
        radius=mass**0.6 # larger stars are less dense
        print("The radius of the star is "+str(radius)+" times the radius of the Sun")

    # To find the volume and surface area of the star, we want values in terms of meters^3 and meters^2
    # So we'll multiply each of the previous values by the Sun's mass and radius.
    mass_actual=mass*1.989E28 # mass in kilograms
    print("The mass of the star in kilograms: "+str(mass_actual))
    radius_actual=radius*696000000 # radius in meters
    print("The radius of the star in meters: "+str(radius_actual))   
    volume=(4/3)*math.pi*(radius_actual**3) # volume in meters cubed
    print("The volume of the star in meters cubed: "+str(volume))    
    sa=4*math.pi*(radius_actual**2) # surface area in meters squared
    print("The surface area of the star in meters squared: "+str(sa))    
    B=random.uniform(3.3,3.7)
    luminosity=mass**B # we are using the value of mass in relation to the sun's mass again for this part
    print("The luminosity of the star is "+str(luminosity)+" times the luminosity of the Sun")
    # Time to find the temperature
    # We need to convert the luminosity of the star into real units now
    l_actual=luminosity*3.828E26
    print("The luminosity of the star in watts: "+str(l_actual))
    temperature=(l_actual/(sa*sb))**(1/4) # Solved Stefan-Boltzmann law for temperature
                                          # sb is the Stefan-Boltzmann constant
    print("The surface temperature of the star in Kelvin: "+str(temperature))

    ####### Time to produce a blackbody radiation curve #######

    T=temperature
    x=np.arange(1,1000)
    Intensity=((2*h*c**2)/(x**5))(1/(exp(((h*c)/(x*k*T)-1))
    plt.plot(x, Intensity, '--', label='Intensity')
    plt.legend(loc='upper right')
    plt.title("Planck Blackbody Radiation Curve")
    plt.xlabel("Wavelength")
    plt.ylabel("Intensity")

    plt.show()

Я пытаюсь построить кривую излучения черного тела для температуры данной звезды. Проблема, которую я получаю, ошибка . Код над функциями Plot полностью корректен и работает, я оставил его для справки. Я уверен, что это легко исправить, но я не знаю, как построить такие функции, как это? Все мои импортные и определенные переменные находятся над кодом.

Answer 1

Есть несколько проблем с

Intensity=((2*h*c**2)/(x**5))(1/(exp(((h*c)/(x*k*T)-1))

Скобки несбалансированы, отсутствует пропущенный * между закрытыми наборами скобок, а функция exp не определена. Поскольку x является np.ndarray, вам нужно использовать np.exp вместо math.exp, поэтому оно должно быть

Intensity = ((2 * h * c**2) / x**5) * (1 / (np.exp((h * c) / (x * k * T)) - 1))

Однако следует также отметить, что диапазон (1, 1000) довольно бесполезен в этой ситуации, поскольку x (длина волны) находится в метрах , а излучение в диапазоне длин волн от 1 м до 1 км не очень интересно. Это должно быть больше похоже на (1e-8, 1e-6) (от 10 нм до 1 мкм), которое может быть сгенерировано с помощью np.linspace(1e-8, 1e-6, 1000) и затем даст нормальные кривые черного тела:

---

Полный пример кода

import matplotlib.pyplot as plt
import numpy as np
import math

h = 6.62607015e-34
k = 1.380649e-23
c = 299792458.00
sb = (2 * math.pi**5 * k**4) / (15 * c**2 * h**3)

def one_star():
    mass = np.random.uniform(0.1, 50)
    radius = mass**0.9 * 6.96e8
    B = np.random.uniform(3.3, 3.7)
    sa = 4 * math.pi * radius**2
    T = (luminosity / (sa * sb))**(1/4)
    x = np.linspace(1e-8, 1e-6, 1000)
    I = ((2 * h * c**2) / x**5) * (1 / (np.exp((h * c) / (x * k * T)) - 1))
    plt.plot(x, I)

for i in range(5):
    one_star()
plt.show()