7 Optical Sources

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.

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 library has a reasonably complete collection of Planck radiator models, both spectral and wide band. A comprehensive collection of physical constants, pertinent to optical radiation is also included.

The first part of this notebook introduces these functions in the pyradi.ryplanck library.

The second part provides a simple model for scene radiance and the electron counts that can be expected in a detector charge well for the scene flux. These calculations are made with the pyradi.rypflux library.

The third part of the notebook provides some programming details and tricks on how the Planck functions are implemented.

We use the material developed in Electro-Optical System Analysis and Design: A Radiometry Perspective, CJ Willers, SPIE, DOI: 10.1117/3.1001964, 2013 (\url{http://spie.org/Publications/Book/2021423}).

Watch out for integer numbers

It seems that Numpy interprets 400 as an integer int32. When this int32 value is raised to the fourth power, the value overflows, resulting in a negative number. This error can easily happen when entering integer temperature values, raised to a high power.

In ryplanck, all temperature input values are always converted to float with .astype(float) (but only in ryplanck).

Physical constants

This module uses the CODATA physical constants to derive the constants required for radiation calculations. The constants can be printed as a full set, or these can be accessed individually. The CODATA constants are part of the SciPy library scipy.constants, a comprehensive set of easy-to-access constants. Where the values are available in scipy.constants these are used, if not the values are calculated in the ryplanck module.

For more details on CODATA see http://physics.nist.gov/cuu/pdf/RevModPhysCODATA2010.pdf.

Planck radiation law functions

(Sec. 3.1) The equations describing the Planck radiation law are available with spectral variables for wavelength, wavenumber or frequency. The exitance calculation is available in radiant units (watts) and photon rate units (quanta per second). The Planck-law temperature derivatives are also available.

All the functions assume wavelength in micrometres ($\mu$m), wavenumber in cm$^{-1}$ and frequency in hertz. All return values are in area units of m$^2$. Values returned are expressed in exitance units [e.g., W/(m$^2$)], radiance units [e.g., W/(m$^2$.sr)] can be obtained by dividing the output by numpy.pi.

The Planck law functions have the form: ryplanck.planck or planckxx where two-letter code xx denotes the type of spectral variable and return type (see below).

The ryplanck.planck(spectral, temperature, type='el') function docstring states as follows:

Calculates the Planck law spectral exitance from a surface at the stated 
temperature. Temperature can be a scalar, a list or an array. Exitance can 
be given in radiant or photon rate units, depending on user input in type.

Args:
    | spectral (scalar, np.array (N,) or (N,1)):  spectral vector.
    | temperature (scalar, list[M], np.array (M,), (M,1) or (1,M)):  
                Temperature in [K]
    | type (string):
    |  'e' signifies Radiant values in [W/m^2.*].
    |  'q' signifies photon rate values  [quanta/(s.m^2.*)].
    |  'l' signifies wavelength spectral vector  [micrometer].
    |  'n' signifies wavenumber spectral vector [cm-1].
    |  'f' signifies frequency spectral vecor [Hz].

Returns:
    | (scalar, np.array[N,M]):  spectral radiant exitance (not radiance) 
              in units selected.
    | For type = 'el' units will be [W/(m^2.um)].
    | For type = 'qf' units will be [q/(s.m^2.Hz)].
    | Other return types are similarly defined as above.
    | Returns None on error.

Raises:
    | No exception is raised, returns None on error.

The first function parameter is a spectral vector, typically with shape (N,) or (N,1) but strictly a colum vector [not a row vector shape (1,N)]. This spectral domain type must be indicated by the appropriate character in the third parameter.

The second parameter is the object temperature in kelvin. The temperature values can be a scalar, a list, or a column vector with shape (M,), (M,1), or (1,M), but not of shape shape (Q,M).

The third parameter is a string denoting the calculation required. The first character in the string must be one of 'e' for radiant units or 'q' for photon rate units. The second character in the string must be one of 'l' for wavelength, 'n' for wavenumber, or 'f' for frequency.

The return value will be in the format as requested in the third parameter:

• 'el' will return exitance in units of W/(m$^2$.$\mu$m), calling ryplanck.planckel
• 'ql' will return exitance in units of q/(s.m$^2$.$\mu$m), calling ryplanck.planckql
• 'en' will return exitance in units of W/(m$^2$.cm$^{-1}$), calling ryplanck.plancken
• 'qn' will return exitance in units of q/(s.m$^2$ cm$^{-1}$), calling ryplanck.planckqn
• 'ef' will return exitance in units of W/(m$^2$.Hz), calling ryplanck.planckef
• 'qf' will return exitance in units of q/(s.m$^2$.Hz), calling ryplanck.planckqf

The ryplanck.planckxx(spectral, temperature) forms have only two parameters: spectral and temperature.

In addition to to the series of planckxx funtions there are a similar series of dplnckxx series and ryplanck.dplanck to calculate the temperature derivative of the Planck law. These functions have the same function parameters as for planckxx. The return values are similar, except that these have an additional unit (1/K).

The ryplanck.planck and ryplanck.dplanck functions can take different data types for the spectral and temperature parameters.

Either or both spectral and temperature can be scalar, in which case the return type is a scalar.

If spectral is a numpy array, it must have shape (N,) or (N,1). Row vectors of shape (1,N) is not allowed

Temperature can be a scalar, a list, or any form of vector with shape (M,), (M,1), or (1,M), but not of shape (Q,M).

The return value of the ryplanck.planck and ryplanck.dplanck functions always have the spectral domain along the row direction (axis=0) and temperature domain along the column direction (axis=1).

The example above also shows how the ryplanck.planck and ryplanck.dplanck functions can be invoked with multiple temperature values, allowing compact code and fast vectorised calculations.

Stefan-Boltzman wideband exitance

The Stefan-Boltzman equation is available as ryplanck.stefanboltzman

Calculates the total Planck law exitance, integrated over all wavelengths,
from a surface at the stated temperature. Exitance can be given in radiant or
photon rate units, depending on user input in type.

Args:
    | temperature (float):  temperature scalar in [K].
    | type (string):  'e' for radiant or 'q' for photon rate exitance.

Returns:
    | (float): integrated radiant exitance in  [W/m^2] or [q/(s.m^2)].
    | Returns a -1 if the type is not 'e' or 'q'

Raises:
    | No exception is raised.

Wien law

The Wien law is so simple that only the constants are provided. The constants are as follows:

ryplanck.pconst.wel  um.K     radiant and wavelength
ryplanck.pconst.wql  um.K     photon rate and wavelength
ryplanck.pconst.wen  cm-1/K   radiant and wavenumber
ryplanck.pconst.wqn  cm-1/K   photon rate and wavenumber
ryplanck.pconst.wef  Hz/K     radiant and frequency
ryplanck.pconst.wqf  Hz/K     photon rate and frequency

Planck law graphs

Example applications

The number of bits required to express colour ratio

Calculate the number of bits required in an image to express the colour ratio between an MTV flare and a relatively cold aircraft fuselage. It is assumed here that the flare and fuselage are resolved, i.e., the individual image pixels are completely filled by the flare and fuselage respectively. First calculate the radiance ratio of aircraft fuselage to MTV flare in 3-5 um band. The flare is modelled as a blackbody with temperature 2200 K and emissivity of 0.15. The aircraft has a temperature of 250 K and an emissivity of 1. The number of bits required is given by $\log_2(L_{\rm flare}/L_{\rm aircraft})$

Cumulative blackbody curves

The following code produces (Fig 3.4) in the book. The top figure shows the spectral radiance of a few real-world objects. The bottom graph shows the normalised cumulative radiance of the same objects. The bottom graph is useful to estimate the percentage of total radiance above/below a wavelength. From this curve it is evident that a tungsten lamp radiates less than 10% of its total flux in the visual spectral band. From this graph it can be concluded that approximately 40% of the sun's output is in the visual spectral band.

Reflected sunlight radiance

(Sec 3.7) derives an equation for the reflected sunlight radiance of a surface: $$ L\lambda= \frac{ \epsilon{\rm s}L{{\rm bb}\lambda}(T{ \rm s})

A_{\rm s}\tau_s\rho_d\cos\thetai}{\pi R{\rm sun}^2}

\psi \epsilon{\rm s}L{{\rm bb}\lambda}(T_{ \rm s})\tau_s \rho_d\cos\thetai, $$ where where $\epsilon{\rm s}\approx 1$ is the emissivity of the sun's surface, $T_s$ is the sun's surface temperature, $\tau_s$ is the atmospheric transmittance between the sun and the surface, and $\thetai$ is the angle between the surface normal and the sun vector. The reflected sun radiance from a perfectly Lambertian surface is then given by $L\lambda=\rhod E{\lambda \rm sun}/\pi$, where $\rhod$ is the surface diffuse reflectance function. The constant $\psi =A{\rm sun}/(\pi R_{\rm sun}^2) =2.1757\times10^{-5}$ [sr/sr] follows from the geometry.

This really simple model does not account for sky radiance or light incident on the surface from other sources, applicable only to situations where the surface is directly illuminated by the sun, and with a object surface normal angle within 70 degrees or so from the sun direction. At larger angles the sun illumination decreases rapidly and other sources such terrain or sky radiance contribute increasingly more.

This section compares the simple model presented above with the results obtained using a Modtran-calculated solar irradiance.

First to the imports and general preparation.

Modtran5 was used to calculate the solar irradiance and transmittance to space for a slant path. Modtran was used in direct solar irradiance mode, with Tropical profile with 23 km Rural visibility. The solar irradiance and transmittance to space is calculated at sea level, for a 45 degree zenith angle, on day 1. The data files are available on the pyradi website. Download, extract and plot the data.

The Modtran tape7 file contents is as follows (the first few lines):

T F 1    3    3    0    1    1    1    1    1    1    1    0    0   0.000   0.00
    1    1    1    0    0    0  23.00000   0.00000   0.00000   0.00000   0.00000
 -99.000 -99.000 -99.000
 -99.00000 -99.00000 -99.00000  4.114728  0.000594 ! H2O & O3 COLUMNS [GM/CM2]
   36TROPICAL MODEL      
   0.00000 100.00000  45.00000 140.37962   0.87815   0.00000    0        0.00000
  -99  -99    1    0
 -99.00000 -99.00000 -99.00000 -99.00000 -99.00000 -99.00000   0.00000 -99.00000
   12495.0   33340.0       1.0       1.0NW        W1    
    0   0    0.000  0    0.000  0    1.000
    FREQ   TRANS     SOL TR  SOLAR
12495.00  0.6706 5.03E-06 7.50E-06                                                              3.995E-01
12496.00  0.4790 3.47E-06 7.23E-06                                                              7.360E-01
12497.00  0.7075 5.07E-06 7.16E-06                                                              3.460E-01
12498.00  0.6852 4.90E-06 7.15E-06                                                              3.780E-01
12499.00  0.5352 3.84E-06 7.17E-06                                                              6.252E-01
...

Extract the data from the file using rymodtran.loadtape7. Then convert the spectral irradiance from W/(cm$^2$.cm$^{-1}$) to W/(m$^2$.$\mu$m) and plot the result.

Consider an object with reflectance of one, with normal vector pointed towards the sun.

It appears that the simple model fits the Modtran data reasonably closely. The shape of the two radiance curves are somewhat different. Considering the simplicity of the model the results are quite acceptable.

(Sec 2.10.3) Next calculate the luminance of the object's surface. The photopic luminance of a source is defined as $$ L_\nu=\int_0^\infty K_{\rm max} V_\lambda L_{e\lambda} d \lambda, $$ where $K_\lambda=K_{\rm max}V_\lambda$ is the spectral photopic efficacy, $V_\lambda$ is the photopic efficiency, $K_{\rm max}=683$ lm/W is the maximum value of photopic efficacy, referenced to a 2042-K blackbody standard source, and $L_{e\lambda}$ is the source's radiance. The spectral shape of photopic luminous efficiency is can be roughly approximated by $$ V_\lambda=1.019\exp(-285.51(\lambda - 0.5591)^2). $$ There are data tables from the CIE standards that describe this shape more accurately but this approximation is sufficient for the present requirement.

According to Wikipedia the sun has a (apparent?) luminance of 1.6e9 cd/m2 at noon, which, if converted with $\psi$=2.1757e-5 gives a value of 34811 lm/(m$^2$.sr). Comparing this to the value of 23440 lm/(m$^2$.sr) calculated above, yields a ratio of 1.48. Considering the fact that the Wikipedia atmosphere is not known and that the simple calculation is for a slant path of 45 degrees, the result from the simple model seems feasible.

Scene radiance

Low light scene radiance

This section develops a simple model for the reflected photon rate flux for low light scenes (thermal emission is not included). The pyradi.rypflux library provides the simple models for scene illumination used here.

Source spectral content

The target object will either radiate by its own thermal self-emittance, reflected incident light, or both. In the visible and near infrared bands the signal is mostly caused by reflection. At wavelengths beyond 6 $\mu$m the signal is mostly caused by thermal self-emission. In the 3--5 $\mu$m spectral band the signal is caused by both sources.

We will be using the colour temperature as an approximation of source spectral content, in the visible and near infrared parts of the spectrum. We will however use the sun's colour temperature for all spectral bands, because the sun's colour temperature is also approximately equal to the surface temperature.

http://sizes.com/units/color_temperature.htm
https://www.firstrays.com/bb1-1286570267-680/start/color_temp.htm
http://www.clarkvision.com/articles/night.photography.image.processing/

We must be careful when using colour temperature, it does not mean that the object has that temperature. It does not even mean that the spectral content of the object can be described by a scaled spectrum of a radiator at that temperature. The colour temperature's only significance is that it very broadly describes the visual colour of an object. From https://en.wikipedia.org/wiki/Color_temperature ''The color temperature of a light source is the temperature of an ideal blackbody radiator that radiates light of comparable hue to that of the light source. Color temperature is a characteristic of visible light that has important applications in lighting, photography, videography, publishing, manufacturing, astrophysics, horticulture, and other fields. In practice, color temperature is only meaningful for light sources that do in fact correspond somewhat closely to the radiation of some black body, i.e., those on a line from reddish/orange via yellow and more or less white to blueish white; it does not make sense to speak of the color temperature of, e.g., a green or a purple light. Color temperature is conventionally stated in the unit of absolute temperature, the Kelvin, having the unit symbol K.''

Source temperature is undefined if the source spectrum deviates significantly from the Planck law spectral shape. For example, for sources with little continuum and high spectral line variation (e.g., sodium and mercury lamps with small fluorescence) colour temperature has no meaning.

We will be using the colour temperature as an approximation of source spectral content, in the visible and near infrared parts of the spectrum. We will however use the sun's colour temperature for all spectral bands, because the sun's colour temperature is also approximately equal to the surface temperature.

The moon's colour temperature:

http://www.cast-lighting.com/landscape-articles/moonlighting-landscape-lighting-design-imitates-nature ''It turns out that the moon absorbs twice as much violet as it does red. For this reason, the moon is slightly reddish. A full moon when it appears directly overhead has a color temperature of about 4,150 K. Compared to incandescent lighting (around 3,000 K), moonlight is slightly blue but not nearly as blue as a bright sunny day (as much as 10 000 K)''

''Johannes Purkinje, a 19th century physiologist, found that at very low light levels, the human eye could no longer perceive the color red, but could still perceive blues and greens. This occurs because the eye's retinal cones (responsible for color perception) require a lot of light. At lower brightness levels, only the retinal rods are activated. These rods (responsible for seeing fine detail and contrast) can only respond to blues and greens.''

''The luminance level at which this perceptual shift occurs is at about .001 candelas/meters. Moonlight has about this same luminance level. This puts moonlight right at the threshold of the Purkinje Shift and this is why the moon appears slightly blue. It should be noted though, that the effect is so slight that the brain easily shifts its perception to judging the moon to be white. If we shift our gaze from the moon to the ground and other objects that it illuminates, the most noticeable effect is not only a blueish color but also an absence of all other colors.''

The rypflux module captures low light level lux values in the rypflux.PFlux.lllux dictionary. The values were originally published in the RCA/Burle Electro-Optics Handbook. Copies of this book floats around the internet.

Relating photometric units with photon rate units

The illuminance table above shows the approximate illuminance levels for different scenarios. The sensor models below require irradiance levels in units of q/(s.m$^s$), not lx/m$^2$. A conversion between the two unit systems is required. This conversion requires two information elements: the illuminance magnitude and spectral content. The magnitude is given in the table, and in the absence of any better information, we would assume the spectral content to be a scaled version of a thermal source at the colour temperature. This use of colour temperature is extending beyond its original purpose and introduces an element of risk. Yet the answer with risk is better than no answer at all.

If the luminance exceeds about $3$ lm/(m$^2$sr) [twilight], the eye is light-adapted, and the cones in the retina are operating. Under light-adapted conditions, the eye's spectral response is called photopic. If the luminance is less than $3\times 10^{−5}$ lm/(m$^2$sr) [overcast night], the eye is dark-adapted. The cones are no longer sensitive, and the rods sense the light.Under dark-adapted conditions, the eye's spectral response is called scotopic.

Under scotopic vision, the eye is not sensitive to color and has no foveal vision. Between photopic and scotopic vision, the eye's response is called mesoptic. Mesopic vision has a spectral response somewhere between photopic and scotopic vision. In the table above the FracPhotop column shows my very rough mixing ratio of photopic and scotopix spectral curves to model mesopic vision. These values are not validated and are purely guesswork.

The eye's photopic spectral response can be roughly approximated by: \begin{equation} V_\lambda=1.019\exp(-285.51(\lambda - 0.5591)^2). \end{equation}

The eye's scotopic spectral response can be roughly approximated by: \begin{equation} V_\lambda^\prime=0.99234\exp(-321.1(\lambda - 0.5028)^2). \end{equation}

There are also data tables from the CIE standards that describe this shape more accurately. The ryutils module has a function luminousEfficiency that returns luminous efficiency data for four different CIE curves as well as the equations shown above.

Extrapolating to other spectral bands

The photopic luminance of a source is defined as \begin{equation} L_\nu=\int_0^\infty K_{\rm max} V_\lambda L_{e\lambda}(T) d \lambda, \end{equation} where $K_{\rm max}V_\lambda$ is the spectral photopic efficacy, $V_\lambda$ is the photopic efficiency, $K_{\rm max}=683$ lm/W is the maximum value of photopic efficacy, referenced to a 2042-K blackbody standard source, and $L_{e\lambda}$ is the source's radiance.

The scotopic luminance of a source is defined as \begin{equation} L_\nu=\int_0^\infty K_{\rm max}^\prime V_\lambda^\prime L_{e\lambda}(T) d \lambda, \end{equation} where $K_{\rm max}^\prime V_\lambda^\prime$ is the spectral scotopic efficacy, $V_\lambda^\prime$ is the scotopic efficiency, $K_{\rm max}^\prime=1700$ lm/W is the maximum value of scotopic efficacy, referenced to a 2042-K blackbody standard source, and $L_{e\lambda}$ is the source's radiance.

The conversion between photometric and photon rate units is based on the calculation of a scaling factor $k$ given by solving the following equation. Under strong light conditions the photopic equation should be used, whereas under very dark conditions the scotopic equation should be used. In the mesopic range a somewhat arbitrary mix of the two is used.

\begin{equation} L_{\nu\textrm{table}} = k L_\nu=k\, \left[ f_P K_{\rm max}\int_0^\infty V_\lambda L_{e\lambda}(T) d \lambda + (1-f_P)K_{\rm max}^\prime\int_0^\infty V_\lambda^\prime L_{e\lambda}(T) d \lambda \right] \end{equation}

where $k$ is the fraction of the expected illuminance given the colour temperature, and $f_P$ is the fraction of photopic flux (scotopic flux will be fraction $1-f_P$.
Calculate $k$ for each entry in the table, using the stated luminance $L_\textrm{table}$ and colour temperature $T$.

Subsequently the wideband integrated photon radiance corresponding to the scenarios in the table can be determined from

\begin{equation} L_{q} = k\int_0^\infty L_{q\lambda}(T) d \lambda \end{equation}

This procedure critically depends on the sources' spectral radiance in the various different spectral bands. For this calculation the approach is taken that for natural scenes the spectral shape can be modelled by a Planck curve at the appropriate colour temperature. City light has a very different spectral character and this derivation does not apply to such lighting conditions.

The closer the spectral bands are to the visible spectral band, the more accurate the results would be. It is questionable whether the colour temperature in the infrared bands are relevant or correct, but there are no better available data.

Quick validation

Of course these simple checks do not validate all the data points in the table, but they do validate the calculation.

Starlight irradiance

The report ''Night sky radiometric measurements during follow-on-evaluation testing of AN/PVS-7 (A,B) at Fort Benning, Ga'' R.J. Stefanik, Center for Night Vision and Electro-Optics, AMSEL-NV-TR-0079, AD-A211273, May 1989, published the following spectral starlight irradiance data

As validation check, simply consider the values at 0.5 $\mu$m and 0.8 $\mu$m and convert to SI units and photon rate. Energy per photon $h\nu=hc/\lambda$, in Joule. $\Phi_q = \Phi_e/(h\nu) = \Phi_e \lambda/(hc)$

1. At 0.5 $\mu$m starlight has an irradiance of 0.8 pW/(cm$^2$.nm) = 0.8e-05 W/(m$^2$.$\mu$m). This corresponds to a photon rate radiance of 1.67e+12 q/(s.m$^2$.sr) over the spectral width in the subsequent calculation. This corresponds closely with the values calculated in the model (1.77e+12 q/(s.m$^2$.sr)).

2. At 0.8 $\mu$m starlight has an irradiance of 2 pW/(cm$^2$.nm) = 2e-05 W/(m$^2$.$\mu$m). This corresponds to a photon rate radiance of 5.144402e+12 q/(s.m$^2$.sr) over the spectral width in the subsequent calculation. This corresponds to a factor of 3.2 of the values calculated in the model (1.61e+12 q/(s.m$^2$.sr)).

Value of $k$ for direct sunlight

The value for k for direct sunlight is calculated in [W, Section 3.7], and found to be 2.17e-5 (when accounting for units differences), which corresponds well with the value of 1.96e-05 calculated here.

Example application

Consider the photon flux for an example scenario (details are given in the code). This example calculates the electron count in the detector charge well, as well as the associated photon noise for the signal.

See the notebook 09b-StaringArrayDetectors.ipynb for theory. The functions to calculate the electrons in the charge well are contained in rystare.py.

Planck law technical implementation detail

The implementation of the Planck functions in Python/Numpy is explained in more detail here. This information is not required to use the functions, but may be useful to Python code developers.

The Planck equations have a form similar to this: $$ 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)} $$

One of the computational issues is that for combinations of short wavelength (high frequency) and low temperature, the value in the exponent becomes very big and the calculation produces NaN or Inf results (values should be near zero or zero). This follows from the range of values that the IEEE floating point format can handle (hence all languages have this problem). This fault condition is tested for and corrected if it does occur (the function correctly returns zeros for those values):

#np.exp() has upper limit in IEEE double range, catch this in Planck calcs
explimit = 709.7

#test value of exponent to prevent infinity, force to exponent to zero
#this happens for low temperatures and short wavelengths
exP =  pconst.c2l / (spectral * temperature)
exP2 = numpy.where(exP<explimit, exP, 1);
p = (pconst.c1el / ( numpy.exp(exP2)-1)) / (spectral ** 5)
#if exponent is exP>=explimit, force Planck to zero
planckA = numpy.where(exP<explimit, p, 0);

In pyradi we try to use vectorised calculations where possible - this means not to use loops where it can be avoided. The objective with this approach is that it is more compact, easier to read and less prone to errors in coding (especially when doing copy-paste operations). The Planck function interface supports a vectorised operations (as detailed above).

The Planck equation does not lend itself to vectorised calculation because of its relative complex form. Attempting to write the Planck equation in vector/matrix form will be quite a challenge. In ryplanck, the spectral variable and temperature variable are accepted in array form, but then 'flattened' into long vectors of matching values of spectral sample and temperature sample - for which the Planck equation can be readily calculated. After the Planck calculation is complete, the flat array is reshaped back, according to the shape of the input arrays.

In a non-vector calculation the Planck equation will require two loops: one for spectral variation and one for the temperature variation. The flatten operation essentially creates combinations of all the input values: each spectral value is paired with each temperature value. The flattened array represents the pairs of values that would be formed by the two nested loops, but instead of looping over time, the flattened array spreads across space (trading looping time for memory). The flattening process is performed such that the spectral/temperature variable structure remains such that it can be reshaped back into the original format. Note that flattening and reshaping has almost no computing overhead, it does not make or change variables, it only changes the view on the data.

There are twelve variations of the Planck equation and its temperature derivative. It was an arduous task to include exactly the same flattening and reshaping code in all twelve functions - and then changing them all when a change was required. Instead there is now only one implementation of the flattening and reshaping code, and it applies to all twelve functions. The code is implemented in a Python decorator (not unique to Python, many other languages also have decorators). The decorator is almost like a wrapping function that is applied to the user function (decorating the original). Using the decorator function it is possible to execute code before the user function is called and then also after the user function is called. In this application, the flattening is done before the function call and the reshaping is done after the function call - keeping only one copy of the code, but applying it to all the functions. Python decorators are also used for many other applications.

The code is as follows. First the user function showing the decorator, indicated here as @fixDimensions. Note that the user function (docstring removed) focuses only on the Planck equation and not on array dimensions. The fixDimensions intercepts all the input and output from this function and process according to the decorator requirements.

@fixDimensions
def planckel(spectral, temperature):

    #test value of exponent to prevent infinity, force to exponent to zero
    #this happens for low temperatures and short wavelengths
    exP =  pconst.c2l / (spectral * temperature)
    exP2 = numpy.where(exP<explimit, exP, 1);
    p = (pconst.c1el / ( numpy.exp(exP2)-1)) / (spectral ** 5)
    #if exponent is exP>=explimit, force Planck to zero
    planckA = numpy.where(exP<explimit, p, 0);

    return planckA

The decorator function (docstring removed) focuses only on manipulation of input and output data formats, and has no idea what is happening inside the Planck functions, which it calls.

The decorator function:

• checks the shapes of the input variables and takes error action if necessary.
• converts all input to numpy arrays, interpreted as float values (if not already so)
• reshapes and creates a meshgrid to provide all the combinations of input sample values.
• flattens the arrays by using the ravel function.

def fixDimensions(planckFun):
  @wraps(planckFun)
  def inner(spectral, temperature):

    #confirm that only vector is used, break with warning if so.
    if isinstance(temperature, numpy.ndarray):
        if len(temperature.flat) != max(temperature.shape):
            print('ryplanck: temperature must be of shape (M,), (M,1) or (1,M)')
            return None
    #confirm that no row vector is used, break with warning if so.
    if isinstance(spectral, numpy.ndarray):
        if len(spectral.flat) != spectral.shape[0]:
            print('ryplanck: spectral must be of shape (N,) or (N,1)')
            return None
    tempIn = numpy.array(temperature, copy=True,  ndmin=1).astype(float)
    specIn = numpy.array(spectral, copy=True,  ndmin=1).astype(float)
    tempIn = tempIn.reshape(-1,1)
    specIn = specIn.reshape(-1,1)

    #create flattened version of the input dataset
    specgrid, tempgrid = numpy.meshgrid(specIn,tempIn)
    spec = numpy.ravel(specgrid)
    temp = numpy.ravel(tempgrid)

At this point the Planck function is called with the flattened arrays.

    #this is the actual planck calculation
    planckA = planckFun(spec,temp) 

After the function returns, the return values are handed back to the decorator function which then proceeds to reshape the array back into a format the user can work with.

    #now unflatten to proper structure again, spectral along axis=0
    if tempIn.shape[0] == 1 and specIn.shape[0] == 1:
        rtnVal = planckA[0]
    elif tempIn.shape[0] == 1 and specIn.shape[0] != 1:
        rtnVal = planckA.reshape(specIn.shape[0],)
    elif tempIn.shape[0] != 1 and specIn.shape[0] == 1:
        rtnVal = planckA.reshape(specIn.shape[0],-1)
    else:
        rtnVal = planckA.reshape(tempIn.shape[0],-1).T

    return rtnVal
  return inner

Python and module versions, and dates