-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathblackbody.py
189 lines (153 loc) · 6.67 KB
/
blackbody.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
'''
blackbody.py - Color of thermal blackbodies.
Description:
Calculate the spectrum of a thermal blackbody at an arbitrary temperature.
Constants:
PLANCK_CONSTANT - Planck's constant, in J-sec
SPEED_OF_LIGHT - Speed of light, in m/sec
BOLTZMAN_CONSTANT - Boltzman's constant, in J/K
SUN_TEMPERATURE - Surface temperature of the Sun, in K
Functions:
blackbody_specific_intensity (wl_nm, T_K) -
Get the monochromatic specific intensity for a blackbody -
wl_nm = wavelength [nm]
T_K = temperature [K]
This is the energy radiated per second per unit wavelength per unit solid angle.
Reference - Shu, eq. 4.6, p. 78.
blackbody_spectrum (T_K) -
Get the spectrum of a blackbody, as a numpy array.
blackbody_color (T_K) -
Given a temperature (K), return the xyz color of a thermal blackbody.
Plots:
blackbody_patch_plot (T_list, title, filename) -
Draw a patch plot of blackbody colors for the given temperature range.
blackbody_color_vs_temperature_plot (T_list, title, filename) -
Draw a color vs temperature plot for the given temperature range.
blackbody_spectrum_plot (T_K) -
Draw the spectrum of a blackbody at the given temperature.
References:
Frank H. Shu, The Physical Universe. An Introduction to Astronomy,
University Science Books, Mill Valley, California. 1982. ISBN 0-935702-05-9.
Charles Kittel and Herbert Kroemer, Thermal Physics, 2nd edition,
W. H. Freeman, New York, 1980. ISBN 0-7167-1088-9.
License:
Copyright (C) 2008 Mark Kness
Author - Mark Kness - [email protected]
This file is part of ColorPy.
ColorPy is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation, either version 3 of
the License, or (at your option) any later version.
ColorPy is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with ColorPy. If not, see <http://www.gnu.org/licenses/>.
'''
import math, numpy, pylab
from . import colormodels
from . import ciexyz
from . import plots
# Physical constants in mks units
PLANCK_CONSTANT = 6.6237e-34 # J-sec
SPEED_OF_LIGHT = 2.997925e+08 # m/sec
BOLTZMAN_CONSTANT = 1.3802e-23 # J/K
SUN_TEMPERATURE = 5778.0 # K
def blackbody_specific_intensity (wl_nm, T_K):
'''Get the monochromatic specific intensity for a blackbody -
wl_nm = wavelength [nm]
T_K = temperature [K]
This is the energy radiated per second per unit wavelength per unit solid angle.
Reference - Shu, eq. 4.6, p. 78.'''
# precalculations that could be made global
a = (PLANCK_CONSTANT * SPEED_OF_LIGHT) / (BOLTZMAN_CONSTANT)
b = (2.0 * PLANCK_CONSTANT * SPEED_OF_LIGHT * SPEED_OF_LIGHT)
wl_m = wl_nm * 1.0e-9
try:
exponent = a / (wl_m * T_K)
except ZeroDivisionError:
# treat same as large exponent
return 0.0
if exponent > 500.0:
# so large that the final result is nearly zero - avoid the giant intermediate
return 0.0
specific_intensity = b / (math.pow (wl_m, 5) * (math.exp (exponent) - 1.0))
return specific_intensity
def blackbody_spectrum (T_K):
'''Get the spectrum of a blackbody, as a numpy array.'''
spectrum = ciexyz.empty_spectrum()
(num_rows, num_cols) = spectrum.shape
for i in range(0, num_rows):
specific_intensity = blackbody_specific_intensity (spectrum [i][0], T_K)
# scale by size of wavelength interval
spectrum [i][1] = specific_intensity * ciexyz.delta_wl_nm * 1.0e-9
return spectrum
def blackbody_color (T_K):
'''Given a temperature (K), return the xyz color of a thermal blackbody.'''
spectrum = blackbody_spectrum (T_K)
xyz = ciexyz.xyz_from_spectrum (spectrum)
return xyz
#
# Figures
#
def blackbody_patch_plot (T_list, title, filename):
'''Draw a patch plot of blackbody colors for the given temperature range.'''
xyz_colors = []
color_names = []
for Ti in T_list:
xyz = blackbody_color (Ti)
xyz_colors.append (xyz)
name = '%g K' % (Ti)
color_names.append (name)
plots.xyz_patch_plot (xyz_colors, color_names, title, filename)
def blackbody_color_vs_temperature_plot (T_list, title, filename):
'''Draw a color vs temperature plot for the given temperature range.'''
num_T = len (T_list)
rgb_list = numpy.empty ((num_T, 3))
for i in range(0, num_T):
T_i = T_list [i]
xyz = blackbody_color (T_i)
rgb_list [i] = colormodels.rgb_from_xyz (xyz)
# note that b and g become negative for low T - MatPlotLib skips those on the semilog plot.
plots.color_vs_param_plot (
T_list,
rgb_list,
title,
filename,
plotfunc = pylab.semilogy,
tight = True,
xlabel = r'Temperature (K)',
ylabel = r'RGB Color')
def blackbody_spectrum_plot (T_K):
'''Draw the spectrum of a blackbody at the given temperature.'''
spectrum = blackbody_spectrum (T_K)
title = 'Blackbody Spectrum - T %d K' % (int (T_K))
filename = 'BlackbodySpectrum-%dK' % (int (T_K))
plots.spectrum_plot (
spectrum,
title,
filename,
xlabel = 'Wavelength (nm)',
ylabel = 'Specific Intensity')
#ylabel = 'Intensity ($W/m^2$)') # with LaTex symbols, the axis text gets too big...
# Create sample figures
def figures ():
'''Create some blackbody plots.'''
# patch plots
T_list_0 = plots.log_interpolate ( 1200.0, 20000.0, 48)
T_list_hot = plots.log_interpolate (10000.0, 40000.0, 24)
T_list_cool = plots.log_interpolate ( 950.0, 1200.0, 24)
blackbody_patch_plot (T_list_0, 'Blackbody Colors', 'Blackbody-Patch')
blackbody_patch_plot (T_list_hot, 'Hot Blackbody Colors', 'Blackbody-HotPatch')
blackbody_patch_plot (T_list_cool, 'Cool Blackbody Colors', 'Blackbody-CoolPatch')
# color vs temperature
blackbody_color_vs_temperature_plot (range(1200, 16000, 50), 'Blackbody Colors', 'Blackbody-Colors')
blackbody_color_vs_temperature_plot (range(10000, 40000, 100), 'Hot Blackbody Colors', 'Blackbody-HotColors')
blackbody_color_vs_temperature_plot (range(950, 1200, 1), 'Cool Blackbody Colors', 'Blackbody-CoolColors')
# spectrum of specific temperatures
blackbody_spectrum_plot (2000.0)
blackbody_spectrum_plot (3000.0) # Proxima Centauri
blackbody_spectrum_plot (SUN_TEMPERATURE) # Sun
blackbody_spectrum_plot (11000.0) # Rigel
blackbody_spectrum_plot (15000.0)