2

I'm trying to graph a blackbody at different remperatures in the same graph, the temperatures are 10, 10^3, 10^5, 10^7, 10^9 and 10^12 kelvin. The 10^12 looks fine but the others look flat. I already tried to put it in log scale but it doesn't seem to help.

Here is my code:

import numpy as np 
import matplotlib.pylab as plt

k = 1.38 * 10**-23  

h = 6.62607015e-34

c = 3e8

def pla(T, nu):
    x = (h*nu) / (k*T)

    c1 = 2 * h * nu**3

    planck = c1 / (np.exp(x) * c**2)

    return planck

nus = np.linspace(0, 4e23)

T10 = pla(10, nus)
Te3 = pla(1e3, nus)
Te5 = pla(1e5, nus)
Te7 = pla(1e7, nus)
Te9 = pla(1e9, nus)
Te12 = pla(1e12, nus)

plt.xlabel("frecuencia")
plt.ylabel("temperatura")

plt.plot(nus, T10, color='blue', label='T10')
plt.plot(nus, Te3, color='red', label='Te3')
plt.plot(nus, Te5, color='orange', label='Te5')
plt.plot(nus, Te7, color='purple', label='Te7')
plt.plot(nus, Te9, color='gray', label='Te9')
plt.plot(nus, Te12, color='green', label='Te12')

plt.legend(loc='upper right')

plt.show()

result:

Hadus
  • 1,551
  • 11
  • 22
  • I suggest plotting one line at a time. Because some of them are flat. Also just simply print them out: `print(Te5)` (It is all zeros) – Hadus Feb 01 '20 at 21:42
  • Thanks, I just print them and they were 0s but I had seen the graphs for the others temperatures so changed nus=np.linspace(0, 4e23) for 4e20, 4e18, etc and the graphs and values were right, I don't get why is that – Jorge Amézquita Feb 02 '20 at 05:23
  • off topic comment but you could add/use colors related to temperatures of black body see [Star B-V color index to apparent RGB color](https://stackoverflow.com/a/22630970/2521214) to visually enhance the graph (either use as background or change the graph geometry completely) – Spektre Feb 02 '20 at 10:14

1 Answers1

0

Here is an implementation I found at: https://stackoverflow.com/a/22417613/6304086

import matplotlib.pyplot as plt
import numpy as np

h = 6.626e-34
c = 3.0e+8
k = 1.38e-23

def planck(wav, T):
    a = 2.0*h*c**2
    b = h*c/(wav*k*T)
    intensity = a/ ( (wav**5) * (np.exp(b) - 1.0) )
    return intensity

# generate x-axis in increments from 1nm to 3 micrometer in 1 nm increments
# starting at 1 nm to avoid wav = 0, which would result in division by zero.
wavelengths = np.arange(1e-9, 3e-6, 1e-9) 

# intensity at 4000K, 5000K, 6000K, 7000K
intensity4000 = planck(wavelengths, 4000.)
intensity5000 = planck(wavelengths, 5000.)
intensity6000 = planck(wavelengths, 6000.)
intensity7000 = planck(wavelengths, 7000.)

plt.plot(wavelengths*1e9, intensity4000, 'r-') 
# plot intensity4000 versus wavelength in nm as a red line
plt.plot(wavelengths*1e9, intensity5000, 'g-') # 5000K green line
plt.plot(wavelengths*1e9, intensity6000, 'b-') # 6000K blue line
plt.plot(wavelengths*1e9, intensity7000, 'k-') # 7000K black line

# show the plot
plt.show()

Blackbody at 4000K, 5000K, 6000K, 7000K as function of wavelength from 1nm to 3000 nm

Community
  • 1
  • 1
Hadus
  • 1,551
  • 11
  • 22
  • or a direction implementation using the RADIS library https://stackoverflow.com/a/65386199/5622825 – erwanp Dec 21 '20 at 00:34