Thin-film approximate point scattered function

We first review the basic theory of the point scattered function, which is based on four basic assumptions:

1. The scattered neutrons are homogeneous.

2. The neutrons are scattered only once in the sample.

3. The source and scattered neutrons have the same attenuation coefficient in the sample.

4. The sample width is large enough to ensure that the scattered neutrons penetrate the bottom of sample (not the sides) and reach the detector.

These are reasonable assumptions. Because neutrons are electrically neutral, each should be homogeneous after scattering. The vast majority of neutrons also scatter only once in the sample and therefore the need for multiple scattering corrections is not significant.

Figure 1 shows a schematic of neutron radiography and the angular response of the detector.

Fig. 1Full-screen View

Schematic diagram of neutron radiography and angular response of the detector

In Fig. 1, r is the coordinate parameter of the detector, d is the thickness of the detector converter, θ is the scattering angle, D is the distance from the sample to the detector, T is the thickness of the sample, t is the location where scattering occurs in the sample, and R and l are both geometric parameters. The following relationship exists between the geometric parameters. $\cos \theta =\frac{D+t}{R}\mathrm{,}R=\sqrt{{\left(t+D\right)}^2+r^2}\mathrm{,}l=\frac{tR}{t+D}$ (1) Then, for the detector coordinate r, the incident neutron flux Φ_sc (E) and the scattered neutron flux of the detector Φ_st (E) can be expressed in the following form. $\begin{array}{l}\text{d}{\text{Φ}}_{\text{sc}}\left(E\right)=\frac{\text{d}{\text{Φ}}_{\text{st}}\left(E\right)}{4\pi R^2}\cdot {\text{e}}^{-{\Sigma }_{\text{t}}\left(E\right)l}\\ \text{d}{\text{Φ}}_{\text{st}}\left(E\right)={\text{Φ}}_0\left(E\right)\cdot {\text{e}}^{-{\Sigma }_{\text{t}}\left(E\right)\left(T-t\right)}\cdot {\Sigma }_{\text{s}}\left(E\right)\text{d}t\end{array}$ (2) Here, Φ₀ (E) is the neutron flux of the source, ∑_s(E) and ∑_t(E) are the scattering macroscopic cross section and the absorbing macroscopic cross section of the sample, respectively. ∑_det(E) is the absorption cross section of the detector material.

We have to take into account the neutron attenuation effect, where a single neutron is detected with a probability of $P=1-{\text{e}}^{-\Sigma d/\cos \theta }\mathrm{.}$ (3) The neutron flux at the detector surface is $\begin{array}{ll}\text{d}{I}_{\text{det}}\hfill & ={\displaystyle {\int }_E}P\cdot \text{d}{\text{Φ}}_{\text{sc}}\left(E\right)\text{d}E\hfill \\ \hfill & ={\displaystyle {\int }_E}{\text{Φ}}_0\left(E\right)\frac{D+t}{4\pi R^2}\cdot {\text{e}}^{-{\Sigma }_t\left(E\right)l}{\text{e}}^{-{\Sigma }_t\left(E\right)\left(T-t\right)}\hfill \\ \hfill & {\Sigma }_s\left(E\right)\cos \theta \left(1-{\text{e}}^{-{\Sigma }_{det}\left(E\right)d/\cos \theta }\right)\text{d}t\text{d}E\hfill \\ \mathrm{.}\hfill & \hfill \end{array}$ (4) Then, according to the definition of the point scattered function we have $\begin{array}{c}\text{PScF}\left(r\right)=I_{\text{det}}/{\displaystyle {\int }_E}P\cdot {\text{Φ}}_0\left(E\right)\text{d}E\\ =\frac{{\displaystyle {\int }_{E\mathrm{,}t}}{\text{Φ}}_0\left(E\right)\frac{D+t}{4\pi R^{3/2}}{\text{e}}^{-{\Sigma }_t\left(E\right)\left(\frac{t}{t+D}R+T-t\right)}\left(1-{\text{e}}^{-\frac{{\Sigma }_{det}\left(E\right)\cdot d}{\cos \theta }}\right){\Sigma }_s\left(E\right)\text{d}t\text{d}E}{{\displaystyle {\int }_E}{\text{Φ}}_0\left(E\right)\cdot \left(1-{\text{e}}^{-{\Sigma }_{det}\left(E\right)\cdot d}\right)\text{d}E}\mathrm{.}\end{array}$ (5) This function is quite complicated and inconvenient to use. Many studies related to scattering correction have been conducted to approximate this function.

For a very thin detector, it can be considered that $d\to 0$. Because the series expansion of the exponential function ${\text{e}}^{-{\displaystyle \sum d}}=1-{\displaystyle \sum d}+\mathcal{O}\left[2\right]$ involves a second-order infinitesimal term $\mathcal{O}\left[2\right]$, we can use this approximation by considering that the main and first-order terms occupy 99% of the exponential function (i.e., d satisfies 1-∑d=0.99e^-∑d). This implies that the following approximation can be applied. $\begin{array}{ll}\hfill & 1-{\text{e}}^{-\frac{{\Sigma }_{det}\left(E\right)d}{\cos \theta }}\to \frac{{\Sigma }_{\text{det}}\left(E\right)d}{\cos \theta }\hfill \\ \hfill & 1-{\text{e}}^{-{\Sigma }_{det}\left(E\right)d}\to {\Sigma }_{\text{det}}\left(E\right)d\hfill \end{array}$ (6) Using the geometric relations and the adjustment order, and writing the new function as TPScF (r), we obtain $\begin{array}{l}\text{TPScF}\left(r\right)=\frac{{\displaystyle {\int }_{E\mathrm{,}t}}{\text{Φ}}_0\left(E\right)\cdot {\Sigma }_{\text{det}}\left(E\right)\cdot {\Sigma }_{\text{s}}\left(E\right)\cdot g\left(E\mathrm{,}t\right)\text{d}t\text{d}E}{{\displaystyle {\int }_E}{\text{Φ}}_0\left(E\right)\cdot {\Sigma }_{\text{det}}\left(E\right)\text{d}E}\\ g\left(E\mathrm{,}t\right)=\frac{1}{4\pi \sqrt{R}}{\text{e}}^{-{\Sigma }_{\text{t}}\left(E\right)\left(\frac{tR}{t+D}+T-t\right)}\mathrm{.}\end{array}$ (7) The integral in the numerator of Eq. 7 with respect to t can be written as $\begin{array}{ll}\hfill & F\left(E\mathrm{,}r\right)={\displaystyle {\int }_0^T}f\left(E\mathrm{,}r\mathrm{,}t\right)\text{d}t\hfill \\ \hfill & ={\displaystyle {\int }_0^T}\frac{1}{4\pi \sqrt{R\left(r\mathrm{,}t\right)}}{\text{e}}^{-{\Sigma }_t\left(E\right)\left(\frac{tR\left(r\mathrm{,}t\right)}{t+D}+T-t\right)}\text{d}t\mathrm{.}\hfill \end{array}$ (8) This integral can be calculated numerically using the Newton-Cotes method after the associated parameters are given. Finally, the thin-film approximation of the model is $\text{TPScF}\left(r\right)=\frac{{\displaystyle {\int }_E}{\text{Φ}}_0\left(E\right)\cdot {\Sigma }_{\text{det}}\left(E\right)\cdot {\Sigma }_s\left(E\right)\cdot F\left(E\mathrm{,}r\right)\text{d}E}{{\displaystyle {\int }_E}{\text{Φ}}_0\left(E\right)\cdot {\Sigma }_{\text{det}}\left(E\right)\text{d}E}\mathrm{.}$ (9) Here, the integral of E is related to the specific energy spectrum.

Pre-processing of the model

To better apply this model, we need to pre-process the model, as described in this section. We focus on Eq. 8 and apply the Newton-Cotes method to obtain a numerical solution [14] using $\begin{array}{c}F\left(E\mathrm{,}r\right)\approx \frac{h}{3}(f\left(E\mathrm{,0}\right)+f\left(E\mathrm{,}T\right)\\ +4\cdot {\displaystyle \sum _{n=1}^m}f\left[E\mathrm{,}\left(2n-1\right)h\right]\\ +2\cdot {\displaystyle \sum _{n=1}^{m-1}}f\left[E\mathrm{,2}nh\right]),\end{array}$ (10) where h represents the step size and m is obtained by solving the following equation. $\left(2n-1\right)h_{\text{m}}=T\mathrm{,}m=\left[h_{\text{m}}\right]$ (11) Here, [x]represents the Greatest Integer Function, which returns the value of the greatest integer, which is less than or equal to x. The error that results from Eq. 10 is given by the following equation. $\text{Error}\left(x\right)=\frac{T{h}^4}{180}\frac{{\partial }^4}{\partial t^4}f\left(E\mathrm{,}\xi \right)\mathrm{,}\xi \in \left(\mathrm{0,}T\right)$ (12) Eq. 10 can be calculated very quickly once the determined thickness T of the sample and the absorption cross section ∑t(E) are known.

In addition, the integration of TPScF (r) is difficult. Therefore, the approach we adopt is to divide the integration interval into many small intervals as follows. $\text{TPScF}\left(r\right)=\frac{{\displaystyle {\sum }_{i=0}^{E_{\text{m}}}{\text{Φ}}_0\left(E_i\right){\Sigma }_{\text{det}}\left(E_i\right){\Sigma }_{\text{s}}\left(E_i\right)F\left(E_i\mathrm{,}r\right)}}{{\displaystyle {\sum }_{i=0}^{E_{\text{m}}}{\text{Φ}}_0\left(E_i\right){\Sigma }_{\text{det}}\left(E_i\right)}}$ (13) here, $E\in \left(E_0\mathrm{,}{E}_{\text{m}}\right)$ is the range of values of the energy spectrum of the neutron source.

Thus, we consider the model from a numerical perspective, making it less abstract and more computable in concrete terms. Next, we use examples to demonstrate the validity of the model and then perform an error analysis.

Numerical simulation of scattering correction for water and polyethylene

We construct a simple neutron radiography experiment. The neutron converter is made of 100 μm-thick ⁶Li material, and cylindrical metallic aluminum with a radius of 0.6 cm and a height of 1.1 cm is used as the sample, which is located 0.8 cm away from the detector. The experimental setup is shown in Fig. 2.

Fig. 2Full-screen View

Neutron radiography experiment setup.

As the neutron source in the experiment, we use thermal neutrons ranging in energy from 0.01 eV to 1 eV and whose energy spectrum satisfies Maxwell’s spectrum. The neutron distribution on the converter can be obtained by importing the relevant parameters into MCNP. The results are represented as density maps for PScF (r) and TPScF (r), in Fig. 1.

In the figure, r indicates the distance from the detection point to the center of the detector plane. Fig. 3(a) represents the simulation results using PScF (r) as a point scattered function, and Fig. 3(b) represents TPScF (r) for comparison purposes. The red areas in Fig. 3 indicate heavy neutron scattering, whereas the blue areas indicate very light neutron scattering. The boundary of the red area lies roughly at 0.6 cm, which is generally in line with the actual condition of the sample. Note that the boundaries of the red areas in each figure are not particularly clear. The causes behind this are the idealization of the neutron energy spectrum and the haphazard nature of the scattering in the sample. However, there are well-established methods to solve this problem, mainly by correcting the energy spectrum and by using computer algorithms to make the edges clearer, thus making neutron imaging more representative of the actual experimental conditions [15, 16]. Comparing Fig. 3(a) with (b), it can be seen that TPScF (r) gives an overall smaller result than PScF (r), which indicates that an overall error exists between them. The Relative Error Function is therefore defined as $\sigma \left(r\right)=\frac{\left|\text{TPScF}\left(r\right)-\text{PScF}\left(r\right)\right|}{\text{PScF}\left(r\right)}\mathrm{.}$ (14) For the experiment described in Fig. 2, the image of the relative error function is mapped in Fig. 4

Fig. 3Full-screen View

(Color online) Results of (a) PScF (r) and (b) TPScF(r) for the simulation experiment

Fig. 4Full-screen View

(Color online) Results of the relative error distribution between PScF (r) and TPScF (r).

This result indicates that the error is very small in the central area of the screen and is larger along the edge areas. The average error over the entire detection screen is approximately 0.306. However, this model is valid because of the similarity between Figs. 3(a) and (b), so, we further review the validity of the model through error analysis below.

Error analysis

The error function presented in Fig. 4 resembles an exponential function in the radial direction, so we can introduce an exponential function to adjust for the error. First, we define $\text{NPScF}\left(r\right)=\text{TPScF}\left(r\right)+C{\text{e}}^{-ar}$ (15) as the new approximate point scattered function. By fitting the relative error data, the values of the parameters C, a can be obtained, and then we have $\text{NPScF}\left(r\right)=\text{TPScF}\left(r\right)+8.7646\times {10}^{-5}{\text{e}}^{0.3980r}\mathrm{.}$ (16) Then, to study the error relative to the original scattered function, the distribution is plotted in Fig. 5.

Fig. 5Full-screen View

(Color online) Relative error distribution between PScF (r) and NPScF (r).

The figure indicates that the relative errors are at a very low level in most regions, and a maximum value of 17% and a mean value of 2.5% exist across the entire region. This result indicates that a further improvement can be made to the approximation model, which validates that NPScF (r) is a very effective point scattered function for the thin-film approximation.

Quantitative neutron photographic analysis method

Owing to the nature of the interaction between neutrons and matter, the results of neutron imaging are often very inaccurate in actuality. We can reduce the error by using more advanced experimental equipment, and/or we can improve the imaging technique using algorithms to correct an image. One such method, the Quantitative Neutron Photographic Analysis Method uses algorithms to improve imaging technology.

Quantitative neutron photographic analysis uses the exponential neutron decay pattern, $I_{\text{sample}}\left(i\mathrm{,}j\right)=I_{\text{o}}\left(i\mathrm{,}j\right){\text{e}}^{-\left({\Sigma }_{\text{A}}T\left(i\mathrm{,}j\right)\right)}\mathrm{,}$ (17) where I is the intensity of the neutron, (i, j) is the prime point in the image, and ∑_A is the thickness of the sample. The distribution of substances can be determined using Eq. 17. $T\left(i\mathrm{,}j\right)=-\frac{\ln \left(\frac{I_{\text{sample}}}{I_{\text{o}}}\right)}{{\Sigma }_{\text{A}}}$ (18) This is the principle behind quantitative neutron photographic analysis method.

Numerous experimental studies have shown that quantitative neutron photographic analysis is influenced by three main factors.

(1). Collimation Ratio

The collimation ratio is the ratio of the distance from the neutron source to the sample with respect to the diameter of the output aperture. Because neutrons are electrically neutral, they cannot be made parallel by focusing the neutron beam through the electromagnetic field, and therefore it is not possible to achieve full collimation of the neutron beam using a collimator.

(2). Scattered Neutrons

The presence of substances with high scattering cross section (e.g., water and polyethylene) in the sample can cause a large number of scattered neutrons to reach the detector, thus distorting and blurring the neutron image. The solution is generally to increase the distance between the sample and the detector.

(3). Neutron Energy Spectrum

Because incident neutrons are not mono-energetic and neutrons having different energies interact with matter having different cross sections, the images and results obtained using a single cross section parameter can produce large errors. In addition, the energy spectrum of the neutrons scattered by the sample may have changed from that of the neutron source owing to attenuation during propagation.

The focus of this study is the correction of neutron images by applying a correction for scattered neutrons. To further investigate the optimization of the proposed model with respect to neutron imaging, the current iterative algorithm for the quantitative analysis of scattered samples is presented in Fig. 7.

Fig. 7Full-screen View

Neutron radiography quantitative analysis iterative algorithm flow chart

In the algorithm, I_su is the normalized neutron intensity, I_sc is the scattered neutron distribution function, and k is the number of iterations. The criterion equation $\frac{{\displaystyle {\sum }_i{\displaystyle {\sum }_j\left|I_{\text{su}}^{k+1}\left(i\mathrm{,}j\right)-I_{\text{su}}^k\left(i\mathrm{,}j\right)\right|}}}{{\displaystyle {\sum }_i{\displaystyle {\sum }_j{I}_{\text{su}}^k\left(i\mathrm{,}j\right)}}}⩽\epsilon$ (19) is used to select a convergence coefficient ε that terminates the iterative calculations. The above criterion actually represents the relative change in the transmitted neutron intensity calculated for two adjacent iterations. As the algorithm converges, the value of the change should become increasingly smaller, gradually converging toward zero. This algorithm has been tested in a large number of neutron scattering experiments and is highly reliable. This iterative algorithm can be applied to the thin-film point scattering function model developed above.

Optimization algorithm based on thin-film approximation models

It should be noted that in the iterative algorithm (Fig. 7), PScF (r) is calculated based on the sample thickness, which is later brought into the calculation of the neutron intensity prior to the iterative calculations. However, if we use NPScF (r) instead of PScF (r), the iterative convergence can be greatly accelerated, as NPScF (r) is equivalent to an approximate numerical integration of PScF (r). When we applied the iterative algorithm to the experiments in Section 3.1; using this substitution only three iterations were required to obtain satisfactory results, as indicated in Fig. 6.

Fig. 6Full-screen View

(Color online) Scattering imaging results using ((a) NPScF (r), and (b) PScF (r))

The results obtained show the reliability and validity of the iterative algorithm on the one hand, and the faster correction of imaging results using NPScF (r) instead of PScF (r) on the other. The iterative calculation using PScF (r) requires five iterations to converge to something closer to actuality and takes approximately three times as long to run as the calculation using NPScF (r). Comparing result (a) using the NPScF (r) function with result (b) using the conventional PScF (r), it can be observed that (b) produces poor imaging results, is slower to converge, and exhibits greater distortion at the edges, whereas (a) produces good imaging results and slight distortion at the edges.

We followed up with similar simulations for trapezoidal objects and cones and found that using NPScF (r) instead of PScF (r) for iterative calculations not only improved the simulation, but also improved the speed of the algorithm.

Moreover, in recent years, a “3D detector” has been introduced for use in neutron photographic techniques [18, 19]. It adds 3D microstructures such as trenches or holes to a planar detector and fills them with a neutron converter, which can significantly improve the detection efficiency. For this class of detectors, the model proposed in this paper is also applicable, and together with the algorithm in Fig. 7, yields optimized results. However, the relationship between the geometric parameters of the planar detector (Fig. 1) must be recalculated. That is, to apply the proposed model to the 3D detector, the relationship between geometric parameters must be solved again. A mathematically feasible approach is to use the symmetry of the 3D detector microstructure to solve by parts. For example, if the 3D microstructure is axisymmetric, then it can be solved by the radius with the center of symmetry being the center of a circle, determining the geometric relationship at a fixed radius, and then solving using the method given in Sect. 2.

Effects of energy spectrum hardening

For the thermal neutron spectra considered in this study, the cross sections interacting with matter approximately obey the 1/v law. As neutrons penetrate the sample, low-energy neutrons decay more rapidly than high-energy neutrons, ultimately increasing the proportion of higher-energy neutrons in the neutron energy spectrum, a phenomenon known as energy spectrum hardening.

The current solution to this problem is to define an "effective neutron cross section to correct for the energy spectrum. Considering that the macroscopic cross section is no longer a constant, the neutron transmittance is no longer simply exponential in relation to the thickness of the sample, even in the absence of scattered neutrons. Therefore, the energy response of the detector should also be considered. We can define an effective cross section based on the ratio of the incident neutron intensity to the transmitted neutron intensity as $\frac{I_{\text{su}}}{I_{\text{o}}}={\text{e}}^{-{\Sigma }_{\text{eff}}\left(T\right)\cdot T}\mathrm{,}$ (20) where -∑_eff(T) represents the effective cross section for a sample of thickness T.

The intensity of the neutrons detected by a detector is typically expressed by the following equation: $I_{\text{o}}={\displaystyle {\int }_E}\text{Φ}\left(E\right)\left(1-{\text{e}}^{-{\Sigma }_{\det }\left(E\right)\cdot d}\right)\text{d}E\mathrm{,}$ (21) and the transmitted neutron intensity is $I_{\text{su}}={\displaystyle {\int }_E}\text{Φ}\left(E\right){\text{e}}^{-\Sigma \left(E\right)\cdot T}\left(1-{\text{e}}^{-{\Sigma }_{\det }\left(E\right)\cdot d}\right)\text{d}E\mathrm{.}$ (22) According to the 1/v law, which is approximately satisfied by thermal neutrons, the conventional neutron cross section can be expressed as $\Sigma \left(E\right)={\Sigma }_{\text{th}}\frac{\sqrt{E_{\text{th}}}}{\sqrt{E}}\mathrm{,}$ (23) where $E_{\text{th}}\approx 0.025\text{eV}$ is the energy of a thermal neutron. Considering the detector energy response as a whole, $\epsilon \left(E\right)=1-{\text{e}}^{{\displaystyle {\sum }_{\det }\left(E\right)\cdot d}}\approx {\displaystyle {\sum }_{\det }\left(E\right)\cdot d}$, and combining the above equations yields ${\Sigma }_{\text{eff}}\left(T\right)=-\frac{1}{T}\ln \left(\frac{{\displaystyle {\int }_E}\text{Φ}\left(E\right){\text{e}}^{-{\Sigma }_{\text{th}}\frac{\sqrt{E_{\text{th}}}}{\sqrt{E}}}\epsilon \left(E\right)\text{d}E}{{\displaystyle {\int }_E}\text{Φ}\left(E\right)\epsilon \left(E\right)\text{d}E}\right)\mathrm{.}$ (24) However, owing to the difficulties of performing integration operations and in obtaining relevant data, use of a simplified approximation of the above equation is desirable. The approach taken in [17] is to apply a linear fit, but we believe that the error in doing so would be very large. Instead, we can adopt a numerical integration approach, as described in Section 2.2 to approximate the above equation. Again applying the Newton-Cotes method, let ${\Sigma }_{\text{eff}}\left(T\right)\approx \frac{1}{T}\ln \left(\frac{4\cdot {\displaystyle {\sum }_{n=1}^m\text{Φ}\left(\left(2n-1\right)h\right)}{\text{e}}^{-{\Sigma }_{\text{th}}\frac{\sqrt{E_{\text{th}}}}{\sqrt{E}}}\epsilon \left(\left(2n-1\right)h\right)+2\cdot {\displaystyle {\sum }_{n=1}^{m-1}\text{Φ}\left(2nh\right)}{\text{e}}^{-{\Sigma }_{\text{th}}\frac{\sqrt{E_{\text{th}}}}{\sqrt{E}}}\epsilon \left(2nh\right)}{4\cdot {\displaystyle {\sum }_{n=1}^m\text{Φ}\left(\left(2n-1\right)h\right)}\epsilon \left(\left(2n-1\right)h\right)+2\cdot {\displaystyle {\sum }_{n=1}^{m-1}\text{Φ}\left(2nh\right)}\epsilon \left(2nh\right)}\right)\mathrm{,}$ (25) where h is the set step size, and note that Φ(E)=0 at both E=0 and $E=\infty$. With the introduction of the effective cross section, the use of an effective cross section to correct the energy spectrum can make the decomposition lines of the spectrum more visible, so that the edges of the sample or large variations in density are more clearly reflected in the image.

However, it would be preferable to make certain modifications to NPScF (r), i.e., Eq. 15, to achieve an improved energy spectrum correction.

Maxwell’s spectrum can be formally defined as Φ (E) = aE exp (-E/b), and the hardening of the energy spectrum is essentially an excess retention of fast neutrons and an excess depletion of slow neutrons, and is therefore reflected in an increase in the parameter b. Considering the specific form of Maxwell’s spectrum, an increase in parameter b is equivalent to an increase in the background temperature. Thus, we can perform the following transformations to account for these factors. $\text{Φ}\left(E\right)\to {\text{Φ}}_{\lambda }\left(E\right)=\frac{2\pi n}{{\left(\pi kT\right)}^{3/2}}\sqrt{\frac{2}{m}}E{\text{e}}^{-\frac{E}{k\left(T+\lambda \right)}}$ (26) The temperature-corrected energy spectrum is then used in the step calculations. Parameter λ can then be determined by fitting the specific energy spectrum to Eq. 26.

Substituting the corrected energy spectrum into Eq. 13 and replacing the mean scattering cross section with the effective scattering cross section in Eq. 25 yields $\text{TPScF}\left(r\right)=\frac{{\displaystyle {\sum }_{i=0}^{E_{\text{m}}}{\text{Φ}}_{\lambda }}\left(E_i\right){\Sigma }_{\text{det}}\left(E_i\right){\Sigma }_{\text{eff}}\left(E_i\right)F\left(E_i\mathrm{,}r\right)}{{\displaystyle {\sum }_{i=0}^{E_{\text{m}}}{\text{Φ}}_{\lambda }}\left(E_i\right){\Sigma }_{\text{det}}\left(E_i\right)}\mathrm{.}$ (27) Subsequently, we can correct for the error by applying error analysis, as discussed in Section 3.2, by choosing a suitable exponential function based on the symmetry of the spatial distribution.