Consistent Algorithms Marching Along Characteristics for the Numerical Solution of the Boltzmann Equation

,


Introduction
A line propagation of radiation and particles is observed for low-density noncharged particles as neutron in nuclear reactors and X-rays or γ-rays in nuclear medicine and engineering nondestructive test applications.For charged particles, this collimation is found in accelerator and spallation target technologies for ADS nuclear plants.This natural behavior of collimated radiation and particles introduces in the numerical computation a new paradigm of operator, the one-speed Boltzmann transport operator, which has been extensively studied in the nuclear community Vladimirov  or no spectral dependency in the case of radiation, the following formulation is obtained from the linear Boltzmann equation, Vladimirov 1 ,Kaper et al. 3 , Dautray and Lions 5 , which is usual for the mathematical modeling of the interaction of the radiation with the participating medium.Based on these initial assumptions, we consider the following boundary problem: for a given function g − , find the constant velocity radiation intensity φ ∈ L 2 S n−1 × Ω such that ω • ∇φ ω, x σ t x φ ω, x where Σ ± { ω, σ ∈ S N−1 × Γ; ±ν σ • ω > 0} are, respectively, the phase space surfaces of the incident and emergent radiation on the physical surface Γ ∂Ω, which is the smooth boundary of the bounded, possibly convex region Ω ⊂ R n .
The coefficients σ t and f in 2.1 are, respectively, the total extinction coefficient absorption outscattering and the scattering and fission coefficient.We assume that σ t ∈ L ∞ Ω and, since ω • ω cos θ 0 , with θ 0 being the angle between the direction of incident radiation ω and the emergent scattered radiation ω, we also assume that Under these assumptions, L 2 Σ ± ; |ν σ •ω|dσ dω become the appropriated spaces for traces of the functions in the boundary value problem 2.1 and 2.2 .For a more detailed discussion, see Cipolatti 16 and the references therein.

Mathematical Preliminaries
The stationary transport 2.1 is composed of the following operators.The operator where derivatives are in the sense of distributions.The space H L 2 S n−1 × Ω is the closure of C S n−1 × Ω with respect to the norm where dω denotes the measure on S n−1 associated with the Lebesgue measure R n .The space W is defined from H as W φ ∈ H; Aφ ∈ H .

3.3
W is a Hilbert space for the norm

3.4
The second operator is a bounded operator which has continuous inverse for σ t x > 0, x ∈ Ω.The total macroscopic cross-section is, in applications, piecewise smooth piecewise and infinitely differentiable , bounded, and strictly positive on Ω with 0 < inf{σ t x ; x ∈ Ω} ≤ σ t x ≤ c 0 sup{σ t x ; x ∈ Ω}.

3.6
The streaming operator is the sum

3.7
The third operator describes the anisotropic scattering and fission process.Kernel f can be expanded in an absolutely and uniformly convergent series where P k−1 is the surface harmonic of degree k − 1 Legendre polynomials .The coefficients σ k , k ≥ 1 are nonnegative, and the first coefficient verifies the inequality and 1/λ σ f is the nonnegative fission cross section divided by the criticality parameter.With this, the collision operator takes its usual form

3.11
The operator C with this kernel is compact from and this is the main mathematical problem in this engineering problem.
The trace operator 12 is continuous and surjective, where W is a Hilbert space with the norm Given the incident flux by the boundary condition 2.2 , g − ∈ L 2 Σ − , |ω • ν σ |dσ dω , the following prolongation of g − inside the phase space is defined: where τ ± ω, x sup{t ≥ 0 : x ± tω ∈ Ω}are the distance to the emergent and incident boundary, respectively.Note that this definition implies that J σ is the solution operator of the problem L σ φ 0, φ| Σ − g − .These operators can be used to formulate the first-order boundary non-homogeneous stationary direct problem given 2.1 and 2.2 as the following equivalent homogeneous problem: to find ϕ φ − J g − ∈ W ∩ dom A such that where dom A {ϕ ∈ W : ϕ 0 on Σ − }.
Theorem 3.1.For each h ∈ H, the boundary value problem 3.14 has a unique solution witch is given by and verifies the following estimate: where d is the diameter of Ω.
Proof.See Chapter III and Lemma 3.1 of Vladimirov 1 .
In Sections IV thought VII of Vladimirov 1 , operators for more general Sobolev spaces.We can easily adapt them to our Hilbert space framework.

Lemma 3.2. The operator
This result is fundamental for the reformulation of boundary value problem given by 2.1 and 2.2 in terms of an integral equation based on a compact perturbation of the identity operator used in this work: to find φ ω, x such that 3.17 Remark 3.3.The study of the numerical behavior of compact operators is an very well known subject.The spectral properties of compact operators are very similar to that of finite-dimensional operators, fact that make analysis very easy.For this kind of operator, the subspaces associated with a eigenvalue are finite dimensional.It is based on this behavior our interest in model this problem using integral equation 3.17 .A situation similar also occurs in models based on the Laplacian operator paradigm, which fortunately, for bounded domains, is an operator with compact resolvent.Since the one-velocity transport equation has a more complex structure than the Laplacian, it is necessary to do the composition of the inverse of the streaming operator, L −1 σ , and the collision operator, C, to obtain a compact operator.Note that L −1 σ presupposes and line integrals through the ray paths of particles inside the domain, and is this fact that will have a remarkable influence on the numerical behavior of the problem.We will see that this connects compactness and characteristic basis.

Relation between the Integral Operator and the Characteristic Basis
Due to its very nature, the transport equation for the angular flux may be viewed as a system of ordinary differential equations coupled by the collision operator.Each angular field is a different field which has its own direction of propagation as its characteristic directions.This means what characteristic means: the discontinuities of each angular flux field will separately propagate in its own direction.This observation introduces a natural way to solve the transport equation, Askew 17 , that will be computationally more expensive but will produce more accurate calculations.For this, consider the mapping in the angular phase space This is an isometry that rotates the axis x to a new direction x ω coincident with the propagation direction.For each fixed direction ω ∈ S n−1 , it introduces a new coordinate system x ω , x ω ⊥ with the second coordinate in the projection Π ω Ω of the physical domain Ω in a subspace of R n made with directions ω ⊥ and evidences the ordinary character of each angular flux in the equation subject to the same incident boundary condition Note that this incident boundary condition can be interpreted as an initial value for this system of ordinary differential equations coupled through the collision operator.For simplicity, we have used the same notation φ ω, x ω , x ω ⊥ for the rotated angular flux Q φ ω, x ω , x ω ⊥ .As pointed out in Section 3., the inversion of the operator L σ and its composition with special class of collision operator found in nuclear applications gives a compact operator Of course this important property will pose the direct problem involving recovery of the angular flux from appropriate parameter data σ, σ s , and σ f .Some more functionals defined in terms of the rotated coordinates systems need to be introduced in order to make comprehensive the approximations that will be done to solve the problem.First of all, the optical distance from interior point x to a point where ω x − x 0 / x − x 0 is the direction that aligns x and x 0 , and the extension for points in the possibly bare medium outside the physical domain Ω may be done by assuming σ x 0, x ∈ R n \ Ω.When x 0 is a boundary point, it is given by x ± ω, x x ± τ ± ω, x ω which is the point in the boundary where the line passing x with direction ω crosses the emergent or the incident boundary, respectively.It may be used to define the attenuation from the incident boundary We also define the right inverse of the trace operator, the extension operator as and the extension with attenuation from the incident boundary By noting that the attenuation is the inverse of the integrating factor for the system 4.3 , we may integrate from the incident boundary to the interior point x to obtain an integral equation which is an explicit expression of 3.17 : The integration may also be done between two arbitrary points, inside or outside the domain Ω, x, and x 0 .In this case we must extend the cross-sections outside the domain as zero.By defining ω Remark 4.1.Problems given by the equations 1 2.1 and 2.2 , 2 3.17 , are alternative formulations to the one-velocity transport equation model and express different aspects of this equation.Note that the initial value character of the problem will be accomplished if we construct the solution by using a method that marches along each one of the characteristics.
Remark 4.2.The consistence that will be appropriated to the application of the Lax equivalence theorem framework that has been formulated for time transient problems will be accomplished by the adoption of a rotated pixel basis in the next section.The compact operator L −1 σ C is subcritical or critical, which means the inverse of its eigenvalue is not less than one.So this gives stability to the class of problems been investigated.
As a consequence of consistency, the Lax equivalence theorem applies, Valtchev and Roberty 18 , and may be stated as follows.
Theorem 4.3.For a well-posed initial value problem, a consistent numerical scheme is convergent if and only if it is stable.

The Natural Partition
The integral form of the transport equation 4.11 has been written in a different coordinate system for each direction of propagation ω.In order to attempt the regularities requirements of each angular flux φ, we must introduce a grid with one direction as ω.This is impossible unless we choose a set of discrete ordinates as representative of all angular directions.This must be done a priori.Let S be the surface of the unity sphere in R 3 and where θ and ψ are the usual angles of the spherical coordinates system.To be more specific, let us consider the two-dimensional equation embedded in R 3 .So, we choose a set of of equally spaced angular angles ψ j j − 1 /2J 2π and at least symmetrical with the angle θ π/2 of 2I collocations directions out of the plane generated by the vectors cos ψ , sin ψ , that is, out of the plane xy.We now may use these directions to introduces grids, one for each direction, or in the case of a two-dimensional problem, only in the plane directions Π j {cos ψ j , sin ψ j }.Note that in the two dimensional problem, since the problem is embedded in the three-dimensional space, angular flux is defined for all directions in the S 2 sphere but has spatial variations only in the plane xy.This introduces, for each direction out of the xy plane, the necessity of a sin θ i factor correction.This gives the possibility of extending this formulation to investigate axis symmetric problems, but, from now on, we will fix in the two-dimensional problem.We have introduced 2J coordinate systems.Without loss of generality, let us consider each one with a rectangular equally spaced grid.These coordinates systems are supposed to be equally rotated with respect to the others.We have chose an even number for J, so, directions of propagations appear in groups of four, to complain with the four quadrants.
We know that it is necessary to combine the different directions in order to calculate the collision operator for groups of pair the directions.Since the cross sections are properties common for all angular flux, the complete intercessions of the 2J rotated grid is the appropriate partition consistent with the numerical solution of the transport equation.To fix ideas, let us introduce the appropriate notation in R 2 .The direction ω j contains one strip

5.3
Journal of Applied Mathematics which may be intercept with a strip in one direction in ω ⊥ j χ m j x, y to form a rotated pixel 0, else.

5.5
The regular spacing h x j , h y j 1 has been adopted without loss of generality for simplicity.Note that in the direction ω j cos ψ j e x sin ψ j e y , the integrated 4.11 must be collocated for points ω j , x 0ω j , x ⊥ ω j and ω j , x ω j , x ⊥ ω j that will be referenced by the indexes j, m j , n j and j, m j 1, n j , respectively.This is a consequence of the fact that ω j aligns the points and x ⊥ 0ω j x ⊥ ω j .As the marching direction is ω j , the equation index m j will grow from its value corresponding to a point x ω j in the incident boundary, where the boundary condition is imposed, to the corresponding point in the emergent boundary.This happens for each transversal index n j in each direction ω j .We now face the problem of computing the collision integral in 4.11 with a consistent discretization.Note that the measure dω dX ω will corresponds to a volumetric Lebesgue measure when the two integrals interchanges, that is, when the Fubini theorem is applicable.The natural way to obtain this in the discretized approximation is by considering the partition presented by Roberty 19 that makes a complete intersection of the pixels P jm j n j .This may be done with the complete interception of the strips defining the pixels to form the natural pixels χ n j x, y .

5.7
This interception produces a domain partition which is consistent in the sense that if we increase the number of directions and reduce the size of the spacing h j 1 by using a more refined mesh, the error in the approximation to the exact stationary transport equation will approache zero accordingly.A particle or photon collected in one interception can only come from pixels aligned with the current position thought the discrete directions considered in the model.The approximation made can answer the question "from where the particle came?", in a way consistent with the exact equation.Calculations of the vertex, area, centroid, and all geometric information are a classical computational geometry problem of domain decomposition, O'Rourke 20 , which may present a great complexity in dimensions greater than two.
As an optimization problem, the mesh generation is the determination of the convex polygonal sets limiting the maximum area that follows linear restrictions of the type: Ax ≤ b, where x is a vector in R 2 , b is a vector in some R N b related with the strips n j , and A is a 2 × N b matrix related with the direction angle θ j .The search process of the elements in the mesh executes the following steps: 1 decomposition of the domain in search subdomains 2 determination of the minimum and maximum strips for the search subdomain 3 search for the true intersections, 4 calculus of elements vertices, edges, area, and pertinence index for all directions 5 assemble subdomains to form the natural pixel partition.

Decomposition of the Domain in Search Subdomains
The greater majority of the combinations n 1 , . . ., n 2J are false intersections and can be a prior avoided by restricting the search for parameters on regions for which we can expect true intersections.We call these regions search subdomains, and they can be a pixel for some prescribed direction j, an intersection of some given strip with a sector, or any other type of polygonal subdomain.With them, we avoid a search with a large number of false intersections in the step 3 .As an example in Figure 1 we show a detail of the mesh J64M128 for huge numbers of direction 64 and pixels 128 × 128 , which is typical of image processing and reconstruction problems.Since the rotational symmetry allows us to calculate the interception only in a J/π sector, we may have less computation than the case of doing only two interceptions, as illustrated by Figure 2. Finally, in Figure 3, we present the complete mesh used to simulate our solution with test and benchmark problems.

Determination of the Minimum and Maximum Strips for the Search Subdomain
This is a complementary procedure for minimizing the determination of false intersection.

Search for the True Intersections
In each domain of search, every strip combinations must be tested until the area of the domain is totally filled with natural pixels.

Calculus of Element Vertexes, Edges, Area, and Pertinence Index for all Directions
This calculation is done simultaneously with the search for true intersections.Note that not all intersections of lines limiting the strips are vertexes of the polygons and not all lines across a vertex are a polygonal edge.

Assemble Subdomains to Form the Natural Pixel Partition
In each subdomain, the pertinence indexes are calculated independently and are assembled to form the matrix of pertinence index for the entire mesh.For more information, please see Roberty 19 and in the references there.In Figure 4, the dependence of the numbers of polygonal with the numbers of views and slices is presented.Note that the almost exponential growth observed for growing numbers of slice and views introduces memories problems in the computational implementation of the present methodology.In Figure 5, we present a comparison of mesh J6M3 with the classical Delaunay triangulation made for the same vertex number.Note that it is inconsistent with the natural partition, that is, it is not a triangular decomposition of it.

The Discrete Equation for the Two-Dimensional Problem
As we mention before, we may consider the integral form of the transport 4.11 for direction ω j , j 1, 2J, in the pixel P jm j n j , n j 1, 2M and m j 1, 2M.Since we are focusing on the two dimensional problem, it necessary introduces the sin θ i correction for each out of the xy plane direction.Note that the xy plane corresponds to the case sin θ i sin π/2 1.There is by hypotheses no variation of the angular flux in the z xy ⊥ direction, but it still remains dependent of azimuthal angle.Let us rebuild the notation to make this behavior explicit.
Since ω 3 • ∇ x φ ∂φ/∂z 0 and all other functional variations in this direction are also zero, we rewrite ω ji sin ψ i • ω j and obtain by straightforward calculations, the optical distance for points with projection xy along the ω 3 direction as α xy, xy 0 , ψ i xy•ω j xy 0 •ω j σ xy ω , xy ω ⊥ xy, xy 0 dx y ω j sin ψ i .

6.1
Since all other parameters, as the cross sections, and sources are linear measures of the change in the particle population in the direction ω ji , it must be corrected accordingly.
As a consequence of the divergence theorem, the angular flux in direction ω ji has average value on all lateral segments surfaces in three dimension , that is, the segments that are transversal to the strip χ m j x, y , parallel to ω j for all m j , are equal to zero.We complement this by introducing the following hypotheses: 1 that angular flux is well represented by its average value in each transversal section to the strip χ n j x, y , 2 the cross sections are known by its average values in the natural pixel 3 the average cross sections inside each rotated rectangular pixel P jm j n j are the area volume in three dimensions average weight of each natural pixel cross section value, 4 to accomplish with high cross section values, we integrate the exponentials inside the rotated pixel P jm j n j , 5 the mean value of angular flux inside the pixel P jm j n j is used in the calculation of scattering and fission reaction rate.

The Method of Successive Approximations
With respect to the intensity of the extinction cross section σ t , we observe two important class of problems: 1 there exists a transillumination of the medium Ω, which means that the extinction cross sections not sufficient to produce an optical path that shields the radiation flux.In this class of problems, a flux radiation in the incident boundary propagates through the medium and can be collected in the emergent boundary.Its modification will be used to interpret the optical properties of the medium inside Ω, Since the extinction coefficients may present high values, we considered exponential flux attenuation inside the rotated pixel and calculated the reaction rates using the average flux expressed in terms of the incident and emergent flux in the rotated pixel.Equation 4.11 becomes φ m j 1,n j ,j φ m j ,n j ,j E 1,m j ,n j ,j q m j ,n j ,j E 2,m j ,n j ,j for m j , n j 1, 2M; j 1, 2J. 7.1 Here, the source with contribution of the scattering and which forms fission inside the rotated pixel is , E 4,m j ,n j ,j,i −1 σ m j ,n j ,j exp −σ m j ,n j ,j σ m j ,n j ,j 1 − exp −σ m j ,n j ,j , σ m j ,n j ,j e∈P jm j n j a e σ e .

7.2
Remark 7.1.Note that the correction for extending this methodology for axis symmetric problems is σ e,i σ e sin ψ i 7.3 as mentioned in the introduction of the natural partition in Section 5.
The localization of the natural element e n 1 ,...,n 2J inside the rotated pixel P m j ,n j ,j is made by the incident matrix.As we mentioned before, given one pair e, j , one may calculate the rotate integer coordinates m e, j , n e, j , j by extracting the integer part of the natural element centroid in each rotate coordinate system.This is a procedure similar to the retro projection that is done in image reconstruction algorithms.The function relating the two integer coordinate systems is a discrete version of transformation 4.2 .It is one-to-one between e, j → m, n, j .It is not convenient to make it onto since the implementation of an procedure for the localization of the boundary of Ω will be computationally expensive.

Journal of Applied Mathematics
We will start the march in the algorithm in the position labeled with m j 1, which may correspond to a pixel outside the domain, and march in one of the directions j until the pixel labeled 2M, which also may be outside the domain, is arrived.This introduces no problem if we take care in extending outside the domain the cross sections value by zero, the flux by its boundary value and calculates the exponentials accordingly.
As mentioned before, these equations are consistent and introduce no spatial error when the properties are constant; in that case only, the angular discretization error will be presented.The simulation problems that we will solve in Section 8.3 have material parameters very close to this situation.
As we pointed out before, the marching scheme equation 7.1 is a discretization of the first-order form of the first order system representing the stationary transport equation that has initial value in the incident boundary.The fact that the total cross-section is bigger than the scattering assures contractility of the compact operator related to the scattering kernel, and from this, we will obtain the stability and as a consequence a convergent scheme.

Algorithms Based on Successive Approximation Method
Properties of collision operator equation 3.11 and of the boundary conditions allow us to distinguish the basic class of problems for the stationary transport equation: 1 The collision operator and the external source localized inside Ω are negligible.In this case, the equation is integrable, and we have the direct problem for the X-ray transform.
2 The fission part of the collision operator L −1 σ C fission has the first positive eigenvalue greater than one, and the domain Ω is subcritical external source problem.The incident flux may be zero or not.There are two important cases of nonfission domain, that is, when there is no fissile material inside the domain and σ f is properly zero.The external source inside the domain problem, with q / 0 and g − 0 and transillumination problem with S 0 and g − / 0. The transillumination problem has a greater interest for tomography and nondestructive testing.
3 The fission part of the collision operator L −1 σ C fission is defined with support inside Ω, and we want to find his first positive eigenvalue and its positive eigenfunction.Once we find the eigenvalue, we give to κ the same value and obtain a new fission operator with first eigenvalue equal one, which means that we have found the critical material composition for the domain Ω.

The External Source Inside the Domain Problem
The algorithm for the angular flux calculation in the external source inside the domain problem follows the following steps.
1 Establish an error tolerance criteria.
2 Input material properties and source: σ t , σ s , S. it with the fission operator.
14 For inner iterations equal 1 to maximum number of inner iterations do.
15 Take care to zero the scattering reaction rate, but not the reaction fission rate, since it will not change inside the internal loop.
16 Calculate scattering reaction rate based on the previous inner iteration calculated flux.
17 Calculate the new internal iteration flux with 7.1 with zero external source and zero incident flux and external loop value of the fission reaction rate.
18 Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the inner iteration before the maximum number of inner iterations.Note that we may alternatively use the fission reaction rate error as criteria.
19 End inner iteration do.
20 Take the convergent inner loop value as the outer loop flux.
21 Calculates the new eigenvalue value as a ratio between the fission rate based on the actual outer iteration flux and the fission rate based on the previous outer iteration flux.
22 Verify if the difference between the old eigenvalue and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the outer iteration before the maximum number of outer iterations.
23 End outer iteration do.
Remark 7.2.The first and the inner part of the second algorithm are based on the contracting behavior of the scattering collision operator, that is, we have a compact operator with spectral radius less than one and the rate of convergence will be proportional to this ratio.Perhaps other rates of convergence may be found with other algorithms.The outer loop of the second algorithm is also based on compactness.The problem is to determine the first eigenvalue in the pencil problem is, a perturbation of the identity by a compact operator in a compact weigh operator pencil problem.Its convergence with a decreasing sequence of eigenvalues may be proved by using the fact that the multiplicative region inside Ω has compact support and the composite operator is positive and compact.The fact that the first eigen angular flux is positive almost everywhere in S n−1 × Ω is a classical result based on Jentzsch's theorem for compact operators which leaves nonnegative those functions that are non negative.
Once again the consistence has a positive consequence, since numerical artifacts have not been introduced in the approximation, and good properties of the continuous equation are preserved.

The Transillumination Test Problem
In order to start the evaluation of the method proposed, we first study a subcritical transillumination test problem with no fission and no external sources.The flux for test is φ test exp −0.0125 • y 2 0.05 • x − M for j 24 and zeros for other views.

8.1
The algorithm is used to solve problem 2.1 with source A high-contrast inclusion problem as shown in Figure 6 in a 24 views with 30 × 30 rectangular pixels/view model is solved with the algorithm proposed here which is based on 7.4 accordingly.Material properties are given in Table 1.As can be seen from Figure 7 comparing the given and the reconstructed solution showa that they are almost equal.From Figure 8 showing the absolute error for two views 9 and 24, the error is of order 10 −8 , bigger in the region with the high-contrast inclusion, but with the same magnitude.The absolute error variation with J varying between 6 to 16 and M from 10 to 24 is showed in Figure 9.

A Transillumination Problem
The same transillumination problem with no fission and no external sources is now investigated.The incident radiation boundary condition is g n 2J , 2J 1, for n 2J 1, . . ., 2M and g n j , j 0 for all others views.Figure 10 shows selected particles radiation flux for j 9 and j 24.We can see the depleted effect in the main direction of the incident flux j 24 due to absorption and outscattering and the buildup effect in the other direction j 9 due to the scattered flux.Also, it is possible to note the differences in magnitudes in the flux for the incident direction, j 24, with magnitude order of 1 and the other directions where there is no incident radiation, with magnitude order of 0, 02, one up to 50-times the other.Another important information that can be determinated by the solver presented here characterizes the graph for the albedo operator the incident to emergent flux mapping  In the inverse problem, we ask if it is possible to determine the coefficients functions σ t and f from the a priori knowledge of the albedo operator, that is, data such as that found in Figure 11 for different incident radiation g − .Finally, the convergent sequence {φ iter ; iter 1, . . ., 31, . ..} is tested as a Cauchy sequence in the complete Hilbert space H. Figure 12 shows that before the eleventh iteration, we achieved convergence of order 0, 008.If we claim for 9, 8e − 008, it is achieved in the model problem in the 31 iteration.

Reactor with the EIR-2A and EIR-2B Data
The EIR-2A and EIR-2B are two problems found in Mordant 21 with data about values on cross sections, source and geometric domain definition presented in Table 2 and represented in Figure 13.The total and scattering cross-section is the same for the two problems.The source values is for the EIR-2A problem, and in this case the fission cross-section is zero.The   2. fission cross-section is for the EIR-2B problem, and in this case the source term is zero.Results obtained are in accordance with the theory and are showed in figures.In all cases, the number of axial directions is I 1. Figure 14 shows the angular flux for direction j 1 in the source problem. in the first eigenvalue function, where a not-so-high convergence rate is observed.Finally, Figure 17 shows the first eigen angular flux for direction j 1.These results showing a monotonically decreasing sequence of eigenvalues with the smallest accumulation point at 1 and a nonnegative angular flux are in complete agreement with theoretical results pointed out in Remark 7.2.We must finally note that the continuity in the direction perpendicular to ω 1 has been artificially introduced by the graphic algorithm used to present the results.In reality, the calculations does not impose this restriction on the function.

Conclusions
The present methodology utilizes the natural basis concept, originally developed for image reconstruction, in the numerical solution of very classical and technological problems found  in reactor physics calculations.The algorithms presented here may be considered as giving almost exact solutions in case where material properties and sources are constants in a small number of distinct regions inside the domain, as observed in the simulations EIR-2A and EIR-2B implemented.It does not array system matrices, and as a consequence few bigger problems signify more processing time.The algorithms are deduced by introducing, as less as possible, artifacts related with preprocessing the transport equation by doing inconsistent averages and simplifications hypotheses.It respects the basic characteristic of the transport equation, that is, no continuity requirements in the directions transverse to the propagation direction.It is also based on a scheme consistent with the exact equation.The algorithm implementation also shows compatibility with the present status of computer development.

Figure 1 :
Figure 1: A detail in the J32M64 mesh.

Figure 2 :
Figure 2: Sector mesh generated strip by strip in the j 16 and M 15 problem.

5 HFigure 4 :Figure 5 :
Figure 4: Numbers of polygonal dependence with the numbers of views and slices.

Figure 7 :Figure 8 :
Figure 7: Given and reconstructed view 24 angular test flux in the model.

Figure 9 :
Figure 9: Maximum error variation for J and M.

Figure 10 :
Figure 10: Selected radiation particles flux for illumination on view 24 in a high-contrast nonhomogeneous problem.

Figure 11 :
Figure 11: Emergent flux for illumination on view 24 in the high-contrast model problem.

Table 1 :
Absorption and scattering coefficient in the transillumination model.

Table 2 :
Cross sections and sources values for EIR problems.