7

I am trying to write a simple python code for a plot of intensity vs wavelength for a given temperature, T=200K. So far I have this...

import scipy as sp
import math
import matplotlib.pyplot as plt
import numpy as np
pi = np.pi
h = 6.626e-34
c = 3.0e+8
k = 1.38e-23

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

I don't know how to define wavelength(wav) and thus produce the plot of Plancks Formula. Any help would be appreciated.

dr jimbob
  • 17,259
  • 7
  • 59
  • 81
  • Presumably you want a range of wavelengths. Do you know what range? Look into [`range`](http://docs.python.org/2.7/library/functions.html#range). – jonrsharpe Mar 14 '14 at 23:22
  • 1
    Usually the wavelength is just even spaced vector, you can use `np.arange` or `np.linspace`. Remember the indent in your planck function. – Daniel Thaagaard Andreasen Mar 14 '14 at 23:24

3 Answers3

11

Here's a basic plot. To plot using plt.plot(x, y, fmt) you need two arrays x and y of the same size, where x is the x coordinate of each point to plot and y is the y coordinate, and fmt is a string describing how to plot the numbers.

So all you need to do is create an evenly spaced array of wavelengths (an np.array which I named wavelengths). This can be done with arange(start, end, spacing) which will create an array from start to end (not inclusive) spaced at spacing apart.

Then compute the intensity using your function at each of those points in the array (which will be stored in another np.array), and then call plt.plot to plot them. Note numpy let's you do mathematical operations on arrays quickly in a vectorized form which will be computationally efficient.

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

And you see:

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

You probably will want to clean up the axes labels, add a legend, plot the intensity at multiple temperatures on the same plot, among other things. Consult the relevant matplotlib documentation.

dr jimbob
  • 17,259
  • 7
  • 59
  • 81
  • 1
    why `math.e**b` instead of just `np.exp(b)` – ev-br Mar 14 '14 at 23:54
  • Because I took the function from `user3421850`'s question and they used `math.e**b`. – dr jimbob Mar 15 '14 at 04:04
  • @Zhenya - Also now that I look at it, their formula was incorrect, since they have an extra pi in the formula for `a`. (Note h is Planck's constant, not h-bar (Planck's reduced constant ħ = h/2pi).) – dr jimbob Mar 15 '14 at 04:27
  • Re units: it's anyway better to always use dimensionless units rather than (arbitrary) SI units. Write the formula on a piece of paper, use substitutions (here, x = \hbar * \omega / k T), express both rhs and lhs in terms of dimensionless functions of dimensionless variables, plot those (and use axis labels for stating the exact variables). – ev-br Mar 15 '14 at 10:52
  • @Zhenya - Why are you commenting to me about this? The `planck` equation was taken from OP's question only trivially modified initially (add indentation and remove parenthesis on return statement). There are many situations it makes more sense to plot in SI units; e.g., when engineering or when you have familiarity with wavelengths (e.g., visible is ~400nm (violet) - 700nm (red)). Furthermore you should never drop or unnecessarily add dimensionless constants like pi as you can never add them back in through dimensional analysis (bundling with another variable is fine). – dr jimbob Mar 15 '14 at 17:04
  • @Zhenya - The original question wanted intensity vs wavelength not frequency, so your substitution is something quite different. Note its B_\nu d\nu = - B_\lambda d\lambda, so the Blackbody formula plotted against frequency and wavelength will for example have peaks at different points (peak frequency != c/peak wavelength). Finally, T varies independently of wavelength (and frequency) and wavelength (and frequency) appear in other parts of the formula, so it doesn't make sense to bundle them together. (Doing so will make ugly/bizarre units.) But, feel free to make a plot of your own. – dr jimbob Mar 15 '14 at 17:14
  • `plt.hold(True)` functionality was [deprecated](https://github.com/matplotlib/matplotlib/issues/12337) in 2.1, and intentionally removed in 3.0. It will return the following error: `AttributeError: module 'matplotlib.pyplot' has no attribute 'hold'`. – DevonDahon Oct 05 '20 at 12:24
  • 1
    @DevonDahon - Matplotlib 2.1 with the deprecation was from October 2017. This comment predates that by three years. Will edit though. I believe just need to remove call to `hold()` as it's the default behavior. – dr jimbob Oct 06 '20 at 14:04
2

You may also want to use the RADIS library, which allows you to plot the Planck function against wavelengths, or against frequency / wavenumber, if needed !

from radis import sPlanck

sPlanck(wavelength_min=135, wavelength_max=3000, T=4000).plot()
sPlanck(wavelength_min=135, wavelength_max=3000, T=5000).plot(nfig='same')
sPlanck(wavelength_min=135, wavelength_max=3000, T=6000).plot(nfig='same')
sPlanck(wavelength_min=135, wavelength_max=3000, T=7000).plot(nfig='same')

Planck function calculated with RADIS for T=4000 K, 5000 K, 6000 K, 7000 K

erwanp
  • 1,262
  • 10
  • 19
0

Just want to point out that there seems to be an equivalent of what OP wants to do in astropy:

https://docs.astropy.org/en/stable/api/astropy.modeling.physical_models.BlackBody.html

Unfortunately, it is not very clear to me yet how to get wavelength vs frequency based expression.

Tobbey
  • 469
  • 1
  • 4
  • 18