In Python, you can manipulate your data in various ways. There are many, different functions that you can use in order to manipulate your data.
At the end of this module, we will be generating the Planck spectrum, which describes the intensity of a blackbody.
Before we generate it, we need to know how to read a file, and plot its content
There are a lot of modules to read in data from a file:
genfromtxt
loadfromtxt
read_csv
read_table
We will use "genfromtxt
" for this exercise
## First, you need to import your modules
%matplotlib inline
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
We can also define the path to the directory that contains all of our files
data_path = '../data/'
If you have questions about the function, you can add a "?" at the end of your function.
help(np.genfromtxt)
Help on function genfromtxt in module numpy.lib.npyio: genfromtxt(fname, dtype=<class 'float'>, comments='#', delimiter=None, skip_header=0, skip_footer=0, converters=None, missing_values=None, filling_values=None, usecols=None, names=None, excludelist=None, deletechars=None, replace_space='_', autostrip=False, case_sensitive=True, defaultfmt='f%i', unpack=None, usemask=False, loose=True, invalid_raise=True, max_rows=None, encoding='bytes') Load data from a text file, with missing values handled as specified. Each line past the first `skip_header` lines is split at the `delimiter` character, and characters following the `comments` character are discarded. Parameters ---------- fname : file, str, pathlib.Path, list of str, generator File, filename, list, or generator to read. If the filename extension is `.gz` or `.bz2`, the file is first decompressed. Note that generators must return byte strings in Python 3k. The strings in a list or produced by a generator are treated as lines. dtype : dtype, optional Data type of the resulting array. If None, the dtypes will be determined by the contents of each column, individually. comments : str, optional The character used to indicate the start of a comment. All the characters occurring on a line after a comment are discarded delimiter : str, int, or sequence, optional The string used to separate values. By default, any consecutive whitespaces act as delimiter. An integer or sequence of integers can also be provided as width(s) of each field. skiprows : int, optional `skiprows` was removed in numpy 1.10. Please use `skip_header` instead. skip_header : int, optional The number of lines to skip at the beginning of the file. skip_footer : int, optional The number of lines to skip at the end of the file. converters : variable, optional The set of functions that convert the data of a column to a value. The converters can also be used to provide a default value for missing data: ``converters = {3: lambda s: float(s or 0)}``. missing : variable, optional `missing` was removed in numpy 1.10. Please use `missing_values` instead. missing_values : variable, optional The set of strings corresponding to missing data. filling_values : variable, optional The set of values to be used as default when the data are missing. usecols : sequence, optional Which columns to read, with 0 being the first. For example, ``usecols = (1, 4, 5)`` will extract the 2nd, 5th and 6th columns. names : {None, True, str, sequence}, optional If `names` is True, the field names are read from the first line after the first `skip_header` lines. This line can optionally be proceeded by a comment delimeter. If `names` is a sequence or a single-string of comma-separated names, the names will be used to define the field names in a structured dtype. If `names` is None, the names of the dtype fields will be used, if any. excludelist : sequence, optional A list of names to exclude. This list is appended to the default list ['return','file','print']. Excluded names are appended an underscore: for example, `file` would become `file_`. deletechars : str, optional A string combining invalid characters that must be deleted from the names. defaultfmt : str, optional A format used to define default field names, such as "f%i" or "f_%02i". autostrip : bool, optional Whether to automatically strip white spaces from the variables. replace_space : char, optional Character(s) used in replacement of white spaces in the variables names. By default, use a '_'. case_sensitive : {True, False, 'upper', 'lower'}, optional If True, field names are case sensitive. If False or 'upper', field names are converted to upper case. If 'lower', field names are converted to lower case. unpack : bool, optional If True, the returned array is transposed, so that arguments may be unpacked using ``x, y, z = loadtxt(...)`` usemask : bool, optional If True, return a masked array. If False, return a regular array. loose : bool, optional If True, do not raise errors for invalid values. invalid_raise : bool, optional If True, an exception is raised if an inconsistency is detected in the number of columns. If False, a warning is emitted and the offending lines are skipped. max_rows : int, optional The maximum number of rows to read. Must not be used with skip_footer at the same time. If given, the value must be at least 1. Default is to read the entire file. .. versionadded:: 1.10.0 encoding : str, optional Encoding used to decode the inputfile. Does not apply when `fname` is a file object. The special value 'bytes' enables backward compatibility workarounds that ensure that you receive byte arrays when possible and passes latin1 encoded strings to converters. Override this value to receive unicode arrays and pass strings as input to converters. If set to None the system default is used. The default value is 'bytes'. .. versionadded:: 1.14.0 Returns ------- out : ndarray Data read from the text file. If `usemask` is True, this is a masked array. See Also -------- numpy.loadtxt : equivalent function when no data is missing. Notes ----- * When spaces are used as delimiters, or when no delimiter has been given as input, there should not be any missing data between two fields. * When the variables are named (either by a flexible dtype or with `names`, there must not be any header in the file (else a ValueError exception is raised). * Individual values are not stripped of spaces by default. When using a custom converter, make sure the function does remove spaces. References ---------- .. [1] NumPy User Guide, section `I/O with NumPy <http://docs.scipy.org/doc/numpy/user/basics.io.genfromtxt.html>`_. Examples --------- >>> from io import StringIO >>> import numpy as np Comma delimited file with mixed dtype >>> s = StringIO("1,1.3,abcde") >>> data = np.genfromtxt(s, dtype=[('myint','i8'),('myfloat','f8'), ... ('mystring','S5')], delimiter=",") >>> data array((1, 1.3, 'abcde'), dtype=[('myint', '<i8'), ('myfloat', '<f8'), ('mystring', '|S5')]) Using dtype = None >>> s.seek(0) # needed for StringIO example only >>> data = np.genfromtxt(s, dtype=None, ... names = ['myint','myfloat','mystring'], delimiter=",") >>> data array((1, 1.3, 'abcde'), dtype=[('myint', '<i8'), ('myfloat', '<f8'), ('mystring', '|S5')]) Specifying dtype and names >>> s.seek(0) >>> data = np.genfromtxt(s, dtype="i8,f8,S5", ... names=['myint','myfloat','mystring'], delimiter=",") >>> data array((1, 1.3, 'abcde'), dtype=[('myint', '<i8'), ('myfloat', '<f8'), ('mystring', '|S5')]) An example with fixed-width columns >>> s = StringIO("11.3abcde") >>> data = np.genfromtxt(s, dtype=None, names=['intvar','fltvar','strvar'], ... delimiter=[1,3,5]) >>> data array((1, 1.3, 'abcde'), dtype=[('intvar', '<i8'), ('fltvar', '<f8'), ('strvar', '|S5')])
Now we can go ahead and read in the data, and save it to two arrays, i.e. x1
and y1
:
x1, y1 = np.genfromtxt('../data/dataset1.txt',
unpack=True,
dtype=np.float)
We now use unpack
to tell Python to throw out the two columns and we caught them with arrays x
and y
, but we could have just captured whatever came out, then it just would be a merged array:
data = np.genfromtxt('../data/dataset1.txt', dtype=np.float)
You can now examine the output from genfromtxt
:
print(data.shape)
(100, 2)
print(data[:,0])
[-0.04826417 0.25723241 0.46084084 0.64903566 0.71677476 0.96776548 1.29983937 1.43387289 1.62135006 1.76059405 2.03743367 2.14958662 2.42060419 2.715836 2.9275557 2.9653225 3.19679942 3.4614748 3.54848877 3.76910315 3.98784208 4.19339365 4.4279586 4.62071517 4.79496323 5.00002693 5.23539664 5.37992987 5.60752551 5.93871238 6.08304306 6.3391359 6.56212174 6.70258117 6.84234717 7.01380281 7.22471216 7.42548426 7.58640829 7.93911613 8.08523514 8.27281698 8.50502081 8.58780099 8.97544342 9.06840145 9.32994708 9.4832187 9.7143636 9.9669831 10.05587739 10.24745965 10.50448738 10.64142437 10.97846567 11.13011018 11.30068135 11.57165861 11.65991901 11.94990211 12.18747631 12.4178335 12.55826718 12.82237611 12.85147113 13.04777332 13.42264142 13.46737614 13.81423191 13.97654878 14.1939695 14.37205959 14.53569893 14.66694151 14.91467005 15.22140463 15.36485266 15.59336199 15.69060866 15.89937533 16.23907517 16.36178348 16.54728757 16.73726968 17.04225118 17.24279728 17.33569323 17.56503833 17.70515795 18.02824083 18.23112027 18.33248662 18.62369577 18.77999396 18.89922539 19.10448361 19.3775723 19.56626863 19.83302568 20.05474468]
x1
array([-0.04826417, 0.25723241, 0.46084084, 0.64903566, 0.71677476, 0.96776548, 1.29983937, 1.43387289, 1.62135006, 1.76059405, 2.03743367, 2.14958662, 2.42060419, 2.715836 , 2.9275557 , 2.9653225 , 3.19679942, 3.4614748 , 3.54848877, 3.76910315, 3.98784208, 4.19339365, 4.4279586 , 4.62071517, 4.79496323, 5.00002693, 5.23539664, 5.37992987, 5.60752551, 5.93871238, 6.08304306, 6.3391359 , 6.56212174, 6.70258117, 6.84234717, 7.01380281, 7.22471216, 7.42548426, 7.58640829, 7.93911613, 8.08523514, 8.27281698, 8.50502081, 8.58780099, 8.97544342, 9.06840145, 9.32994708, 9.4832187 , 9.7143636 , 9.9669831 , 10.05587739, 10.24745965, 10.50448738, 10.64142437, 10.97846567, 11.13011018, 11.30068135, 11.57165861, 11.65991901, 11.94990211, 12.18747631, 12.4178335 , 12.55826718, 12.82237611, 12.85147113, 13.04777332, 13.42264142, 13.46737614, 13.81423191, 13.97654878, 14.1939695 , 14.37205959, 14.53569893, 14.66694151, 14.91467005, 15.22140463, 15.36485266, 15.59336199, 15.69060866, 15.89937533, 16.23907517, 16.36178348, 16.54728757, 16.73726968, 17.04225118, 17.24279728, 17.33569323, 17.56503833, 17.70515795, 18.02824083, 18.23112027, 18.33248662, 18.62369577, 18.77999396, 18.89922539, 19.10448361, 19.3775723 , 19.56626863, 19.83302568, 20.05474468])
You can check that the 1st column of data
is the same as x1
with np.array_equal
:
np.array_equal(x1, data[:,0])
True
Another nice thing about genfromtxt
is that it can read dta from a URL:
## Setting up path to remote file
remote_file = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.data.csv'
## Extracting data from file, and saving it as to variable `A`
A = np.genfromtxt(remote_file, unpack=True, delimiter=',')
Now A
has the shape of 9 columns by 768 rows
A.shape
(9, 768)
As you can see, the shape of A
is different than the one in the URL.
To fix it, you can transpose
the array:
A.T
array([[ 6. , 148. , 72. , ..., 0.627, 50. , 1. ], [ 1. , 85. , 66. , ..., 0.351, 31. , 0. ], [ 8. , 183. , 64. , ..., 0.672, 32. , 1. ], ..., [ 5. , 121. , 72. , ..., 0.245, 30. , 0. ], [ 1. , 126. , 60. , ..., 0.349, 47. , 1. ], [ 1. , 93. , 70. , ..., 0.315, 23. , 0. ]])
print(A.T.shape)
(768, 9)
Now that we've read the data, we can use to fit a straight line to it.
The steps to follow are:
myline
myline
¶def myline(x, m, b):
"""
Functional form of a straight line
Parameters
-----------
x : `float`, `int`, `np.ndarray`, list
Variable that tells you how far along
m : `float`, `int`
Slope or gradient
b : `float`, `int`
Y-intercept
Returns
---------
y : `float`, `int`, `np.ndarray`, list
Value for how far up on the y-axis
"""
y = (m * x) + b
return y
We can now fit a line to the data, and find the parameters (m
, and b
)
that best describe our data:
## Import `curve_fit` function from scipy
from scipy.optimize import curve_fit
## Calculating best-fit parameters
bestfit, covar_matrix = curve_fit(myline, x1, y1, p0 = [1.0, 1.0])
print(bestfit)
[3.12950289 3.72975489]
print("Best-fit parameters: m = {0} and b = {1}".format(*bestfit))
Best-fit parameters: m = 3.129502893478391 and b = 3.7297548918690517
In this example:
myline
: Is the function used to fit the datax1
, y1
: x- and y-valuesp0
: Initial guesses for the two parameters. This variable is optional.You can read more about curve_fit
by typing:
help(curve_fit)
Help on function curve_fit in module scipy.optimize.minpack: curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(-inf, inf), method=None, jac=None, **kwargs) Use non-linear least squares to fit a function, f, to data. Assumes ``ydata = f(xdata, *params) + eps`` Parameters ---------- f : callable The model function, f(x, ...). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments. xdata : An M-length sequence or an (k,M)-shaped array for functions with k predictors The independent variable where the data is measured. ydata : M-length sequence The dependent data --- nominally f(xdata, ...) p0 : None, scalar, or N-length sequence, optional Initial guess for the parameters. If None, then the initial values will all be 1 (if the number of parameters for the function can be determined using introspection, otherwise a ValueError is raised). sigma : None or M-length sequence or MxM array, optional Determines the uncertainty in `ydata`. If we define residuals as ``r = ydata - f(xdata, *popt)``, then the interpretation of `sigma` depends on its number of dimensions: - A 1-d `sigma` should contain values of standard deviations of errors in `ydata`. In this case, the optimized function is ``chisq = sum((r / sigma) ** 2)``. - A 2-d `sigma` should contain the covariance matrix of errors in `ydata`. In this case, the optimized function is ``chisq = r.T @ inv(sigma) @ r``. .. versionadded:: 0.19 None (default) is equivalent of 1-d `sigma` filled with ones. absolute_sigma : bool, optional If True, `sigma` is used in an absolute sense and the estimated parameter covariance `pcov` reflects these absolute values. If False, only the relative magnitudes of the `sigma` values matter. The returned parameter covariance matrix `pcov` is based on scaling `sigma` by a constant factor. This constant is set by demanding that the reduced `chisq` for the optimal parameters `popt` when using the *scaled* `sigma` equals unity. In other words, `sigma` is scaled to match the sample variance of the residuals after the fit. Mathematically, ``pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N)`` check_finite : bool, optional If True, check that the input arrays do not contain nans of infs, and raise a ValueError if they do. Setting this parameter to False may silently produce nonsensical results if the input arrays do contain nans. Default is True. bounds : 2-tuple of array_like, optional Lower and upper bounds on independent variables. Defaults to no bounds. Each element of the tuple must be either an array with the length equal to the number of parameters, or a scalar (in which case the bound is taken to be the same for all parameters.) Use ``np.inf`` with an appropriate sign to disable bounds on all or some parameters. .. versionadded:: 0.17 method : {'lm', 'trf', 'dogbox'}, optional Method to use for optimization. See `least_squares` for more details. Default is 'lm' for unconstrained problems and 'trf' if `bounds` are provided. The method 'lm' won't work when the number of observations is less than the number of variables, use 'trf' or 'dogbox' in this case. .. versionadded:: 0.17 jac : callable, string or None, optional Function with signature ``jac(x, ...)`` which computes the Jacobian matrix of the model function with respect to parameters as a dense array_like structure. It will be scaled according to provided `sigma`. If None (default), the Jacobian will be estimated numerically. String keywords for 'trf' and 'dogbox' methods can be used to select a finite difference scheme, see `least_squares`. .. versionadded:: 0.18 kwargs Keyword arguments passed to `leastsq` for ``method='lm'`` or `least_squares` otherwise. Returns ------- popt : array Optimal values for the parameters so that the sum of the squared residuals of ``f(xdata, *popt) - ydata`` is minimized pcov : 2d array The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters use ``perr = np.sqrt(np.diag(pcov))``. How the `sigma` parameter affects the estimated covariance depends on `absolute_sigma` argument, as described above. If the Jacobian matrix at the solution doesn't have a full rank, then 'lm' method returns a matrix filled with ``np.inf``, on the other hand 'trf' and 'dogbox' methods use Moore-Penrose pseudoinverse to compute the covariance matrix. Raises ------ ValueError if either `ydata` or `xdata` contain NaNs, or if incompatible options are used. RuntimeError if the least-squares minimization fails. OptimizeWarning if covariance of the parameters can not be estimated. See Also -------- least_squares : Minimize the sum of squares of nonlinear functions. scipy.stats.linregress : Calculate a linear least squares regression for two sets of measurements. Notes ----- With ``method='lm'``, the algorithm uses the Levenberg-Marquardt algorithm through `leastsq`. Note that this algorithm can only deal with unconstrained problems. Box constraints can be handled by methods 'trf' and 'dogbox'. Refer to the docstring of `least_squares` for more information. Examples -------- >>> import numpy as np >>> import matplotlib.pyplot as plt >>> from scipy.optimize import curve_fit >>> def func(x, a, b, c): ... return a * np.exp(-b * x) + c Define the data to be fit with some noise: >>> xdata = np.linspace(0, 4, 50) >>> y = func(xdata, 2.5, 1.3, 0.5) >>> np.random.seed(1729) >>> y_noise = 0.2 * np.random.normal(size=xdata.size) >>> ydata = y + y_noise >>> plt.plot(xdata, ydata, 'b-', label='data') Fit for the parameters a, b, c of the function `func`: >>> popt, pcov = curve_fit(func, xdata, ydata) >>> popt array([ 2.55423706, 1.35190947, 0.47450618]) >>> plt.plot(xdata, func(xdata, *popt), 'r-', ... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt)) Constrain the optimization to the region of ``0 <= a <= 3``, ``0 <= b <= 1`` and ``0 <= c <= 0.5``: >>> popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5])) >>> popt array([ 2.43708906, 1. , 0.35015434]) >>> plt.plot(xdata, func(xdata, *popt), 'g--', ... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt)) >>> plt.xlabel('x') >>> plt.ylabel('y') >>> plt.legend() >>> plt.show()
We can now overplot the best-fit line to the data:
# Initializing figure (optional)
fig = plt.figure(figsize=(10,10), facecolor='white')
ax = fig.add_subplot(111, facecolor='white')
# Plotting values
plt.plot(x1, myline(x1, *bestfit), 'r--', linewidth=2, label='Best fit')
plt.plot(x1, y1, 'bo', label='Data')
# Setting up limits
plt.xlim((-1, 21)) # Limit the x-range of our plot
# Axis labels
plt.xlabel('This is the x-label', fontsize=20)
plt.ylabel('This is the y-label', fontsize=20)
# Maybe add a title
plt.title('You can also add a title with color',
fontsize=20,
color='blue')
# And finally, a legend:
plt.legend(loc='best')
<matplotlib.legend.Legend at 0x151e747128>
The final script looks like this:
%load ../scripts/fitting_line.py
Now let's plot a distribution of values
## Importing modules
import numpy as np
import scipy
import matplotlib.pyplot as plt
Now we need to define the mean and standard deviation of a normal distribution, and create an array of random values:
# Mean and standard deviation
mu, sigma = 100, 15
# Array of random values
x = mu + (sigma * np.random.randn(10000))
# Printing out values of `x`
print(x)
[ 77.35324464 113.82141577 100.03311945 ... 111.53419448 135.45606453 93.43486279]
We can also define a function for the PDF of the distribution:
# Function for the PDF of distribution
def normpdf(x, mu, sigma):
"""
PDF of a distribution with a mean and standard deviation
Parameters
-----------
x : `np.ndarray`, list
List/Array of values of the distribution
mu : `float`, `int`
Mean of the distribution
sigma : `float`, `int`
Standard deviation of the distribution
Returns
--------
y : `np.ndarray` or list
List/array of the normalized PDF values
"""
u = (x-mu)/np.abs(sigma)
y = (1/(np.sqrt(2*np.pi)*np.abs(sigma)))*np.exp(-u*u/2)
return y
Now we construct a histogram with plt.hist
:
# Initialize figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, facecolor='white')
# Creating histogram
n, bins, patches = plt.hist(x,
bins=50,
density=True,
histtype='stepfilled',
facecolor='green',
alpha=0.75,
label='Normal distr.')
# Normalized PDF
y_pdf = normpdf(x, mu, sigma)
plt.plot(x, y_pdf, 'ro', linestyle='', label='PDF')
# Adding labels and title
plt.title(r'Histogram of IQ: $\mu = %s, \sigma = %s$' %(mu, sigma),
fontsize=20)
# Setting up axis
plt.axis([40, 160, 0, 0.03])
# Adding a grid
plt.grid(True)
# Adding legend
plt.legend(loc='best', prop={'size': 15})
<matplotlib.legend.Legend at 0x15202ffdd8>
The final script for this looks like:
%load ../scripts/histogram_pdf.py
The next step is to write a script that generates the Planck spectrum (wavelength and intensity at that wavelength for many wavelenghts)
Create an equation in Python syntax such that for temperature T=300 K, and wavelenth ($\lambda = 1cm$) it finds the intensity of a Planck Spectrum
Planck Spectrum:
$$ I = \frac{2hc^2}{\lambda^5}\frac{1}{e^{hc/(\lambda\ k T)} - 1} $$# Write your answer here
This method uses does the following:
## First import your modules
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
## Dictionary of constants
def const_dict_func():
"""
Dictionary of constants
Returns
---------
const_dict : `dict`
Dictionary of constants
"""
const_dict = {}
const_dict['c'] = 3.0e8 # Speed of light
const_dict['h'] = 6.62e-34 # Planck's constant
const_dict['k'] = 1.38e-23 # Boltzmann's constant
return const_dict
This will let you call the dictionary and get the values for each of the constants.
Keep in mind that the we had to know the units beforehand. There are other ways to do this, but we'll discuss them at the end.
def planck_spectrum(T, wavelen):
"""
Computes the Planck spectrum for a given Temperature and wavelength
Parameters
-----------
T : `float`, `int`
Temperature used. In units of `Kelvin`
wavelen : `float`, `int`
Wavelengths to evaluate
Returns
-----------
I : `np.ndarray`, `float`, `int`
Intensity of the Power spectrum at given temperature and wavelength
"""
## Constants
# Calling function `const_dict_func` and saving output as a dictionary
const_dict = const_dict_func()
# Saving values of constants as new variables
c = const_dict['c']
h = const_dict['h']
k = const_dict['k']
## Computing the Planck spectrum or 'radiance'
# Radiance
I = (2 * h * c ** 5) / (wavelen**5)
I *= 1./(-1. + np.exp((h * c)/(wavelen * k * T)))
return I
This function will plot the spectrum for the Planck spectrum as function of wavelength $\alpha$ at fixed Temperature $T$
## Plotting data
def plot_planck(data, T):
"""
Plots the Planck spectrum
Parameters
------------
data : `np.ndarray`
Data containing wavelength and planck spectrum information.
Shape is (N, 2)
T : `float`, `int`
Temperature, at which the Planck spectrum is analyzed.
"""
# Clearing previous figures
plt.clf()
plt.close()
# Initializing figure
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, facecolor='white')
# Plotting spectrum
plt.plot(data[:, 0], data[:, 1], marker='o', color='b', linestyle='--',
label='T={0} K'.format(T))
#
# Axis labels
ax.set_xlabel(r'$\lambda$ [m]', fontsize=20)
ax.set_ylabel('Intensity', fontsize=20)
# Legend
ax.legend(loc='best', prop={'size':20})
#
# Showing figure
plt.show()
def main():
"""
Main function
"""
#
# Temperature as an input parameter
T = 300 ## In units Kelvin
#
# Wavelenght
wavelen_min = 1.e-7 # Minimum wavelength in meters
waveln_max = 6.e-5 #Maximum wavelength in meters
samples = 100 # Number of steps to output
#
### List of wavelengths - Computing the wavelengths
# Creating array of number of samples
wavelen = np.arange(samples)
# Populating the array of wavelengths
wavelen = wavelen_min + (waveln_max - wavelen_min) * wavelen / float(samples)
#
# Computing the radiance
radiance = planck_spectrum(T, wavelen)
#
# Stacking the data
data = np.column_stack((wavelen, radiance))
# Sorting data from smallest to largest wavelength
data_sort_idx = np.argsort(data[:,0])
data_sort = data[data_sort_idx]
#
## Saving data to file
# Definiing output path and checking if it exists
output_path = '../datafiles'
if not (os.path.exists(output_path)):
os.makedirs(output_path)
# Saving data to file
np.savetxt('{0}/spectrum.dat'.format(output_path), data)
#
# Plotting data
plot_planck(data, T)
## Calling main function
main()
With all of this defined, you can put it in a single script.
This code is found under ../scripts/planck_spectrum.py
%load ../scripts/planck_spectrum.py
Now we can try to the Planck spectrum for multiple temperatures $T$.
We will have to modify the some of the functions to include $T$ as part of an argument.
We'll modify both the main
and plot_planck
functions.
We'll modify this function to:
T_arr
## Importing new modules
from matplotlib.pyplot import cm
## Plotting data
def plot_planck_2(data, T_arr):
"""
Plots the Planck spectrum
Parameters
------------
data : `np.ndarray`
Data containing wavelength and planck spectrum information.
Shape is ((N, 2), M)
T_arr : `np.ndarray`
Array of the different temperatures.
Shape, (M, ), i.e. it has `M` number of elements.
"""
# Defining the colormap
color_arr=cm.rainbow(np.linspace(0, 1, len(T_arr)))
# Clearing previous figures
plt.clf()
plt.close()
# Initializing figure
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, facecolor='white')
# Plotting spectrum
for kk, T_kk in enumerate(T_arr):
plt.plot(data[kk][:, 0],
data[kk][:, 1],
marker='o',
color=color_arr[kk],
linestyle='--',
label='{0} K'.format(T_kk))
#
# Axis labels
ax.set_xlabel(r'$\lambda$ [m]', fontsize=20)
ax.set_ylabel('Intensity', fontsize=20)
#
# Legend
ax.legend(loc='best', prop={'size':14})
#
# Showing figure
plt.show()
def main_2():
"""
Main function
"""
#
# Temperature as an input parameter
T_arr = np.arange(300, 1000, 100) ## In units Kelvin
#
# Wavelenght
wavelen_min = 1.e-7 # Minimum wavelength in meters
waveln_max = 6.e-5 #Maximum wavelength in meters
samples = 100 # Number of steps to output
#
### List of wavelengths - Computing the wavelengths
# Creating array of number of samples
wavelen = np.arange(samples)
# Populating the array of wavelengths
wavelen = wavelen_min + (waveln_max - wavelen_min) * wavelen / float(samples)
#
## Computing the radiance
# Defining array for `data`
data = [[] for x in range(len(T_arr))]
for kk in range(len(T_arr)):
radiance_kk = planck_spectrum(T_arr[kk], wavelen)
# Stacking the data
data_kk = np.column_stack((wavelen, radiance_kk))
# Sorting data from smallest to largest wavelength
data_kk_sort_idx = np.argsort(data_kk[:,0])
data_kk_sort = data_kk[data_kk_sort_idx]
# Saving to array
data[kk] = data_kk_sort
#
# Plotting data
plot_planck_2(data, T_arr)
main_2()
Now you see how the Planck spectrum changes as a function of wavelength $\lambda$!