# to prepare the python environment

import numpy as np
import pandas as pd

from IPython.display import display
from IPython.display import Image
from IPython.display import HTML

import pyradi.rydetector as rydetector
import pyradi.ryplot as ryplot
import pyradi.ryplanck as ryplanck


"""
Explanation: 9 Detector Modelling
End of explanation
"""


# to calculate the semiconductor properties
tempDet = 80  # detector temperature in K
varshni = pd.DataFrame({ 'Madelung' : pd.Series([0.24, 6e-4, 500, tempDet], index=['Eg(0)','A','B','T']),
           'Piprek' : pd.Series([0.235, 3.2e-4, 170, tempDet], index=['Eg(0)','A','B','T']),
           'Dhanaraj' : pd.Series([0.235, 3.1e-4, 452, tempDet], index=['Eg(0)','A','B','T'])})
eg = pd.DataFrame(varshni.apply(lambda x: rydetector.EgVarshni(x['Eg(0)'], x['A'], x['B'], x['T']), axis=0)).T
lambdaC = 1.24 / eg
varshni = pd.concat([varshni, eg], axis=0).rename({0:'Eg'})
varshni = pd.concat([varshni, lambdaC], axis=0).rename({0:'$\lambda_c$'})
Eg = varshni['Madelung']['Eg']
print("Madelung's value for Eg {} is used in the rest of this analysis.".format(Eg))
print("lambda_c is {}.".format(1.24/Eg))
varshni

"""
Explanation: This notebook forms part of a series on computational optical radiometry
These notebooks are constantly revised and updated, please revisit from time to time.  
The source for this notebook is available at 
https://github.com/NelisW/ComputationalRadiometry/blob/master/09a-DetectorModelling.ipynb
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
The pyradi rydetector  module models the primary functionality of optical detectors.  This notebook demonstrates the development of a simple single element InSb detector model suitable for high-level system modelling.  The performance of the detector is described in terms of its figures of merit. The theoretical background to this model is covered the text book.  
Concept
The detector is modelled as a bulk material device, where a fraction of the incident photons are absorbed in the material and then converted to an electrical signal.  
Analysis
Detector material properties
Calculate the bandgap $E_g$ of the InSb material at the detector operating temperature (80 K), using 
(Eq 5.3) imlemented in rydetector.EgVarshni
$$
    E_g = E_{g}(0) -\left(\frac{A T^2}{T + B}\right).
$$ 
The Varshni model (Eq 5.3) is an approximated empirical fit to predict a semiconductor bandgap at a given temperature. The approximation employs three parameters: the reference bandgap at a
given temperature (usually 0 K), plus two parameters to predict the bandgap at other temperatures (relative to the reference bandgap). Several sets of different Varshni parameters were proposed (Table A.6) by Madelung, Piprek, and Dhanaraj, for the references see my book) using different reference bandgap values, and also based on different experimental data. The bandgap calculation considers three different sets of feasible Varshni coefficients.
The detector cutoff wavelength is calculated from $\lambda_c = 1.24/E_g$.
The data storage in the following cell is done in a Pandas dataframe, simply because it is easier to manipulate and display.  The detector parameters are entered as entries into the dataframe and then the Varshni equation is used to calculate Eg in a new row. Next the Eg values are used to calculate the cutoff wavelength $\lambda_c$.  All the data are neatly then displayed in the HTML table.
End of explanation
"""


# to define detector parameters
# wavelength in micrometers, remember to scale down in functions.
wavelength = np.linspace(1, 5.5, 200)    

a0 = 1.9e4 * 100   # absorption coefficient [m-1], Eq3.5 & 3.6 Dereniak
a0p = 800 * 100    # absorption coefficient [m-1] Eq3.5 & 3.6 Dereniak
d1 = 0                # start of the detector depletion layer in meter
d2 = 5e-6             # end of the detector depletion layer in meter
n1 = 1.0              # refraction index of the air
n2 = 3.42        # refraction index of the semiconductor being analyzed
theta1 = 0.0        # radiation incident angle in radians

absorption = rydetector.Absorption(wavelength / 1e6, Eg, tempDet, a0, a0p)
quantumEffic = rydetector.QuantumEfficiency(absorption, d1, d2, theta1, n1, n2)
responsivity = rydetector.Responsivity(wavelength / 1e6,quantumEffic)

# to plot absorption coefficient
absFig = ryplot.Plotter(1,figsize=(9,4))
absFig.semilogY(1,wavelength,absorption,'Absorption Coefficient',\
    r'Wavelength [$\mu$m]',r'Absorption Coefficient [m$^{-1}$]');

# to plot quantum efficiency
QE = ryplot.Plotter(1,figsize=(9,4))
QE.plot(1,wavelength,quantumEffic,'Spectral Quantum Efficiency',\
    r'Wavelength [$\mu$m]','Quantum Efficiency');

# to plot spectral response
Respons = ryplot.Plotter(1,figsize=(9,4))
Respons.plot(1,wavelength,responsivity,'Spectral Responsivity',\
    r'Wavelength [$\mu$m]','Responsivity [A/W]');

"""
Explanation: The spectral absorption coefficient is calculated by using (Eq 5.78) and (Eq 5.79)
$$
\alpha_{\nu} \propto \sqrt{h\nu-E_g} \approx \alpha_{\lambda_c} +  \alpha_0\sqrt{h\nu-E_g} \;\; {\rm for} \;\; h\nu>E_g.
$$
and 
$$
\alpha_{\nu} = \alpha_{\lambda_c}\exp\left[\frac{h\nu-E_g}{kT}\right] \;\; {\rm for} \;\; h\nu<E_g.
$$
the data in (Table A.6). 
The quantum efficiency is calculated by using (Eq 5.87 and 5.88)
$$
\frac{J}{q E_q}=\eta=(1-\rho)(1-e^{- \alpha d})
$$
and
$$
\rho=
\left(
\frac{n_2-n_1}{n_2+n_1}
\right)^2,
$$
with a refractive index of 3.42 for the detector and 1.0 for air, a depletion layer thickness of 5 $\mu$m, and the spectral absorption coefficient calculated above. 
The spectral responsivity is calculated using the spectral quantum efficiency and (Eq. 5.5)
$$
{\cal R}{e\lambda}
=\frac{i_d}{\Phi{e}} 
= \frac{i_d}{\frac{hc}{\lambda}\Phi_{q}}
=\left{
\begin{array}{ll}
\frac{\eta q\lambda G}{hc}
&{\rm for\ \ } \lambda \leq \lambda_c\
0.0
&{\rm for\ \ } \lambda >\lambda_c\
\end{array}
\right.
$$
The gain of the photovoltaic detector is 1.0.
End of explanation
"""


# to prepare for detector signal current calculation
tempSource = 2000       # source temperature in K
emisSource = 1.0          # source emissivity
distance = 0.1      # distance between source and detector
areaSource = 0.000033     # source area in m2
solidAngSource = areaSource / (distance ** 2)
print ("Source temperature             = {0} K".format(tempSource))
print ("Source emissivity              = {0} ".format(emisSource))
print ("Source distance                = {0} m".format(distance))
print ("Source area                    = {0} m2".format(areaSource))
print ("Source solid angle             = {0} sr".format(solidAngSource))
EsourceQL =(emisSource * ryplanck.planck(wavelength,tempSource,'ql') \
    * solidAngSource) / (np.pi )
EsourceEL = (emisSource * ryplanck.planck(wavelength, tempSource,'el')*\
    solidAngSource) / (np.pi )


"""
Explanation: Calculate the spectral photon and radiant irradiance on the detector for a test source blackbody source with unity emissivity and temperature $T_s = 2000$ K.  The source area $A_s = 33$ mm$^2$, located a distance $R = 100$ mm from the detector. The irradiance is calculated as 
$$
E_{s\lambda}=\epsilon_{\lambda} M_\lambda(T_s)\omega/\pi=M(T_s)A_s/(\pi R^2),
$$
where the spectral exitance is given by (Eq. 3.1 and 3.3) for radiant and photon exitance, respectively.
$$
M_{e\lambda}(T)
=\frac{2\pi h c^2}{\lambda^5 \left(e^{hc/(\lambda kT)}-1\right)} 
=\frac{c_{1e\lambda}}{\lambda^5 \left(e^{c_{2\lambda}/(\lambda T)}-1\right)}
$$
and
$$
M_{q\lambda}(T)
=\frac{2\pi c}{\lambda^4 \left(e^{hc/(\lambda kT)}-1\right)}
=\frac{c_{1q\lambda}}{\lambda^4 \left(e^{c_{2\lambda}/(\lambda T)}-1\right)}
$$
End of explanation
"""


# to prepare for detector background current calculation 
tempBkg = 280.0              # background temperature in K
emisBkg = 1.0                # background emissivity
solidAngBkg = np.pi       # background is hemisphere

print ("Background temperature         = {0} K".format(tempBkg))
print ("Background emissivity          = {0} ".format(emisBkg))
print ("Background solid angle         = {0} sr".format(solidAngBkg))

#spectral irradiance for test setup, for both source and background
# in units of photon rate q/(s.m2)
EbkgQL = (emisBkg * ryplanck.planck(wavelength, tempBkg, 'ql') * \
    solidAngBkg ) / (np.pi )
#in radiant units W/m2
EbkgEL = (emisBkg * ryplanck.planck(wavelength, tempBkg, 'el') * \
    solidAngBkg) / (np.pi )


"""
Explanation: Calculate the spectral photon and radiant irradiance on the detector for a hemispherical background blackbody source with unity emissivity and temperature of 280 K.  The calculation is done as above but for an hemispherical solid angle.
End of explanation
"""


# to calculate the detector currents
import scipy.constants as const
areaDet = (200e-6) ** 2   # detector area in m2

iSourceE = np.trapz(EsourceEL * areaDet * responsivity, wavelength)
iBkgE = np.trapz(EbkgEL * areaDet * responsivity, wavelength)
iSourceQ = np.trapz(EsourceQL * areaDet * quantumEffic * const.e,wavelength)
iBkgQ = np.trapz(EbkgQL * areaDet * quantumEffic * const.e, wavelength)
iTotalE = iSourceE + iBkgE
iTotalQ = iSourceQ + iBkgQ

print ("\nDetector current source        = {0:.3e} A and {1:.3e} A".\
    format(iSourceQ,iSourceE))
print ("Detector current background    = {0:.3e} A and {1:.3e} A".\
    format(iBkgQ,iBkgE))
print ("Detector current total         = {0:.3e} A and {1:.3e} A".\
    format(iTotalQ,iTotalE))

"""
Explanation: Calculate the current through the detector for source and background irradiance using a simplified version of (Eq. 6.17):
$$
 i_{\cal R}=
A_{d} \widehat{\cal R}
\int_{0}^{\infty} E_{s\lambda}{\cal R}_\lambda\,
d\lambda,
$$
using the spectral responsivity calculated  above.  The following calculation is performed for both the radiant irradiance and photon irradiance, and should yield the same resulting current.
The detector area is 200 $\mu$m by 200 $\mu$m.
End of explanation
"""


# to calculate the reverse saturation current
eMob = 100.0    # electron mobility - m2/V.s
hMob = 1.0    # hole mobility - m2/V.s
tauE = 1e-8    # electron lifetime - s
tauH = 1e-8    # electron lifetime - s
me = 0.014 * const.m_e    # electron effective mass
mh = 0.43 * const.m_e     # hole effective mass
na = 1e16        # doping = acceptor concentration - m-3
nd = 1e16        # doping = donor concentration - m-3
Isat,ni = rydetector.Isaturation(eMob, tauE, hMob, tauH, me, 
                                 mh, na, nd, Eg, tempDet, areaDet)
print ("I0               = {0} [A] ".format(Isat))
print ("ni               = {0} [m-3] ".format(ni))

"""
Explanation: Calculate the detector reverse saturation current $I_{\rm sat}$. First, calculate the carrier diffusion lengths $L_e$ and $L_h$ by using (Eq. 5.125 and 5.126)
$$
L_h = \sqrt{kT\mu_h \tau_h/q}
$$
and
$$
L_e = \sqrt{kT\mu_e \tau_e/q}
$$
Second, calculate the intrinsic carrier concentration $n_i$ using (Eq. 5.74)
$$
n_i^2 = n p = 4 \left(\frac{2\pi k T}{h^2}\right)^3 \exp\left[\frac{-E_g}{kT}\right]\left(m^_h m^e\right)^{3/2}
$$
The donor concentration is calculated using (Eq. 5.76 and 5.77)
$$
p_n = n^2_i/n_n = n^2_i/n_d.
$$
and 
$$
n_p = n^2_i/p_p = n^2_i/n_a.
$$
to determine the concentration of carriers.
Finally, calculate the reverse-saturation current with (Eq. 5.131)
$$
I{\rm sat}=q A_d \left(
\frac{n_i^2 L_e}{n_a \tau_e} +    \frac{n_i^2 L_h}{n_d \tau_h}
\right)
$$ 
The values used in this calculation for the InSb detector are as follows:
electron mobility 100.0 cm$^2$/(V$\cdot$s);
electron lifetime  $1\times10^{-8}$ s;
hole mobility  1.0  cm$^2$/(V$\cdot$s);
hole lifetime  $1\times10^{-8}$ s;
electron effective mass 0.014 of electron mass;
hole effective mass 0.43 of electron mass;
acceptor concentration $1\times10^{16}$ m$^{-3}$;
donor concentration $1\times10^{16}$ m$^{-3}$;
$E_g$ as calculated  above; and 
detector temperature of 80 K.
End of explanation
"""


# to plot the IV curves
IVbeta = 1.7     # 1 = dominanr diffusion current, 2 = g-r dominant
Vbias = np.linspace(-200e-3, 100e-3, 200) # bias voltage range
ixvDark = rydetector.IXV(Vbias, IVbeta, tempDet, 0, Isat)
ixvBackgnd = rydetector.IXV(Vbias, IVbeta, tempDet, iBkgE, Isat)
ixvTotalE = rydetector.IXV(Vbias, IVbeta, tempDet, iTotalE, Isat)
IXVplot = ryplot.Plotter(1,2,1,figsize=(9,8))
IXVplot.semilogY(1,Vbias,np.abs(ixvDark) * 1e6, plotCol=['g'],label=['Dark'])
IXVplot.semilogY(1,Vbias,np.abs(ixvBackgnd) * 1e6, plotCol=['r'],label=['Background'])
IXVplot.semilogY(1,Vbias,np.abs(ixvTotalE) * 1e6,'IxV Curve',\
    'Voltage [V]','|Current| [uA]',plotCol=['b'],label=['Test source'])
IXVplot.plot(2,Vbias,ixvDark * 1e6, plotCol=['g'],label=['Dark'])
IXVplot.plot(2,Vbias,ixvBackgnd * 1e6, plotCol=['r'],label=['Background'])
IXVplot.plot(2,Vbias,ixvTotalE * 1e6,'IxV Curve',\
    'Voltage [V]','|Current| [uA]',plotCol=['b'],label=['Test source']);

"""
Explanation: Calculate the current-voltage relationship using (Eq 5.128)
$$
I=
I_{\rm sat}
e^{qV/(kT\beta)}
-
I_{\rm sat}
 -I_{ph},
$$
The saturation current $I_{\rm sat}$ is calculated as explained  above.  The photocurrent is determined from (Eq. 5.4 and 5.5)
$$
{\cal R}{q\lambda}=\frac{i_d}{\Phi{q}} 
=\left{\begin{array}{ll}
\eta q G
&{\rm for\ \ } \lambda \leq \lambda_c\
0.0
&{\rm for\ \ } \lambda >\lambda_c\
\end{array},
\right. 
$$
and 
$$
{\cal R}{e\lambda}
=\frac{i_d}{\Phi{e}} 
= \frac{i_d}{\frac{hc}{\lambda}\Phi_{q}}
=\left{
\begin{array}{ll}
\frac{\eta q\lambda G}{hc}
&{\rm for\ \ } \lambda \leq \lambda_c\
0.0
&{\rm for\ \ } \lambda >\lambda_c\
\end{array}
\right.
$$
The value of $\beta$ is determined as explained in (Sec. 5.9.2).  The value for $\beta=1.7$ is used here because InSb operation is dominated by g-r current (for references see my book).
The current-voltage relationship can be calculated for different background flux values, such as 'dark' ($I_{ph}=0$), or for different irradiance values, such as those calculated above.
The following figure shows the current-voltage relationship for three cases: (a) the dark current (i.e., no background flux), (b) hemispherical background flux at 280 K, and (c) the test source flux.  
Note that the y axis is logarithmic to show detail: negative values are shown positive because the logarithm of a negative number is not defined.
Note that the dark current I-V curve is close to the background current I-V curve, while the test source curve is shifted down by almost two orders of magnitude; this is caused by the photon current resulting from the test source flux.
End of explanation
"""


# to calculate R0A product
deltaFreq = 100.0      # measurement or desirable bandwidth - Hertz
iNoisePsd, R0, johnsonPsd, shotPsd = rydetector.Noise(tempDet, IVbeta, Isat, iBkgE)
iNoise = iNoisePsd * np.sqrt(deltaFreq)

print('R0=              = {0} Ohm'.format(R0))
print('NBW=             = {0} Hz'.format(deltaFreq))
print("iNoiseTotal      = {0} A".format(iNoise))
print("iNoise Johnson=  = {0} A".format(johnsonPsd * np.sqrt(deltaFreq)))
print("iNoise shot      = {0} A".format(shotPsd * np.sqrt(deltaFreq)))

"""
Explanation: (Eq. 5.140) is used to calculate the detector's dynamic resistance under zero-bias voltage, 
$$
R_0=-\left.\left(\frac{dI}{dV}\right)^{-1}\right\vert_{V=0} = \frac{\beta k T}{q I_{\rm sat}}
$$
 using the reverse-saturation current calculated  above.
The Johnson noise is calculated using (Eq 5.141) 
$$
i_n^2 
= \frac{4kT}{R_0}
=\frac{4q I_{\rm sat}}{\beta }
=\frac{2(2q I_{\rm sat})}{\beta },
$$
The noise PSD for the three current components: (a) diffusion, (b) thermal excitation, and (c) photocurrent is calculated using (Eq. 5.143).
$$
i^2_n=2q  I_{\rm sat} e^{qV/(kT\beta)}  + 2q  I_{\rm sat}  + 2q   I_{ph}
$$
The total noise is found by adding all four components in quadrature.
End of explanation
"""


# to calculate NEP
NEPower, detectivity = rydetector.NEP(iNoise, responsivity)
NEPfig = ryplot.Plotter(1,figsize=(9,4))
NEPfig.plot(1,wavelength,NEPower,'Spectral Noise Equivalent Power',\
    r'Wavelength [$\mu$m]','NEP [W]',\
    pltaxis=[np.min(wavelength), np.max(wavelength), 0,NEPower[0]]);

"""
Explanation: The NEP is determined from  (Eq. 5.25)
$$
{NEP}\lambda=\frac{i_d \sqrt{\Delta{}f}}{{\cal R}\lambda}
$$
with the spectral responsivity calculated and the noise current calculated above.  The NEP is calculated and shown next.
End of explanation
"""


# to calculate Dstar
#use either of the following two forms:
#dStar = rydetector.DStar(areaDet, deltaFreq, iNoise, responsivity) # units in cm
dStar = rydetector.DStar(areaDet, 1, iNoisePsd, responsivity) # units in cm

DStarfig = ryplot.Plotter(1,figsize=(9,4))
DStarfig.plot(1,wavelength,dStar,'Spectral normalized detectivity',\
    r'Wavelength [$\mu$m]',r'D* [cm$\sqrt{{\rm Hz}}$/W]');

"""
Explanation: The $D^\ast$ is calculated from the area of the detector, the noise bandwidth, the noise, and the responsivity, by (Eq. 5.26)
$$
D^\ast_\lambda
=\frac{\sqrt{\Delta{}f A_d}}{{NEP}\lambda}
=\frac{{\cal R}\lambda \sqrt{\Delta{}f A_d} }{i_d \sqrt{\Delta{}f}  }
$$
End of explanation
"""


# to plot software versions
# https://github.com/rasbt/watermark
try:
    import pyradi.ryutils as ryutils
    print(ryutils.VersionInformation('matplotlib,numpy,pyradi,scipy,pandas'))
except:
    print("pyradi.ryutils not found")

"""
Explanation: Python and module versions, and dates
End of explanation
"""


