This notebook forms part of a series on computational optical radiometry. The notebooks can be downloaded from Github. These notebooks are constantly revised and updated, please revisit from time to time.

This notebook was last tested with Python 3.6.1. Consider moving to Python 3.6 if you are not already using it.

Feedback is appreciated: neliswillers at gmail dot com.

```
In [1]:
```from IPython.display import display
from IPython.display import Image
from IPython.display import HTML

Consider the mathematical sine function $s(t) = A \sin(2\pi f t)$. Here $f$ is the frequency, $t$ is time and $A$ is the amplitude. As a mathematical function all three these parameters determine the value of the function. A sine wave generator on an engineer's workbench creates a temporal sine wave, i.e., a wave shape of sinusoidal characteristics, that varies with time. The engineer can adjust the properties of the signal by changing the frequency $f$ or amplitude $A$ by turning a knob on the instrument. The engineer cannot adjust the time $t$, because time is not a free variable in the sine wave generator, time flows as part of the human reality. We can say that the sine wave *depends* only on $f$ and $A$, but it *varies* with time $t$.

Imagine a pipe with water flowing through it. One can measure the volume of water flowing past a fixed point in a given time period; the longer the time period, the larger the water volume. It makes sense to define a temporal *flow rate* expressed e.g., in cubic metres per second (at given time $t$). Now generalise this notion to some form of *density* with units of [something/s], [something/m$^3$], or [something/$\mu$m]
These variables can be integrated over time/volume/spectral density to arrive at a cumulative quantity. The key observation here is that a quantity can be expressed as a density over time or volume, or in the radiometry case a quantity over a spectral band.

The spectral domain is similar to the time domain: it is not a free variable that can be adjusted, but rather provides a domain of observation against which an object's properties can be calculated or measured. These properties are said to *vary* with the spectral variable. Some spectral properties *depend* on free variables, e.g., Planck radiance *depends* on the temperature of the object, but *varies* with spectral variable.

Spectral domains (Sec 2.3.3) are commonly expressed in wavelength $\lambda$, wavenumber $\tilde{\nu}$ (the number of waves that will fit into a 1-cm length, written as one word) or frequency $f$. In the optical domain wavelength and wavenumber are more commonly used than frequency. Wavelength is commonly given in units of nanometres (nm) or micrometres ($\mu$m), wavenumber in (cm$^{-1}$) and frequency in [Hz].

Spectral variables come in two varieties: variable without density (equivalent to the sine wave) and density variables (equivalent to the water in the pipe). Variables without density include spectral transmittance of a filter, the atmosphere or emissivity. Variables with density include Planck radiation. Spectral density values have units of [something/$\mu$m], [something/cm$^{-1}$] or [something/Hz].

The total quantity of a spectral variable over a spectral range is determined by integrating the spectral variable over the spectral range, similar to the water example above.

Wavelength and wavenumber are related by $\lambda = 10^4/\tilde{\nu}$ (Sec 2.3.3). The conversion of spectral densities requires the derivative of the previous relation $d \tilde{\nu} = -d\lambda \,10^4/\lambda^2=-d\lambda\, \tilde{\nu}^2 /10^4$.

```
In [2]:
```%matplotlib inline
import numpy as np
#create a micrometre wavelength spectral range with 101 samples
#from 400 nm to 800 nm
wl1 = np.linspace(0.4, 0.8, 101) #rank (N,)
wl2 = wl1.reshape(-1,1) #convert to rank (N,1)
wl3 = wl2.reshape(1,-1) #convert to rank (1,N)
wl4 = wl2.reshape(-1,) #convert back to rank (N,)
print('Wavelength vector as rank one, has shape {}'.format(wl1.shape))
print('Wavelength vector as rank two, has shape {}'.format(wl2.shape))
print('Wavelength vector as rank two, has shape {}'.format(wl3.shape))
print('Wavelength vector as rank one, has shape {}'.format(wl4.shape))

```
```

```
In [3]:
```wl1 = np.linspace(0.4, 0.8, 101) #rank (N,)
wl2 = wl1.reshape(-1,1) #convert to rank (N,1)

`pyradi.ryutils.convertSpectralDomain`

and `pyradi.ryutils.convertSpectralDensity`

.

```
In [4]:
```import pyradi.ryutils as ryutils
import pyradi.ryplot as ryplot
wl = np.linspace(1, 14, 101).reshape(-1,1) # wavelength
wn = ryutils.convertSpectralDomain(wl, type='ln') # wavenumber
f = ryutils.convertSpectralDomain(wn, type='nf') # frequency
p = ryplot.Plotter(1,1,3,figsize=(12,4));
p.plot(1,wl,wn,'','Wavelength $\mu$m', 'Wavenumber cm$^{-1}$');
p.plot(2,wl,f,'','Wavelength $\mu$m', 'Frequency Hz');
p.plot(3,wn,f,'','Wavenumber cm$^{-1}$', 'Frequency Hz',maxNX=4);

```
```

`pyradi.ryplanck.planck`

is used to calculate the spectral density of Planck radiance (Sec 3.1, Sec D.4.1) for a temperature of 1000 K as a wavelength spectral density. This value is then converted to a wavenumber spectral density. Finally, the Planck radiance is calculated in the wavenumber domain, which should give the same result as the converted spectral radiance.

```
In [5]:
```import pyradi.ryplanck as ryplanck
wl = np.linspace(1, 14, 101)#.reshape(-1,1) # wavelength
wn = ryutils.convertSpectralDomain(wl, type='ln') # wavenumber
radiancewl = ryplanck.planck(wl,1000.0, 'el') / np.pi # W/(m2.sr.um)
_,radiancewn1 = ryutils.convertSpectralDensity(wl,radiancewl, 'ln')
radiancewn2 = ryplanck.planck(wn,1000, 'en') / np.pi
p = ryplot.Plotter(2,1,3,figsize=(12,4))
p.plot(1,wl, radiancewl,'Planck radiance','Wavelength $\mu$m',
'Radiance W/(m$^2$.sr.$\mu$m)',maxNX=4)
p.plot(2,wn, radiancewn1,'Planck radiance','Wavenumber cm$^{-1}$',
'Radiance W/(m$^2$.sr.cm$^{-1}$)',maxNX=4)
p.plot(3,wn, radiancewn2,'Planck radiance','Wavenumber cm$^{-1}$',
'Radiance W/(m$^2$.sr.cm$^{-1}$)',maxNX=4);

```
```

See the section Performing spectral integrals below on performing integrals on spectral data.

SciPy has a number of interpolation options in `scipy.interpolate`

. The one-dimensional interpolation scheme enables the definition of the lookup table as a function by using `scipy.interpolate.interp1d`

. The user can define the interpolation scheme to be used - the default is `kind='linear'`

, to do cubic interpolation use `kind='cubic'`

. In the example below observe the difference between the two schemes.

The example below is is an ill-posed problem: normally one would have a much denser input data set. The purpose with using this poor data set it to illustrate the dangers of interpolation. Do not use a potentially unstable interpolation scheme such as cubic interpolation on a poor data set such as the one shown below

The input data set has only four samples. The ordinate (y) data has a slight oscillation. The abscissa (x) data are defined in the wavelength space. Equivalently, the abscissa (x) data are also defined in the wavenumber space, because $\tilde{\nu}=10^4/\lambda$. This is the only information available, there is no requirement on the values between the stated abscissa.

Important: *by selecting the interpolation scheme (linear, cubic, etc.), you are implicitly making a modelling assumption*. Linear interpolation in the wavelength space means a non-linear interpolation in wavenumber space and vice versa (because wavelength and wavenumber are nonlinearly related). Stated differently: if your input data set was compiled to be interpolated linearly in wavelength space (blue lines), then linear interpolation in wavenumber space will give the wrong results (red lines).

Also, by *selecting between the linear and cubic (or any other scheme) you are also making a modelling assumption*. Compare the difference between the results for the linear interpolation (red/blue) and the cubic interpolation (green/cyan) lines.

It is evident in the cubic interpolation graphs (green/cyan) that cubic/high-order interpolation schemes do not take kindly to large variations in abscissa increments. For example, the large increment from 10 to 20 in the wavelength space (green line) resulted in a huge variation in the spline curve. Likewise, the large wavenumber increment from 5000 to 10000 (or 1 to 2 $\mu$m) (cyan line) causes a large variation.

From this experiment we can draw the following conclusions:

Clearly state the space (wavelength/wavenumber) where the input data set applies. Do interpolation

*only*in the prescribed space.Ensure that the input data set samples at similar/even density throughput the abscissa domain. If no input data is available in some subranges in the abscissa domain, construct synthetic data by some modelling assumptions, to ensure even density. Such an even spread in abscissa values will limit wide swings in the high-order interpolation schemes.

Thanks to Riana Willers for pointing out the subtle matter of interpolation in wavelength or wavenumber space.

```
In [6]:
```from scipy import interpolate
import pyradi.ryplot as ryplot
import numpy as np
%matplotlib inline
# create input data in wavelength and wavenumber domain
maxx = 24
wl = np.asarray([1, 2, 5, 10, maxx])
wn = 1e4 / wl # wl and wn values precisely align
fn = np.asarray([2,3,1, 2, -1])
# create the plotting samples - much tighter than input data
wli = np.linspace(1,maxx,400)
wni = 1e4 / wli
#create the lookup functions
fwll = interpolate.interp1d(wl,fn)
fwlc = interpolate.interp1d(wl,fn, kind='cubic')
fwnl = interpolate.interp1d(wn,fn)
fwnc = interpolate.interp1d(wn,fn, kind='cubic')
p = ryplot.Plotter(1,2,2,figsize=(12,8))
p.plot(1,wli,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(1,wli,fwnl(wni),label=['wn scale interpolation - linear'],plotCol=['r'], maxNX=5)
p.resetPlotCol()
p.plot(2,wli,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(2,wli,fwlc(wli),label=['wl scale interpolation - cubic'],plotCol=['g'])
p.plot(2,wli,fwnc(wni),label=['wn scale interpolation - cubic'],plotCol=['c'], maxNX=5)
p.resetPlotCol()
p.plot(3,wni,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(3,wni,fwnl(wni),label=['wn scale interpolation - linear'],plotCol=['r'], maxNX=5)
p.resetPlotCol()
p.plot(4,wni,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(4,wni,fwlc(wli),label=['wl scale interpolation - cubic'],plotCol=['g'])
p.plot(4,wni,fwnc(wni),label=['wn scale interpolation - cubic'],plotCol=['c'], maxNX=5);

```
```

In most real-world system analyses, filter spectral responses will be measured and used in tabular formats. However, in some instances it is convenient to create spectral filter shapes according to a parameter set. The pyradi library has a function `pyradi.ryutils.sfilter`

to assist in calculating a variety of spectral filter shapes. This function has the form (Sec D.4.2):

This function takes a spectral variable $\lambda$, the centre wavelength $\lambda_c$, the spectral width $\Delta\lambda$, the sharpness of cutoff $s$ (if $s=2$, the curve has a Gaussian shape, and if $s=\infty$, the curve has a square top-hat shape), transmittance in the pass band $\tau_p$ and stop band $\tau_s$ and the type of filter. The use of the function is best demonstrated by graphical means.

```
In [7]:
```wavelength = np.linspace(0.1, 1.3, 350).reshape(-1, 1)
#demonstrate the different filter types
width = 0.4
center = 0.9
filterExp=[2, 2, 4, 6, 8, 12, 500] # shaprness
filterPass=[0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99] # tau pass
filterSupp = np.asarray(filterPass) * 0.1 # tau stop
filterType=['bandpass', 'lowpass', 'highpass', 'bandpass', 'lowpass',
'highpass', 'bandpass']
filterTxt = [r's={0}, $\tau_p$={2}, {1}'.format(s,y,z) for s,y,
z in zip(filterExp, filterType, filterPass) ]
filters = ryutils.sfilter(wavelength,center, width, filterExp[0],
filterPass[0], filterSupp[0], filterType[0])
for i,exponent in enumerate(filterExp[1:]):
tau=ryutils.sfilter(wavelength,center, width, filterExp[i+1],filterPass[i+1],
filterSupp[i+1], filterType[i+1])
filters = np.hstack((filters,tau))
#plot sample filters
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4));
smpleplt.plot(1, wavelength, filters,
r"Optical filter for $\lambda_c$={0}, $\Delta\lambda$={1}".format(center,width),
r'Wavelength [$\mu$m]', r'Transmittance', \
['r', 'g', 'y','g', 'b', 'm','c'],label=filterTxt);

```
```

```
In [8]:
```#all passband, different shapes
width = 0.5
center = 0.7
filterExp=[2, 4, 6, 8, 12, 500] # sharpness
filterTxt = ['s={0}'.format(s) for s in filterExp ]
filters = ryutils.sfilter(wavelength,center, width, filterExp[0], 0.8, 0.1)
for exponent in filterExp[1:]:
filters = np.hstack((filters, ryutils.sfilter(wavelength,center, width,
exponent, 0.8, 0.1)))
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4));
smpleplt.plot(1, wavelength, filters,
r"Optical filter for $\lambda_c$=0.7,$\Delta\lambda$=0.5,$\tau_{s}$=0.1, $\tau_{p}$=0.8",
r'Wavelength [$\mu$m]', r'Transmittance', \
['r', 'g', 'y','g', 'b', 'm'],label=filterTxt);

```
```

Detector spectral responsivity values are also normally measured or taken from detector data sheets. In the absence of such information the shape can be calculated from a detailed bulk detector model, such as provided in `pyradi.rydetector`

, but sometimes a simplified shape is required. Such a simple photon detector responsivity can be calculated with the
`pyradi.ryutils.responsivity`

function. The function takes the form (Sec D.4.3):

This function takes a wavelength spectral variable $\lambda$, the peak wavelength value $\lambda_c$, the sharpness of cuton (lower wavelength range $a$, the sharpness of cutoff (longer wavelength range) $n$, and an amplitude scaling factor $k$. Best results are obtained for $0\leq a \leq 5$ and $5\leq n \leq 50$. The use of the function is best demonstrated by graphical means.

```
In [9]:
```lwavepeak = 1.2
params = [(0.5, 5), (1, 10), (1, 20), (1, 30), (1, 1000), (2, 20)] #cut-on and cut-off
parameterTxt = ['a={0}, n={1}'.format(s[0], s[1]) for s in params ]
responsivities = ryutils.responsivity(wavelength,lwavepeak, params[0][0], params[0][1], 1.0)
for param in params[1:]:
responsivities = np.hstack((responsivities, ryutils.responsivity(wavelength,lwavepeak, param[0], param[1], 1.0)))
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4))
smpleplt.plot(1, wavelength, responsivities, "Detector Responsivity for $\lambda_c$=1.2 $\mu$m, k=1", r'Wavelength [$\mu$m]',\
r'Responsivity', \
['r-', 'g-', 'y-','g--', 'b-', 'm-'],label=parameterTxt);

```
```

```
In [10]:
```# here we simply multiply the responsivity and spectral filter spectral curves.
# this is a silly example, but demonstrates the spectral integral.
filtreps = responsivities * filters
parameterTxt = [str(s)+' & '+str(f) for (s, f) in zip(params, filterExp) ]
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4))
smpleplt.plot(1, wavelength, filtreps, "Filtered Detector Responsivity",
r'Wavelength $\mu$m',\
r'Responsivity', \
['r-', 'g-', 'y-','g--', 'b-', 'm-'],label=parameterTxt,legendAlpha=0.5);

```
```

`linspace`

function), but plotted in both wavelength and wavenumber domains. Note that the data was not resampled to wavenumber domain before plotting: the wavenumber plots are not at constant intervals. The wavelength domain was merely recalculated to wavenumber domain. This is illustrated in the graph below.

```
In [11]:
```wlc = np.linspace(0.5, 8.5, 17).reshape(-1,1) # wavelength
wnc = ryutils.convertSpectralDomain(wlc, type='ln') # wavenumber
p = ryplot.Plotter(2,1,3,figsize=(18,2))
for i in range(wlc.shape[0]):
p.plot(1,np.array([wlc[i],wlc[i]]), np.array([0,1.]), 'Constant wavelength intervals',
'Wavelength $\mu$m', plotCol=['r'], markers=['o'], maxNY=1)
p.plot(2,np.array([wnc[i],wnc[i]]), np.array([0,1.]), 'Non-constant wavenumber intervals',
'Wavenumber cm$^{-1}$',plotCol=['r'], markers=['o'], maxNX=5, maxNY=1, pltaxis=[0, 10000, 0, 1]);

```
```

The data and the plot are equally 'valid' or true in either domain, because the domain merely samples the data at different spectral intercepts. So we can load and process data in either domain and plot in either domain, without losing any truth value. Plot in whatever domain illustrates your argument the best.

You should however take care when dealing with spectral densities: it is possible and mathematically allowable to plot a wavenumber density with units of W/(m$^2$.sr.cm$^{-1}$) against a wavelength domain with units of $\mu$m, but it does not make much practical sense to do so.

`numpy.diff`

. Note that the intervals are constant and positive for the wavelength scale, but negative and non-constant for the wavenumber scale as converted from the wavelength scale. Recall from Section Spectral-variables above, that $d \tilde{\nu} = -d\lambda \,10^4/\lambda^2=-d\lambda\, \tilde{\nu}^2 /10^4$. In other words the spectral intervals change in the opposite sense: the wavelength domain scale increases (positive intervals), but when converted to wavenumber domain the scale decreases (negative intervals).

```
In [12]:
```p = ryplot.Plotter(3,1,3,figsize=(18,2))
p.plot(1,(wlc[1:]+wlc[:-1])/2, np.diff(wlc,axis=0), 'Constant wavelength intervals',
'Wavelength $\mu$m')
p.plot(2,(wnc[1:]+wnc[:-1])/2, np.diff(wnc,axis=0), 'Non-constant wavenumber intervals',
'Wavenumber cm$^{-1}$');

```
```

`ryfiles.loadColumnTextFile`

).

A repeating pattern in computational radiometry is the spectral integral: integrating the sum of a spectral variable between an upper and a lower limit. For example, this calculation is required to determine the total inband radiance over a wide spectral range. We use the mathematical integral notation to indicate such integrals, but in practice, these are calculated using numerical summation. $$ L = \int_{\lambda_0}^{\lambda_1} L_\lambda d\lambda = \sum_{i=0}^{N}L(i)\Delta\lambda(i) $$ where the sample $i=0$ corresponds with the sample at $\lambda_0$ and the sample $i=N$ corresponds with the sample at $\lambda_1$ The sum above is a general statement and does not make any assumption about the spectral domain sampling interval. If the spectral domain sampling interval has fixed increments $\Delta\lambda(i)=\Delta\lambda$ then $$ L = \int_{\lambda_0}^{\lambda_1} L_\lambda d\lambda = \Delta\lambda\sum_{i=0}^{N}L(i) $$

In pyradi the spectral integrals are typically (and naively) calculated using the `numpy.trapz`

trapezoidal integration function. This function has the signature

`numpy.trapz(y, x=None, dx=1.0, axis=-1)`

where `y`

is the value to be integrated (i.e., spectral density $L_\lambda$), `x`

is the domain where the value `y`

is sampled(i.e., wavelength $\lambda$), `dx`

is the *constant* $\Delta\lambda$ above, and `axis`

is the axis along which integration is done.

If the spectral domain is sampled at regular intervals the `x`

array is not not required and unless `dx`

is specified, it is assumed that `dx`

=$\Delta\lambda$=1. In this case the spectral sample interval must be accounted for, either by (1) pre-multiplication, (2) by stating the value for `dx`

or by (3) supplying the value for `x`

:

```
In [13]:
```wl = np.linspace(1, 14, 1001)#.reshape(-1,1) # wavelength
temperature = 1000.0 # K
radiancewl = ryplanck.planck(wl,temperature, 'el') / np.pi
print('Integral(1) = {} W/(m^2.sr)'.format( (wl[1]-wl[0]) * np.trapz(radiancewl)) )
print('Integral(2) = {} W/(m^2.sr)'.format( np.trapz(radiancewl, dx=(wl[1]-wl[0]))) )
print('Integral(3) = {} W/(m^2.sr)'.format( np.trapz(radiancewl,wl)) )

```
```

`dx`

or pre-multiplication approach is no longer valid. The only possibility to calculate the integral is to give the value of `x`

, with its non-constant intervals. Continuing from the previous example, convert the wavelength domain to the wavenumber domain, calculate the spectral density in wavenumber units W/(m$^2$.sr.cm$^{-1}$) and calculate the integral using the (non-constant-interval) wavenumber domain as `x`

:

```
In [14]:
```wn = ryutils.convertSpectralDomain(wl, type='ln') # wavenumber
radiancewn = ryplanck.planck(wn, temperature, 'en') / np.pi
print('Integral(wn) = {} W/(m^2.sr)'.format( np.trapz(radiancewn, wn )))

```
```

What happened here? The integral appears to have (practically) the same absolute values, but the sign is wrong!

Recall from Section Spectral-variables above, that $d \tilde{\nu} = -d\lambda \,10^4/\lambda^2=-d\lambda\, \tilde{\nu}^2 /10^4$. We can see this by taking the first few elements of the `numpy.diff`

along the spectral domains. Notice also that the wavelength interval is fixed, but the wavenumber interval varies.

```
In [15]:
```print('First 5 intervals in wl: {}'.format(np.diff(wl)[:5]))
print('First 5 intervals in wn: {}'.format(np.diff(wn)[:5]))
#now calculate the intervals using the equations above:
print('First 5 intervals in wn: {}'.format( - np.diff(wl)[:5] * 1.0e4 / ((wl[:5] + wl[1:6])/2. )**2 ))
print('First 5 intervals in wn: {}'.format( - np.diff(wl)[:5] * ((wn[:5] + wn[1:6])/2. )**2 / 1.0e4))

```
```

```
In [16]:
```print('Integral(wn) = {} W/(m^2.sr)'.format( - np.trapz(radiancewn, wn )))
print('Integral(wn) = {} W/(m^2.sr)'.format( np.abs( np.trapz(radiancewn, wn ))))

```
```

`numpy.flipud`

) both the spectral density and spectral domain variables, to ensure matching elements.

```
In [17]:
```#flip only wavenumber but not the wavenumber spectral density
print('Integral(wn) = {} W/(m^2.sr) - wrong because of a mismatch between x and y ordering'.\
format( np.trapz(radiancewn, np.flipud(wn) )))
#flip both wavenumber and wavenumber spectral density
print('Integral(wn) = {} W/(m^2.sr)'.format( np.trapz(np.flipud(radiancewn), np.flipud(wn) )))

```
```

The relationship between the photon rate radiometric units and radiant units are well defined for a single wavelength. It is best to calculate the spectral photon rate units from first principles and then integrate over the spectral band of interest. However on occasion the information may only be available in wideband radiant units, but the quantity must be expressed in photon rate units. For wideband data the question then arises which wavelength should be used in the conversion? This section determined the appropriate wavelength to be used for the conversion for a given spectral band and source spectrum.

The appropriate conversion wavelength is determined at the wavelength where an accurate spectral calculation gives the same answer as an approximate wideband calculation. The experimentation below indicates that the conversion wavelength is slightly longer than the center of the spectral band - photons at longer wavelengths have lower energy than at shorter wavelengths.

Responsivity is given by \begin{equation} R_e = \frac{\eta q \lambda}{hc} \end{equation} from which quantum efficiency can be determined \begin{equation} \eta = \frac{R_e hc}{q \lambda} \end{equation}

```
In [18]:
```# to investigate the conversion for broadband data between photon rate units and radiant units
import scipy.constants as const
spec = np.loadtxt('data/InSb.txt')
wl = spec[:,0].reshape(-1)
respW = 2. * spec[:,2].reshape(-1)
respQ = respW * const.h * const.c / (wl * 1e-6 * const.e)
wlNan = np.where(respW==0,np.nan,wl)
wlcent = np.mean([np.nanmax(wlNan), np.nanmin(wlNan)])
print('wl min = {}'.format(np.nanmin(wlNan)))
print('wl mid = {}'.format(wlcent))
print('wl max = {}'.format(np.nanmax(wlNan)))
p = ryplot.Plotter(1,2,2,figsize=(10,6));
p.plot(1, wl,respW, 'Radiant domain','Wavelength $\mu$m','Responsivity A/W');
p.plot(2, wl,respQ,'Wavelength $\mu$m','Quantum efficiency');
q = ryplot.Plotter(2,1,1,figsize=(8,4))
for temp in [250.,300.,350.]:
Ve = respW * ryplanck.planck(wl,temp,'el') / np.pi
Vq = respQ * ryplanck.planck(wl,temp,'ql') / np.pi
Vei = np.trapz(Ve,-wl,axis=0)
Vqi = np.trapz(Vq,-wl,axis=0)
# convert from W to q/s
qW = const.h * const.c / (wl * 1e-6)
ratio = ((Vei/qW) / Vqi ) / (np.max(respW)/np.max(respQ))
wlu = np.interp(1.,ratio[::-1],wl[::-1])
q.plot(1, wl,ratio, 'Ratio (normalised responses) as function of conversion wavelength',
'Wavelength $\mu$m', 'Ratio',
label=['{} K, crossover at {:.3f} $\mu$m'.format(temp,wlu)]);
p.plot(3, wl,Ve, 'Radiant domain',
'Wavelength $\mu$m', 'A/(m$^2$.sr.$\mu$m)',
label=['{} K, {:.3e} A/(m$^2$.sr)'.format(temp,Vei)]);
p.plot(4, wl,Vq, 'Quantum domain ',
'Wavelength $\mu$m', 'q/(s.m$^2$.sr.$\mu$m)',
label=['{} K, {:.3e} q/(s.m$^2$.sr)'.format(temp,Vqi)]);
# p.plot(4, wl,Ve/const.e);

```
```

The effective value of a spectral variable is given by (Sec 7.2.2)

$$ {F}_{\rm eff} = \frac {\int_{0}^{\infty}{F} {G} d\lambda} {\int_{0}^{\infty}{G} \lambda}, $$where ${F}$ is the variable in question, and ${G}$ is a weighting function. Note that the effective value of $F$ depends on the shapes of both $F$ and $G$; the effective value of $F$ thus calculated therefore applies only to the specific weighting function $G$ - it is incorrect for any other shapes.

To demonstrate the sensitivity of the effective value normalisation to the spectral variability of the input spectra, consider the effective detector responsivity calculated at different blackbody source temperatures from 500 K to 2000 K. The InSb detector is modelled using the simple detector model in `ryutils.responsivity`

. The source radiance is the weighting function $G$, calculated by using `ryplanck.planck`

. The following graphs show the spectral detector responsivity ($F$ in the above equation) and the effective responsivity as calculated for the range of temperatures. The calculation requires the calculation of a separate quotient and integral for each temperature, yet there are no loops! This magic is achieved by constructing two-dimensional arrays with wavelengths along the row (0'th) axis and temperature along the column (1'st) axis. The first of these arrays is returned by the `ryplanck.planck`

function (because a vector of wavelengths and a vector of temperatures are passed to the function). The second array is obtained by using the `numpy.meshgrid`

array (google to see how this function works). The two arrays are used to calculate the numerator of the quotient. The integrals are implemented using the `numpy.trapz`

function, integrating along the wavelength axis=0 direction. This integration returns a vector of values corresponding to the temperature values along the columns axis=1.

Note the very wide variation in effective responsivity when calculated with different source temperatures. This stems from the fact that the responsivity is weighted differently by the different Planck radiance curves resulting from the different source temperatures. The highest effective responsivity occurs near the wavelength where the peak of the Planck radiance matches the peak of the detector responsivity.

```
In [19]:
```wavelength = np.linspace(0.1, 20, 305).reshape(-1, 1)
detector = ryutils.responsivity(wavelength, 6.1, 1, 15, 1.0)
temperatures = np.linspace(500, 2000, 41)
spectralBaseline = ryplanck.planck(wavelength,temperatures,'el')
_, detVar = np.meshgrid(temperatures, detector)
spectralWeighted = spectralBaseline * detVar
effResp = np.trapz(spectralWeighted, wavelength, axis=0 ) / np.trapz(spectralBaseline, wavelength, axis=0)
p = ryplot.Plotter(2,1,3,figsize=(18,6))
p.plot(1,wavelength,detector,'Spectral responsivity','Wavelength $\mu$m', 'Responsivity A/W',pltaxis=[0.5,7,0,0.8])
p.plot(2,temperatures,effResp,'Effective responsivity','Temperature K', 'Responsivity A/W')

```
Out[19]:
```

Data are easily read in from files. The format most easy to read/write is ASCII coding in a tabular format (also known as a flat file). ASCII tabular files can use commas (csv format) or whitespace (tabs, spaces) to delimit the columns from each other. Python provides specialist CSV readers or just ASCII format readers. The data processing package Pandas also has very powerful data reading and manipulation facilities, especially for data in complex structural formats. In the pyradi context we are more interested in opening files with numpy and this analysis will focus on reading and writing files for Numpy applications.

Numpy (in versions later than 1.6) has powerful functions to read and write data files to and from Numpy arrays. The file formats include ASCII, binary format and compressed binary formats.

ASCII flat files are convenient data stores because these file can be read/written by a text editor and with a little effort with Excel (import data from txt, save as txt). In this example an array is created, written to file and read back from the file. When writing the text file, the string format is set to truncate the number of decimal values - on purpose, because we want to compare the array in memory with the array read in from file. The `numpy.loadtxt`

and `numpy.savetxt`

functions does the work. These functions take as a parameter the filename or file handle (if the file is opened elsewhere). If the filename has the extension `.gz`

, the file is handled as a gzip compressed file.

When reading and writing files, the `with`

context manager is used. The motivation for using `with`

is that if something goes wrong in the opening, reading or writing, the file is closed in an orderly manner. Also, it is more concise and requires less typing, because the `with`

context manager for the file closes the file automatically when leaving the block. Note that the `with`

context manager provides a file handle: the `gz`

extension is lost and the file is not compressed.

In this example the file is written a two-line header, started with the comment symbol. When the file is read again later, the lines with the comment symbol are ignored. The last few lines read only selected columns - useful if you do not need all the data in the file.

```
In [20]:
```arr = np.random.randn(3,4)
print('original data:\n{}'.format(arr))
#write/read uncompressed file
filename = 'arrfile.txt'
with open(filename, 'wb' ) as filehandle:
np.savetxt(filehandle, arr, fmt='%12.5e',
header='# this is a header\n#Line two of header', comments='#')
with open(filename, 'rb' ) as filehandle:
arr2 = np.loadtxt(filehandle, comments='#')
print('read from text file {} data:\n{}'.format(filename, arr2))
#write/read gzip compressed file - does not work with `with`
filename = 'arrfile.gz'
np.savetxt(filename, arr, fmt='%12.5e',
header='# this is a header', comments='#')
arrz = np.loadtxt(filename, comments='#')
print('read from text file {} data:\n{}'.format(filename, arrz))
usecols=[0,3]
arrz = np.loadtxt(filename, usecols=usecols, comments='#')
print('read columns {} from text file {} data:\n{}'.format(usecols,
filename, arrz))

```
```

Data can also be written in binary format or compressed binary format. These files can only be read by the appropriate binary reading files, and not with a text editor. As a rule, these files are smaller than text files. The binary file format does not support headers and is written and read as a single block of data.

In this example the compressed file is not much smaller than the uncompressed file, because the data are random numbers - compression does not help much for this type of data.

```
In [21]:
```arr = np.random.randn(3,4)
print('original data:\n{}'.format(arr))
#write/read uncompressed binary file
filename = 'arrfile.npy'
np.save(filename, arr)
arr2 = np.load(filename)
print('read from binary file {} data:\n{}'.format(filename, arr2))
#write/read binary file
filename = 'arrfile.npz'
np.savez(filename, arr=arr)
arrz = np.load(filename)
print('read from text file {} data:\n{}'.format(filename, arrz['arr']))
#write/read compressed binary file
filename = 'arrfilez.npz'
np.savez_compressed(filename, arr=arr)
arrz = np.load(filename)
print('read from text file {} data:\n{}'.format(filename, arrz['arr']))

```
```

`pyradi.ryplot`

also provides specialist functionality to read files. The use of these files is demonstrated below in the colour example.

The function

`ryfiles.loadColumnTextFile`

was written to load spectral files in pyradi. The first column contains the spectral domain (e.g., wavelength or wavenumber), and the remaining columns contain the spectral variables (e.g.. transmittance). The spectral domain can be scaled; from nanometre to micormetre. The spectral variables can be interpolated to a new spectral domain. This is often required because the spectral domain values used in a calculation is not the same as the spectral domain used in the data file. In short using`pyradi.ryfiles.loadColumnTextFile`

can shorten the code significantly when reading in spectral data files.It is common practice to document the data in a file with a one-line header containing the names of the columns in the file. The function

`ryfiles.loadHeaderTextFile`

reads the column headers from the first line of a text file. The column names in this header file must be comma separated, because some column headings may contain spaces.A function is provided to take an arbitrary string and turn it into a legal filename by removing characters not allowed in filenames. A default set is provided but you may provide your own. See

`ryfiles.cleanFilename`

.Two functions are provided to unzip gzip files

`ryfiles.unzipGZipfile`

and untar tar files`ryfiles.untarTarfile`

. Finally, a function is provided to download a file from the internet, given the URL`ryfiles.downloadFileUrl`

.

(Sec 6.6) describes how to calculate the signal that a given sensor with response $S$ would receive from a source $L$, through some medium with transmittance $\tau_{01\lambda}$. Suppose that this sensor has a spectral filter $\tau_f$ with a narrow (but nonzero) spectral width. This filter is very narrow compared with its central wavelength, say $\Delta \lambda = 0.01 \lambda_c$. Such radiometers are used to determine the spectral radiance of sources or to measure the spectral transmittance of the atmosphere. The apparent irradiance measured by such a system can be written [from (Eq 6.13)] $$ E_{\lambda_c}= k\int_{0}^{\infty} \epsilon_{0\lambda} L_{0\lambda} \,\tau_{a\lambda} \,\tau_{f\lambda} {S}_\lambda \lambda, $$ where $k$ accounts for the geometrical factors such as source area, orientation, and distance. This equation can be written as $$ E_{\lambda_c}= k\int_{\lambda_c - \frac{\Delta\lambda}{2}}^{\lambda_c+\frac{\Delta\lambda}{2}} \epsilon_{0\lambda} L_{0\lambda} \,\tau_{a\lambda} \,\tau_{f\lambda} { S}_\lambda \lambda. $$ By change of variable $\lambda = \lambda_c-x$, $$ E_{\lambda_c}= k\int_{-\frac{\Delta\lambda}{2}}^{+\frac{\Delta\lambda}{2}} \epsilon_{0x} L_{0x} \,\tau_{ax} \,\tau_{f(\lambda_c-x)} { S}_x x. $$ These equations show very clearly that the irradiance measured with the filter centred around wavelength $\lambda_c$ includes source energy from $\lambda_c - \frac{\Delta\lambda}{2}$ to $\lambda_c+\frac{\Delta\lambda}{2}$. Apart from the spectral selection, the filter has an additional effect by smoothing the spectrum being observed because the filter has a nonzero spectral width.

The equation above is called a convolution integral because it describes the convolution between the product $(\epsilon_{0\lambda_c}L_{0\lambda_c} \,\tau_{a\lambda_c} \, { S}_{\lambda_c})$ and $\tau_{f}$. In linear systems terminology, the observed spectral source radiance is being convolved with the filter spectral transmittance. To investigate the effects of this convolution consider the two cases: (1) the observed source has little variation over the filter passband, and (2) the observed source varies significantly over the filter passband.

If the product $(\epsilon_{0\lambda_c} L_{0\lambda_c} \,\tau_{a\lambda_c} \, { S}_{\lambda_c})$ is more or less constant over the filter passband $\Delta\lambda$, the equation can be written to show that the convolution has little effect other than some insignificant amount of smoothing: $$

k\,(\epsilon_{0\lambda*c} L*{0\lambda*c} \,\tau*{a\lambda*c} \, { S}*{\lambda*c})
\int*{-\frac{\Delta\lambda}{2}}^{+\frac{\Delta\lambda}{2}}
\tau_{f(\lambda*c-x)}
x \approx \,
k\,\epsilon*{0\lambda*c} L*{0\lambda*c} \,\tau*{a\lambda*c} \, { S}*{\lambda*c}\, \tau*{f}\,\Delta x.
$$
If the product
$(\epsilon*\lambda L*{\lambda}\;\tau*{a\lambda}\;{ S}*\lambda)$
varies significantly over the filter passband $\Delta\lambda$,
the convolution equation cannot be simplified. In this case, the convolution attenuates and smears out the finer detail in the spectral information. If the actual spectral line is very narrow the measured line will approximate the filter resolution and will be totally erroneous unless the filter spectral smear effect is compensated by deconvolution.

`pyradi.ryutils.convolve`

function is used to do the convolution.

```
In [22]:
```samplingresolution=0.5
wavenum=np.linspace(0, 100, int(100/samplingresolution))
inspectral=np.zeros(wavenum.shape)
inspectral[int(10/samplingresolution)] = 1
inspectral[int(11/samplingresolution)] = 1
inspectral[int(30/samplingresolution)] = 1
inspectral[int(45/samplingresolution)] = 1
inspectral[int(55/samplingresolution)] = 1
inspectral[int(70/samplingresolution)] = 1
inspectral[int(75/samplingresolution)] = 1
inwinwidth=1
outwinwidth=5
outspectral, windowfn = ryutils.convolve(inspectral, samplingresolution, inwinwidth, outwinwidth)
convplot = ryplot.Plotter(1, 1, 1, figsize=(8,4))
convplot.plot(1, wavenum, inspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['r-'],label=['Input'],legendAlpha=0.5)
convplot.plot(1, wavenum, outspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['g-'],label=['Output'],legendAlpha=0.5);

```
```

The Savitzky-Golay filter smooths signals by by fitting successive sub-sets of adjacent data points with a low-degree polynomial by the method of linear least squares.

A Python implementation of this filter is given here. A one-dimensional version of the filter is also included in `pyradi.ryutils.SavitzkyGolay`

.

The example below shows a considerable ringing on the spiky input signal. Clearly this filter is not appropriate for this particular input signal. Another example (atmospheric transmittance) is shown below which better demonstrates the Savitzky-Golay filter performance.

```
In [23]:
```outspectral = ryutils.savitzkyGolay1D(inspectral, window_size=21, order=3, deriv=0, rate=1)
convplot = ryplot.Plotter(1, 1, 1, figsize=(8,4))
convplot.plot(1, wavenum, inspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['r-'],label=['Input'],legendAlpha=0.5)
convplot.plot(1, wavenum, outspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['g-'],label=['Output'],legendAlpha=0.5);

```
```

```
In [24]:
```from scipy.interpolate import interp1d
import pyradi.ryfiles as ryfiles
bunsenfile = ryfiles.downloadFileUrl(url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/bunsenspec.txt')
specRad = np.loadtxt(bunsenfile, comments='%', delimiter=' ')
atmofile = ryfiles.downloadFileUrl(url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/atmotrans5m.txt')
tauAtmo = np.loadtxt(atmofile, comments='%', delimiter=' ' )
wavenum = tauAtmo[:, 0]
tauA = tauAtmo[:, 1]
# convolve transmittance from 1cm-1 to 4 cm-1
tauAtmo4, windowfn = ryutils.convolve(tauA, 1, 1, 4)
# convolve transmittance from 1cm-1 to 40 cm-1
tauAtmo40, windowfn = ryutils.convolve(tauA, 1, 1, 40)
#interpolate bunsen spectrum to atmo sampling
#first construct the interpolating function, using bunsen
bunInterp1 = interp1d(specRad[:,0], specRad[:,1])
#then call the function on atmo intervals
bunsen = bunInterp1(wavenum)
atmoplot = tauA.copy()
atmoplot = np.vstack((atmoplot, tauAtmo4, tauAtmo40))
convplot02 = ryplot.Plotter(1, 1, 2,figsize=(20,5))
convplot02.plot(1, wavenum, atmoplot.T, "Atmospheric Transmittance", r'Wavenumber cm$^{-1}$',
r'Transmittance', ['r-', 'g-','b-'],label=['1 cm-1', '4 cm-1', '40 cm-1' ],legendAlpha=0.5);
convplot02.plot(2, wavenum, atmoplot.T, "Atmospheric Transmittance", r'Wavenumber cm$^{-1}$',
r'Transmittance', ['r-', 'g-','b-'],label=['1 cm-1', '4 cm-1', '40 cm-1' ],legendAlpha=0.5,
pltaxis=[3700, 3900, 0, 1]);

```
```

```
In [25]:
```bunsenPlt = ryplot.Plotter(1,3, 2, figsize=(15,7))
bunsenPlt.plot(1, wavenum, bunsen, "Bunsen Flame Measurement 4 cm-1", r'',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5])
bunsenPlt.plot(2, wavenum, bunsen, "Bunsen Flame Measurement 4 cm-1", r'',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5])
bunsenPlt.plot(3, wavenum, tauA, "Atmospheric Transmittance 1 cm-1", r'',
r'Transmittance', ['r-'])
bunsenPlt.plot(4, wavenum, tauAtmo4, "Atmospheric Transmittance 4 cm-1", r'',
r'Transmittance', ['r-'])
bunsenPlt.plot(5, wavenum, bunsen/tauA, "Atmospheric-corrected Bunsen Flame Measurement 1 cm-1", r'Wavenumber cm$^{-1}$',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5])
bunsenPlt.plot(6, wavenum, bunsen/tauAtmo4, "Atmospheric-corrected Bunsen Flame Measurement 4 cm-1", r'Wavenumber cm$^{-1}$',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5]);

```
```