Plancks Law, Frequency figures

Question

I want to plot the frequency version of planck's law. I first tried to do this independently:

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()

This did not work, I have an issue with scaling somehow. My next idea was to attempt to use this person's code from this question: Plancks Formula for Blackbody spectrum

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()

This didn't work, because I got a divide by zero error. Except that I don't see where there is a division by zero, the denominator shouldn't ever be zero since the exponential term shouldn't ever be equal to one. I chose the frequencies to be the conversions of the wavelength values from the example code.

Can anyone help fix the problem or explain how I can get planck's law for frequency instead of wavelength?

In the second example you don't appear to be plotting the values. May that be why it's failing with divide by zero? — blubberdiblub, Apr 15 '19 at 11:52
If I add a `plt.plot(frequencies, intensity4000)`, a `plt.gca().set_xlim([3e14, 1e15])` and a `plt.show()` to the second example, it works for me, despite some overflow warnings. The warnings indicate that you should have a good look at the ranges your numbers take on, though, as @Asmus correctly points out. — blubberdiblub, Apr 15 '19 at 13:06
I'd imagine part of the problem is that `frequency_range = np.linspace(10**-19,10**19,1000000)` should be `frequency_range = np.logspace(10**-19,10**19,1000000)`. Otherwise you'd be calculating `10**13` points. I general you're going to need to operate in logspace (i.e., find the log of intensity by carrying log through your calculations) becasue you're going to overflow your floats with exponents and cubes of numbers that big (or small). — Daniel F, Apr 15 '19 at 13:06

score 0 · Accepted Answer · answered Apr 15 '19 at 11:58

0

You can not safely handle such large numbers; even for comparably "small" values of b = h*freq/(k*T) your float64 will overflow, e.g np.exp(709.)=8.218407461554972e+307 is ok, but np.exp(710.)=inf. You'll have to adjust your units (exponents) accordingly to avoid this!

Note that this is also the case in the other question you linked to, if you insert print( np.exp(b)[:10] ) within the definition of planck(), you can examine the first ten evaluated b's and you'll see the overflow in the first few occurrences. In any case, simply use the answer posted within the other question, but convert the x-axis in plt.plot(wavelengths, intensity) to frequency (i hope you know how to get from one to the other) :-)

answered Apr 15 '19 at 11:58

Asmus

5,117
1
16
21

However, none of that will lead to a division by zero error with the values in the second example. It's more likely they just forgot to plot them or do further processing on the values that's not shown. – blubberdiblub Apr 15 '19 at 12:02
@blubberdiblub I suppose that may depend on their python/numpy versions, I get a `RuntimeWarning: overflow encountered in exp intensity = a/( (c**2 * (np.exp(b) - 1.0) ))` instead. – Asmus Apr 15 '19 at 12:09
The warning is normal with default settings. But it will happily continue using infinity as the value. Division by infinity is also ok with the default settings and yields 0.0 or -0.0, respectively, depending on the sign of a. – blubberdiblub Apr 15 '19 at 12:14
@blubberdiblub sorry, I should have phrased this differently. In either case, the covered range leads to problems; in the first example in question, we evaluate `brightness = a /(e_term - 1) `, which yields a division by zero due to the *underflow* in `np.exp(0.0) `. In the other case, we have an *overflow* due to `np.exp()` becoming `inf`. – Asmus Apr 15 '19 at 12:32
Yes, you are right, of course. And you're right that it's best to strive for a high quality of the equations and no warnings. My point is separate from that, though. In the question they said *"This didn't work, because I got a divide by zero error"* referring to the second example. However, I cannot see how it wouldn't **work** with the information provided, notwithstanding the overflow warning. The subset of the results that are affected by the overflow may or may not be useless, but it's not failing. I guess you could call it "unsafe" for the 0 results you get, or "safe" if 0 is acceptable. – blubberdiblub Apr 15 '19 at 12:46
@blubberdiblub now I get you, took me long enough :) You're absolutely right, both examples "worked", albeit with warnings. – Asmus Apr 15 '19 at 12:54

Plancks Law, Frequency figures

1 Answers1