<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 12a Flame sensor analysis
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.  
<img src="https
Step1:  The system has the elements in the following picture
Step2:  The sensor has an aperture area of 7.8e-3 m$^2$ and a field of view of 1e-4 sr.
The InSb or PbSe detector has a peak responsivity of 2.5 A/W and spectral
response The preamplifier transimpedance is 10000 V/A.
The sensor filter spectral transmittance is  calculated with pyradi.ryutils.sfilter and parameters  $\tau_s = 0.0001$, $\tau_p = 0.9$, $s = 12$, $\Delta\lambda = 0.8$ $\mu$m, $\lambda_c = 4.3 \mu$m.
The flame area is 1 m$^2$, the flame temperature is 1000 C. The emissivity is
0.1 over most of the spectral band, due to carbon particles in the flame. At
4.3 $\mu$m there is a strong emissivity rise due to the hot CO$_2$ in the flame.
The emissivity can be calculated with pyradi.ryutils.sfilter and  parameters  $\tau_s = 0.1$, $\tau_p = 0.7$, $s = 6$, $\Delta\lambda = 0.45 \mu$m, $\lambda_c = 4.33 \mu$m.
The distance between the flame and the sensor is 1000 m. We use the  Modtran
Tropical climatic model. The path is oriented such that the sensor stares out to space,
at a zenith angle of 88.8 deg
The peak in the flame emissivity and the dip in atmospheric transmittance are
both centered around the 4.3 $\mu$m CO$2$ band. In order to determine the flux
transferred one must perform a spectral calculation taking these strong
spectral variations into account.
From (Eq 6.16), the signal caused by the atmospheric path radiance is given by
\begin{equation}
v{\rm path}=
\widehat{{ R}}\, Z_t \,A_{1} \, \Omega\,
\int_{0}^{\infty}
 L_{{\rm path}\lambda}
 { S}{\lambda}
d\lambda,
\end{equation}
where 
$\widehat{{ R}}$ is the peak detector responsivity, 
$Z_t$ is the proximity electronics (preamplifier) gain, 
$A{1}$ the area of the sensor aperture, 
$\Omega$ is the detector solid angle field of view,
$L_{{\rm path}}$ is the path radiance, 
${ S}=\widetilde{{ R}\lambda} \tau{f\lambda}$,
$\widetilde{{ R}\lambda}$ is the normalised detector response, and 
$\tau{f\lambda}$ is the atmospheric transmittance.  Note that the path radiance term does not have an atmospheric transmittance factor because the radiance is the net effect of the atmosphere see (Eq 4.10).
The signal when the flame is present is given by
\begin{equation}
v_{\rm flame}=
\widehat{{ R}}\, Z_t \,A_{1} \left[
\left(\frac{A_{0} }{R_{01}^2}\right)\,
\int_{0}^{\infty}
 \epsilon_{0\lambda} L_{0\lambda}
\,\tau_{a\lambda} { S}\lambda 
d\lambda
+
\left(\Omega -  \frac{A{0} }{R_{01}^2}  \right)\,
\int_{0}^{\infty}
 L_{{\rm path}\lambda}
 { S}{\lambda}
d\lambda
\right],
\end{equation}
where 
$A{0}$ is the flame area because  $\cos\theta_0=1$, and
$R_{01}^2$ is the distance between the flame and the sensor.
Note that the flame only fills a small portion of the field of view and the rest is filled by sky path radiance.
It is convenient to define the sensor field of view fill factor as $f=\Omega_T/\Omega_P$ where $\Omega_T$ is the target solid angle and $\Omega_P$ is the sensor pixel field of view.  If the target solid angle total fills or exceeds the pixel field of view, the fill factor is unity.  If the target solid angle is smaller than the pixel field of view solid angle the fill factor is smaller than unity.
Step3:  Create the spectral range to be used for the calculation. Create the range in wavenumber domain and immediately create a corresponding wavelength domain. These points correspond value by value, so they can be used interchangeably in conversions and plotting.  Use pyradi.ryutils.convertSpectralDomain for this purpose.
Step4:  Read the input data files.  The relevant files are all gzipped up and available for download from the pyradi website.  Use the ryfiles.downloadUntar function in ryfiles to do the download and unzipping to a local directory in one single function call.
Step5:  Load atmospheric transmittance from  a tape7 file created with Modtran.  The transmittance is in column 1, and the path radiance is in column 9.    The Modtran data is in wavenumber spectral domain with
5 cm$^{-1}$ intervals, but we want to work in wavelength with 2.5 cm$^{-1}$ resolution.  The interpolation is performed as part of the file reading process using ryfiles.loadColumnTextFile.  Modtran calculates path radiance in units of W/(cm$^2$.sr.cm$^{-1}$), but it is required in units of W/(m$^2$.sr.$^{-1}$) - the conversion is done while reading the data from the file.
Step6:  Convert path radiance spectral density from 1/cm$^{-1}$ to 1/$\mu$m, at the sample wavenumber points.
Use the function pyradi.ryutils.convertSpectralDensity for this purpose.
Step7:  Load the detector file in wavelengths, and interpolate to  required spectral domain values.
Step8:  Construct the flame emissivity from parameters, using the pyradi.ryutils.sfilter function. Normally one would use measured emissivity data.
Step9:  Construct the filter spectral transmittance from parameters, using the pyradi.ryutils.sfilter function. Normally one would use measured filter data.
Step10:  Plot the data read in and calculated above.  The visual inspection is a critical tool in validating the calculations.
In this particular instance observe how closely the path radiance matches the blackbody radiance at the same temperature. From this simple observation we can conclude that the atmospheric path emissivity is near unity.
Step11:  The signal when the flame is present is given by
\begin{equation}
v_{\rm flame}=
\widehat{{ R}}\, Z_t \,A_{1} \left[
\left(\frac{A_{0} }{R_{01}^2}\right)\,
\int_{0}^{\infty}
 \epsilon_{0\lambda} L_{0\lambda}
\,\tau_{a\lambda} { S}\lambda 
d\lambda
+
\left(\Omega -  \frac{A{0} }{R_{01}^2}  \right)\,
\int_{0}^{\infty}
 L_{{\rm path}\lambda}
 { S}{\lambda}
d\lambda
\right],
\end{equation}
where 
$A{0}$ is the flame area because  $\cos\theta_0=1$, and
$R_{01}^2$ is the distance between the flame and the sensor.
Note that the flame only fills a small portion of the field of view and the rest is filled by sky path radiance.
Define sensor scalar parameters.
Step12:  Define the  flame properties.
Define the fill factor as the portion of the field of view that is filled by the flame, then the portion fill by background radiance is (1-fill).  The fill fraction is calculated as the solid angle filled by the flame divided by the solid angle of the field of view.
Step13:  Calculate the path radiance signal (absence of the flame). First multiply out the spectral components  - this can be done because all variables are on a common spectral sampling. Then sum the integration by using the numpy.trapz function. Finally multiply the integral with the sensor parameters to determine the output signal.
Step14:  Calculate the  signal with the flame in the field of view. In this case there is both path radiance and flame radiance in the field of view; the two must be added together.  
The flame radiance is calculated by using Planck's law and multiplying with the spectral emissivity. The spectral calculation of the flame signal proceeds similarly to the calculation of the spectral background radiance. The flame signal is then determined by integrating over all wavelengths and then multiplying with the sensor parameters.
The signals from the background and flame are then added together, using the fill factor to determine the fraction contributed by the background.
Step15:  Print out the results
Step16:  It is clear that the flame signal is several times larger than the path radiance
signal, even though the flame only fills 1% of the sensor FOV. The
severity by which the atmosphere attenuates the CO$_2$ exitance from the
flame is shown above. Note the path radiance is much smaller than the flame radiance.
Python and module versions, and dates
<ASSISTANT_TASK:>
Python Code:
from IPython.display import display
from IPython.display import Image
from IPython.display import HTML
Explanation: 12a Flame sensor analysis
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.  
<img src="https://zenodo.org/badge/doi/10.5281/zenodo.9910.png"   align="left"/>
The date of this document and module versions used in this document are given at the end of the file.
Feedback is appreciated: neliswillers at gmail dot com.
Overview
This notebook provides a simple worked example of a small sensor observing a flame
on a smokestack. 
The objective is to calculate the signal of a simple infrared sensor, detecting the
presence or absence of a flame in the sensor field of view. The sensor is
pointed to an area just outside a furnace vent, against a clear sky
background. The sensor must detect a change in signal indicating the presence
of a flame at the vent.  
This problem is taken from my book, Chapter 9.
End of explanation
display(Image(filename='images/sensorsystem.png'))
Explanation: The system has the elements in the following picture:
End of explanation
%matplotlib inline
import numpy as np
import pyradi.ryfiles as ryfiles
import pyradi.ryplot as ryplot
import pyradi.ryplanck as ryplanck
import pyradi.ryutils as ryutils
Explanation: The sensor has an aperture area of 7.8e-3 m$^2$ and a field of view of 1e-4 sr.
The InSb or PbSe detector has a peak responsivity of 2.5 A/W and spectral
response The preamplifier transimpedance is 10000 V/A.
The sensor filter spectral transmittance is  calculated with pyradi.ryutils.sfilter and parameters  $\tau_s = 0.0001$, $\tau_p = 0.9$, $s = 12$, $\Delta\lambda = 0.8$ $\mu$m, $\lambda_c = 4.3 \mu$m.
The flame area is 1 m$^2$, the flame temperature is 1000 C. The emissivity is
0.1 over most of the spectral band, due to carbon particles in the flame. At
4.3 $\mu$m there is a strong emissivity rise due to the hot CO$_2$ in the flame.
The emissivity can be calculated with pyradi.ryutils.sfilter and  parameters  $\tau_s = 0.1$, $\tau_p = 0.7$, $s = 6$, $\Delta\lambda = 0.45 \mu$m, $\lambda_c = 4.33 \mu$m.
The distance between the flame and the sensor is 1000 m. We use the  Modtran
Tropical climatic model. The path is oriented such that the sensor stares out to space,
at a zenith angle of 88.8 deg
The peak in the flame emissivity and the dip in atmospheric transmittance are
both centered around the 4.3 $\mu$m CO$2$ band. In order to determine the flux
transferred one must perform a spectral calculation taking these strong
spectral variations into account.
From (Eq 6.16), the signal caused by the atmospheric path radiance is given by
\begin{equation}
v{\rm path}=
\widehat{{ R}}\, Z_t \,A_{1} \, \Omega\,
\int_{0}^{\infty}
 L_{{\rm path}\lambda}
 { S}{\lambda}
d\lambda,
\end{equation}
where 
$\widehat{{ R}}$ is the peak detector responsivity, 
$Z_t$ is the proximity electronics (preamplifier) gain, 
$A{1}$ the area of the sensor aperture, 
$\Omega$ is the detector solid angle field of view,
$L_{{\rm path}}$ is the path radiance, 
${ S}=\widetilde{{ R}\lambda} \tau{f\lambda}$,
$\widetilde{{ R}\lambda}$ is the normalised detector response, and 
$\tau{f\lambda}$ is the atmospheric transmittance.  Note that the path radiance term does not have an atmospheric transmittance factor because the radiance is the net effect of the atmosphere see (Eq 4.10).
The signal when the flame is present is given by
\begin{equation}
v_{\rm flame}=
\widehat{{ R}}\, Z_t \,A_{1} \left[
\left(\frac{A_{0} }{R_{01}^2}\right)\,
\int_{0}^{\infty}
 \epsilon_{0\lambda} L_{0\lambda}
\,\tau_{a\lambda} { S}\lambda 
d\lambda
+
\left(\Omega -  \frac{A{0} }{R_{01}^2}  \right)\,
\int_{0}^{\infty}
 L_{{\rm path}\lambda}
 { S}{\lambda}
d\lambda
\right],
\end{equation}
where 
$A{0}$ is the flame area because  $\cos\theta_0=1$, and
$R_{01}^2$ is the distance between the flame and the sensor.
Note that the flame only fills a small portion of the field of view and the rest is filled by sky path radiance.
It is convenient to define the sensor field of view fill factor as $f=\Omega_T/\Omega_P$ where $\Omega_T$ is the target solid angle and $\Omega_P$ is the sensor pixel field of view.  If the target solid angle total fills or exceeds the pixel field of view, the fill factor is unity.  If the target solid angle is smaller than the pixel field of view solid angle the fill factor is smaller than unity.
End of explanation
waven = np.arange(2000.0,  3300.0,  2.5).reshape(-1, 1)
wavel= ryutils.convertSpectralDomain(waven,  type='nl')
Explanation: Create the spectral range to be used for the calculation. Create the range in wavenumber domain and immediately create a corresponding wavelength domain. These points correspond value by value, so they can be used interchangeably in conversions and plotting.  Use pyradi.ryutils.convertSpectralDomain for this purpose.
End of explanation
tgzFilename = 'flamesensordata.tgz'
destinationDir = '.'
tarFilename = 'flamesensordata.tar'
url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/'
dlNames = ryfiles.downloadUntar(tgzFilename, url, destinationDir, tarFilename)
print('filesAvailable is {}'.format(dlNames))
Explanation: Read the input data files.  The relevant files are all gzipped up and available for download from the pyradi website.  Use the ryfiles.downloadUntar function in ryfiles to do the download and unzipping to a local directory in one single function call.
End of explanation
tauA = ryfiles.loadColumnTextFile('./path1kmflamesensor.txt',
    [1], abscissaOut=waven, comment='%')
lpathwn = ryfiles.loadColumnTextFile('./pathspaceflamesensor.txt',
    [9], abscissaOut=waven, ordinateScale=1.0e4, comment='%')
Explanation: Load atmospheric transmittance from  a tape7 file created with Modtran.  The transmittance is in column 1, and the path radiance is in column 9.    The Modtran data is in wavenumber spectral domain with
5 cm$^{-1}$ intervals, but we want to work in wavelength with 2.5 cm$^{-1}$ resolution.  The interpolation is performed as part of the file reading process using ryfiles.loadColumnTextFile.  Modtran calculates path radiance in units of W/(cm$^2$.sr.cm$^{-1}$), but it is required in units of W/(m$^2$.sr.$^{-1}$) - the conversion is done while reading the data from the file.
End of explanation
(dum, lpathwl) = ryutils.convertSpectralDensity(waven, lpathwn, type='nl')
Explanation: Convert path radiance spectral density from 1/cm$^{-1}$ to 1/$\mu$m, at the sample wavenumber points.
Use the function pyradi.ryutils.convertSpectralDensity for this purpose.
End of explanation
detR = ryfiles.loadColumnTextFile('./detectorflamesensor.txt',
    [1], abscissaOut=wavel, comment='%')
Explanation: Load the detector file in wavelengths, and interpolate to  required spectral domain values.
End of explanation
emis = ryutils.sfilter(wavel,center=4.33, width=0.45, exponent=6, taupass=0.8,
    taustop=0.1 )
Explanation: Construct the flame emissivity from parameters, using the pyradi.ryutils.sfilter function. Normally one would use measured emissivity data.
End of explanation
#construct the sensor filter from parameters
sfilter = ryutils.sfilter(wavel,center=4.3, width=0.8, exponent=12,
           taupass=0.9, taustop=0.0001)
Explanation: Construct the filter spectral transmittance from parameters, using the pyradi.ryutils.sfilter function. Normally one would use measured filter data.
End of explanation
plotdata = detR
plotdata = np.hstack((plotdata,emis))
plotdata = np.hstack((plotdata,sfilter))
plotdata = np.hstack((plotdata,tauA))
label = ['Detector','Emissivity','Filter','Atmosphere transmittance']
plot1= ryplot.Plotter(1, 1, 2,'Flame sensor',figsize=(12,6))
plot1.plot(1, wavel, plotdata, "Spectral","Wavelength [$\mu$m]",
    "Relative magnitude",
    label=label, maxNX=10, maxNY=10)
LbbTropical = ryplanck.planck(wavel, 273+27, type='el').reshape(-1, 1)/np.pi
plotdata = LbbTropical
plotdata = np.hstack((plotdata,lpathwl))
label2=['300 K Planck Law','Tropical path radiance']
plot1.plot(2, wavel, plotdata, 'Comparing blackbody and path radiance', 
           'Wavelength [$\mu$m]','Radiance [W/(m$^2$.sr.$\mu$m)]', label=label2)
Explanation: Plot the data read in and calculated above.  The visual inspection is a critical tool in validating the calculations.
In this particular instance observe how closely the path radiance matches the blackbody radiance at the same temperature. From this simple observation we can conclude that the atmospheric path emissivity is near unity.
End of explanation
opticsArea = 7.8e-3 # optical aperture area [m2]
opticsFOV = 1.0e-4 # sensor field of view [sr]
transZ = 1.0e4 # amplifier transimpedance gain [V/A]
responsivity= 2.5 # detector peak responsivity =A/W]
Explanation: The signal when the flame is present is given by
\begin{equation}
v_{\rm flame}=
\widehat{{ R}}\, Z_t \,A_{1} \left[
\left(\frac{A_{0} }{R_{01}^2}\right)\,
\int_{0}^{\infty}
 \epsilon_{0\lambda} L_{0\lambda}
\,\tau_{a\lambda} { S}\lambda 
d\lambda
+
\left(\Omega -  \frac{A{0} }{R_{01}^2}  \right)\,
\int_{0}^{\infty}
 L_{{\rm path}\lambda}
 { S}{\lambda}
d\lambda
\right],
\end{equation}
where 
$A{0}$ is the flame area because  $\cos\theta_0=1$, and
$R_{01}^2$ is the distance between the flame and the sensor.
Note that the flame only fills a small portion of the field of view and the rest is filled by sky path radiance.
Define sensor scalar parameters.
End of explanation
flameTemperature = 1000.0 + 273.16    # temperature in [K]
flameArea = 1.  # in [m2]
distance = 1000.  # [m]
fill = (flameArea / distance ** 2) /  opticsFOV # how much of FOV is filled
fill = 1 if fill > 1 else fill # limit target solid angle to sensor FOV
Explanation: Define the  flame properties.
Define the fill factor as the portion of the field of view that is filled by the flame, then the portion fill by background radiance is (1-fill).  The fill fraction is calculated as the solid angle filled by the flame divided by the solid angle of the field of view.
End of explanation
inbandirradiancePath =  lpathwn * detR * sfilter * opticsFOV
totalirradiancePath = np.trapz(inbandirradiancePath.reshape(-1, 1),waven, axis=0)[0]
signalPath = totalirradiancePath * transZ * responsivity * opticsArea
Explanation: Calculate the path radiance signal (absence of the flame). First multiply out the spectral components  - this can be done because all variables are on a common spectral sampling. Then sum the integration by using the numpy.trapz function. Finally multiply the integral with the sensor parameters to determine the output signal.
End of explanation
# get spectral radiance in  W/m^2.sr.cm-1
radianceFlame = ryplanck.planck(waven, flameTemperature,  type='en')\
    .reshape(-1, 1)/np.pi
inbandirradianceFlame = radianceFlame * detR * tauA * emis * sfilter *\
    fill * opticsFOV
totalirradianceFlame = np.trapz(inbandirradianceFlame.reshape(-1, 1),
    waven, axis=0)[0]
totalirradiancePathremainder = totalirradiancePath * (1-fill)
signalFlameOnly = totalirradianceFlame *transZ*responsivity *opticsArea
signalFlame = (totalirradianceFlame + totalirradiancePathremainder ) *\
    transZ*responsivity *opticsArea
Explanation: Calculate the  signal with the flame in the field of view. In this case there is both path radiance and flame radiance in the field of view; the two must be added together.  
The flame radiance is calculated by using Planck's law and multiplying with the spectral emissivity. The spectral calculation of the flame signal proceeds similarly to the calculation of the spectral background radiance. The flame signal is then determined by integrating over all wavelengths and then multiplying with the sensor parameters.
The signals from the background and flame are then added together, using the fill factor to determine the fraction contributed by the background.
End of explanation
print('Optics    : area={0} m^2 FOV={1} [sr]'.format(opticsArea, opticsFOV ))
print('Amplifier : gain={0} [V/A]'.format(transZ))
print('Detector  : peak responsivity={0} [A/W]'.format(responsivity))
print('Flame     : temperature={0} [K] area={1} [m^2] distance={2} [m]'.\
    format(flameTemperature, flameArea, distance))
print('          : fill={0} [-]'.format(fill))
print('Flame only: irradiance={0:9.2e} [W/m^2] signal={1:7.4f} [V]'.\
    format(totalirradianceFlame, signalFlameOnly))
print('Path      : irradiance={0:9.2e} [W/m^2] signal={1:7.4f} [V]'.\
    format(totalirradiancePath, signalPath))
print('Flame+Path: irradiance={0:9.2e} [W/m^2] signal={1:7.4f} [V]'.\
    format(totalirradianceFlame + totalirradiancePathremainder , \
        signalFlame))
(dum, iFlamewl) = ryutils.convertSpectralDensity(waven,  inbandirradianceFlame, type='nl')
(dum, iPathwl) = ryutils.convertSpectralDensity(waven,  inbandirradiancePath, type='nl')
plot1= ryplot.Plotter(1, 1, 2,'Flame sensor',figsize=(12,6))
plot1.plot(1, wavel, iFlamewl, "Irradiance","Wavelength [$\mu$m]",
    "Iradiance [W/(m$^2$.$\mu$m)]",plotCol=['r'], label=['Flame'])
plot1.plot(1, wavel, iPathwl, "Irradiance","Wavelength [$\mu$m]",
    "Iradiance [W/(m$^2$.$\mu$m)]",plotCol=['b'], label=['Path'])
plot1.plot(2, waven, inbandirradianceFlame, "Irradiance","Wavenumber [cm$^{-1}$]",
    "Irradiance [W/(m$^2$.cm$^{-1}$)]",plotCol=['r'], label=['Flame'])
plot1.plot(2, waven, inbandirradiancePath, "Irradiance","Wavenumber [cm$^{-1}$]",
    "Irradiance [W/(m$^2$.cm$^{-1}$)]",plotCol=['b'], label=['Path'], maxNX=10, maxNY=10)
Explanation: Print out the results
End of explanation
print(ryutils.VersionInformation('matplotlib,numpy'))
Explanation: It is clear that the flame signal is several times larger than the path radiance
signal, even though the flame only fills 1% of the sensor FOV. The
severity by which the atmosphere attenuates the CO$_2$ exitance from the
flame is shown above. Note the path radiance is much smaller than the flame radiance.
Python and module versions, and dates
End of explanation
<END_TASK>