Закон Планка, Частота цифры

Question

Я хочу построить частотную версию закона Планка.Сначала я попытался сделать это самостоятельно:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
%matplotlib inline

# Planck's Law
# Constants
h = 6.62607015*(10**-34) # J*s
c = 299792458 # m * s
k = 1.38064852*(10**-23) # J/K
T = 20 # K
frequency_range = np.linspace(10**-19,10**19,1000000)

def plancks_law(nu):
    a = (2*h*nu**3) / (c**2)
    e_term = np.exp(h*nu/(k*T))
    brightness = a /(e_term - 1)
    return brightness

plt.plot(frequency_range,plancks_law(frequency_range))
plt.gca().set_xlim([1*10**-16 ,1*10**16 ])
plt.gca().invert_xaxis()

Это не сработало, у меня проблема с масштабированием.Моей следующей идеей было попытаться использовать код этого человека из этого вопроса: Формула Планка для спектра черного тела

import matplotlib.pyplot as plt
import numpy as np

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

def planck_f(freq, T):
    a = 2.0*h*(freq**3)
    b = h*freq/(k*T)
    intensity =  a/( (c**2 * (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) 
frequencies = np.arange(3e14, 3e17, 1e14, dtype=np.float64) 

intensity4000 = planck_f(frequencies, 4000.)
plt.gca().invert_xaxis()

Это не сработало, потому что я получил ошибку деления на ноль.За исключением того, что я не вижу, где деление на ноль, знаменатель никогда не должен быть равен нулю, поскольку экспоненциальный член не должен быть равен единице.Я выбрал частоты для преобразования значений длины волны из примера кода.

Может кто-нибудь помочь решить проблему или объяснить, как я могу получить закон Планка для частоты вместо длины волны?

Answer 1

Вы не можете безопасно обрабатывать такие большие числа; даже для сравнительно «малых» значений b = h*freq/(k*T) ваш float64 будет переполнен, например, np.exp(709.)=8.218407461554972e+307 в порядке, но np.exp(710.)=inf. Чтобы избежать этого, вам нужно соответствующим образом настроить свои единицы (показатели)!

Обратите внимание, что это также относится и к другому вопросу , с которым вы связались; если вы вставите print( np.exp(b)[:10] ) в определении planck(), вы можете проверить первые десять оцененных b. и вы увидите переполнение в первых нескольких случаях. В любом случае, просто используйте ответ , опубликованный в другом вопросе, но преобразуйте ось x в plt.plot(wavelengths, intensity) в частоту (надеюсь, вы знаете, как добраться от одного к другому) :-)