GALERKIN METHODS FOR PRIMARY ION TRANSPORT IN INHOMOGENEOUS MEDIA

This paper concerns the energy deposition of high-energy (e.g., ≈ 50 − 500 MeV) proton and carbon ions and high-energy electrons (of ≈ 50 MeV), in inhomogeneous media. Our goal is to develop a flexible model incorporated with the analytic theory for ions based on bipartition and FokkerPlanck developments. Both procedures are leading to convection dominated convection diffusion equations. We study convergence for semi-discrete and fully discrete approximations of a such obtained equation, for a broad beam model, using the standard Galerkin and streamline diffusion finite element methods. The analytic broad beam model of the light ion absorbed dose were compared with the results of the modified Monte Carlo (MC) code SHIELDHIT+ and those of Galerkin streamline diffusion approach.

1. Introduction. We, primarily, assume a broad beam of forward-directed ions normally incident at the boundary of a semi infinite medium. As a result of collisions (due to forward-directness assumption), only, a very small portion of the ions is scattered to large angles. These are very few ions with a directional change beyond a certain minimal angle θ m , which have a diffusion-like transport behavior and an, almost, isotropic angular distribution, except at very low energies. The remaining significant portion of the ion particles, deflecting slightly (< θ m ) from their original direction, are convective and referred to as forward-directed ions. The bipartition model is based on a split, of the scattering integral (kernel), through adding and subtracting the diffusion ion source to the diffusion and straightforward equations, respectively, see [15]. A related approach, based on a split of the scattering crosssection into the hard and soft parts is given in [13].
In a previous study [2] we considered the bipartition model applied to solving three types of problems for ion transport in inhomogeneous media: (i) normally incident ion transport in a slab; (ii) obliquely incident ion transport in a semiinfinite medium; (iii) energy deposition of ions in a multilayer medium.
In the present paper we consider the case of, forward directed, ions injected into a medium with large atomic weight. The underlying partial differential equation is therefore the Boltzmann equation within the Fokker-Planck realm. More specifically, we shall study the finite element approximations for a broad beam model.
Neutral (photon, i.e. x-ray, and neutron) and charged (electron and ion) particle beams are extensively used in radiation therapy both for early cancer detection and dose computation algorithms see, e.g. [12] and [15].
An outline of this paper is as follows: In Section 2 a brief description of the ion transport model under continuously slowing down assumption (CSDA) and a derivation for the broad beam model based on Fokker-Planck development are given. In Section 3 both a semi-discrete and a fully discrete scheme are introduced, finite element approximations are applied for a convection-diffusion problem of a broad beam model, where the penetration variable is treated like a time variable. Further, stability and convergence of the semi-discrete problem are proved. Section 4 is devoted to derive optimal error estimates for characteristic Galerkin and streamline diffusion methods for two versions of the broad beam model. Finally in our concluding Section 5 we discuss some simulation results for the bipartition and Galerkin approaches and compare them with results from Monte Carlo simulations.
2. Ion transport models under CSDA. The ion transport describing the actual process of energetic ions in absorbing media is formulated as follows: Let f (x, v) denote the ion distribution function at the position x ∈ R 3 with velocity v ∈ R 3 , and set u := v/ |v| and E = |v| 2 , then f (x, u, E)du dE represents the ion distribution at point x ∈ R 3 , with direction between u and u + du and energy between E and E + dE. Due to the statistical balance principle we may write the following ion transport equation derived from the linear transport equation by Lewis and Miller in [14], where ρ denotes the stopping power and ω is the energy loss straggling factor, κ N ; is a constant depending on N the number of atoms and σ is the elastic cross section. We consider an ion beam of energy E 0 normally incident on the hypersurface of a semi-infinite medium. We let the outward normal to the semi-infinite region on the left to be along the positive x-axis inside the medium (we use the standard notation x = (x, y, z)). In the bipartition strategy the scattering integral is divided into two parts, of which one is the comparatively isotropic diffusion ion source, including almost all of the large-angle scattered ions, the other is the remaining part that spreads mainly in the forward, small-angle, direction. As physical model, we have considered a scattering kernel with strong algebraic fall-off behavior from its peaks at zero angle and zero energy. An example of such a kernel is an inverse power function approximation for the elastic cross-section, viz where C is a constant and is a positive integer which corresponds to the Jacobian of the algebraic fall-offs. See [15] and [16] for details. We shall also assume that 2.1. Ion transport in absorbing isotropic media. Lewis model (2.1) does not contain absorption term. Absorption has usually a regularizing effect and therefore tackling (2.1), with no absorption, is of more mathematical interest. However, in radiation therapy the absorption is of vital importance. Therefore, in what follows, we consider a modified version of (2.1) that yields an absorption term; σ a f below. Then, the wording inhomogeneous media is because σ a may depend on x. To this approach we relay on a Fokker-Planck development based on asymptotic expansions of the linear transport equation in isotropic media [16] (extended to anisotropic media in [5]). From the Fokker-Planck equation we can extract pencil beam or broad beam models. The pencil beam model is studied in an extensive amount of literature for the Fermi equation [11] see, e.g. [3]- [6], [13] and [16]- [17]. Here, we focus on broad beam model using the steady state Fokker-Planck equation where ε ≈ γ ≈ ∆, with ∆ = O(mean free path), are certain smallness parameters, In addition to introducing an absorption term, the justification in considering (2.3) is to circumvent the difficulties with small mean free path. More specifically, using asymptotic expansions the integral scattering operator in (2.1) is replaced by a differential Fokker-Planck operator in (2.3). The effect of this replacement is that the dominant in and out scattering terms cancel, thus effectively increasing the mean free path. Then, it is possible to consider numerical meshes of order O(∆). Physically, the smallness parameter γ is a measure of the peaking of the scattering kernel in angle, and can be roughly thought of as the deviation from the unity of the cosine of characteristic scattering angle. Likewise, the smallness parameter ε is a measure of the peaking of the scattering kernel in energy, and can be thought of as the characteristic value of the fractional energy change due to a single scattering.
For the broad beam model we have chosen the x-axis as the penetration direction. Thus, due to symmetry, the broad beam equation is independent of y, z and ϕ, i.e., (2.5) To proceed we consider (2.5) for forward-peaked beams (µ > 0) in a bounded domain D := {(x, µ, E) : 0 < x < L, 0 < µ < 1, 0 < E < E 0 }, associated with physically relevant boundary conditions. Recalling (2.4), our model problem is now which is of degenerate type with convection in (x, E) and diffusion in (µ, E). Note also that the term Q is actually supplied as the incident source at the boundary and is equivalent to a boundary condition. If we disregard forward-peakedness; µ ∈ [−1, 1] yields a Forward-backward problem (Kolmogorov equation, see [18]). Since collisions often result nonsmooth changes in the energy and direction of the particles on each collision site, therefore the discontinuous Galerkin method, in all three variables x, µ and E, is the most relevant finite element approach for the numerical study of this problem. However, the broad beam equation models very small changes in angle and energy, whence the discontinuities in these variables are not considered to be severe. Hence, in small neighborhoods of the jumps in energy and direction of the beam, the distribution function is approximated by continuous spline-type functions preserving beam configurations with the same inand outward profiles as the original ones. The discontinuity in position variable x, as dose delivery site, has often a more sensitive nature and therefore retained. Now we introduce, the inflow (outflow) boundary and n(x ⊥ ) is the outward unit normal to the boundary Γ at x ⊥ ∈ Γ. The final form of the broad beam model will then be a boundary value problem where the, homogeneous version of, equation (2.6) is associated with the boundary conditions, (2.8) In what follows, due to the nature of Dirac δ function we will be forced to such replacements, otherwise the energy estimates (in L 2 -type norms) will deteriorate. The differential equation in (2.8) is a degenerate type, convection dominated, convection-diffusion equation. Existence, uniqueness and regularity properties for this equation are as for general form of degenerate type equation given in [7].
3.1. The semi-discrete problem. We introduce a finite dimensional function space V h ⊂H 1 (I ⊥ ) with h being the maximal diameter in triangulation of I ⊥ and where for positive integer s, H s and · s denote the L 2 -based Sobolev space and norm, see [1]. An example of such V h is the set of sufficiently smooth piecewise polynomials P (x ⊥ ) of degree ≤ r, satisfying the boundary conditions in (2.8).
To proceed we introduce a bilinear form, A :H 1 (I ⊥ ) ×H 1 (I ⊥ ), defined by Then the continuous variational problem is: find a solution f to (2.8) such that is a smooth approximation for the product of the above two Dirac δ functions. Let nowf ∈ V h be an auxiliary interpolant of the solution f for (2.8) defined by Our objective is to solve the discrete variational problem: whereδ h is assumed to be a finite element approximation ofδ which coincides with the interpolantf (0, To distinguish, we use the inner product notations: (·, ·) ⊥ and (·, ·) Ω , where Ω = [0, L] × I ⊥ := I x × I ⊥ , for integrations over I ⊥ and I x × I ⊥ , respectively. Finally, we assume that the mesh size h is related to the scaling parameter ε by h 2 ≤ ε ≤ h, where in this section we shall use ε ∼ h.

3.2.
Stability. In this part we prove a stability lemma that guarantees the control of both continuous and discrete solutions by the data in the following triple norm Integrating over x ∈ (0,x),x ≤ L yields Consequently ∀x ≤ L, we have that which gives the first assertion of the lemma. Note that the two integrals with positive sign in (3.10) add up to 1 Transferring the negative integrals to the right hand side we obtain the second assertion of the lemma and the proof is complete.
From the above proof we can deduce a control of a quantity involving the norm .
Roughly speaking, (3.11) yields control of a quantity corresponding to the contribution of the ||| · |||-norm at eachx ∈ (0, L). Lemma 3.1 is stated and proved for the continuous problem. Using the same argument we can obtain the semidiscrete version of the stability estimate for standard Galerkin (SG) problem: For convenience, in the sequel we shall use the following boundary integral notation: where α = β or α =β (then n :=ñ), will be obvious from the content.

3.3.
Convergence. In this part we state and prove convergence rates, both in the L 2 -norm and in the triple norm, for the SG-method for the semidiscrete problem with weakly imposed boundary conditions. For the hyperbolic problems with an absorption term of O(1), and f ∈ H r (Ω), the optimal convergence rate for the SG method in the Proof. Subtracting first equations in (3.4) and (3.6) we get, using (3.3), and recalling thatf is an interpolant of f , where in the last step we used (3.5). This can be written in an equivalent form: Let χ = f h −f , then by the same argument as in the stability estimate we have, or equivalently, .
, and the definition of the ||| · |||β norm we have . (3.16) Now the desired result follows from the following interpolation estimate: Proof. Let f ∈ H r (Ω), then by interpolation estimates, see [8], there exists an interpolantf ∈ V h , of f and interpolation constants C 1 and C 2 such that Now recalling the definition of ||| · |||β we have, with · := · L2(Ω) , that , and the proof is complete.
Note that the C 3 -term is the worst one. Otherwise we could get a convergence of order O(h r−1/2 ). This result contains a σ a (E)-weighted L 2 norm estimate as: (3.4) and with f h being the solution of (3.6), there is a constant C = C(Ω, g) such that Proof. Recalling the definition of the triple norm we have that Thus, we obtain immediately (3.21).
Observe that in C = C(Ω, g), the Ω dependence is because of the E depending σ a (E), σ a may depend on x as well (then the proofs require minor modifications), while the g dependence comes from the assumed identity u h (0) :=ũ(0) = g.

3.4.
The fully discrete problem. To extend the semidiscrete algorithm for I ⊥ to I ⊥ × I x and include also discretizations in x, we consider the time discretization schemes as DG, BE and CN, for I x . To this end, we introduce the bilinear forms: where, for f and χ satisfying the boundary conditions of (2.8), both integral terms in (3.23) and (3.24) will vanish. Thus we may rewrite the problem (3.4) as finding a solution f ∈H 1 (I ⊥ ) such that, We subsequently use the finite dimensional subspace V h ⊂H 1 (I ⊥ ) and insert for the discrete solution f h a representation by separation of variables viz: 4. Streamline diffusion method. Since for a forward directed beam min(1 − µ 2 , ε)/ << 1, both the broad beam and Fokker-Planck models are, convectiondominated, convection diffusion equations. To obtain approximate solutions for these types of equations, we may use a certain type of Galerkin method called the Streamline diffusion (Sd) method described below. Because of a lack (weakness) of stability, the SG approximation contains oscillations not present in the exact solution of convection dominated problems. Whence, we need to improve the stability properties of the SG method, without sacrificing accuracy. We consider two ways of enhancing the stability of SG. (a) introducing a weighted least square term; (b) introducing an artificial viscosity. We refer to the Galerkin method with these modifications as the streamline diffusion method. Both modifications enhance stability without a strong effect on the accuracy. We begin by describing the Sd-method for an abstract linear problem for which SG method reads: find F ∈ V h such that where L is a linear operator on a vector space V and V h is a finite dimensional subspace of V . In our problem, L is the convection-absorption operator: The least squares method for (4.1) is to find F ∈ V h that minimizes the residual over V h in an appropriate norm, e.g. the usual L 2 norm, This is a convex minimization problem and the solution F is characterized by We now formulate a Galerkin least squares finite element method for (4.2) by taking a weighted combination of (4.2) and (4.5): compute F ∈ V h such that where δ is a parameter to be chosen. Adding an artificial viscosity modification yields the Sd-method in abstract form: find F ∈ V h such that whereε-term is an artificial viscosity defined in terms of the residual R(F ) = LF −g. We return to equation (2.8) and perform the differentiation with respect to µ, we get a degenerate equation, with convection in the direction of β c := (µ, 2µ, −1) and diffusion in x ⊥ . Unlike (4.3)-(4.6), now we have already a diffusion present in the equation (4.8). We shall consider a less viscous, but symmetric (with the same order of diffusion coefficients in µ and E) problem, where nowβ := (2µ, −1) and we have assumedε = min(1 − µ 2 , ε). Note thatε does not add any artificial viscosity to (4.8). Note also thatβ is different from the one in Section 3. Here, it is the δ-term in the streamline diffusion modification that will automatically add an extra amount of diffusion. Interpreting x > 0 as a time variable both convection and diffusion are in x ⊥ = (µ, E) and hence we have a, non-degenerate, "time dependent", convection-diffusion equation. Problem (4.9) can be stated in unbounded domain with a compact support for f : (4.10) This is most natural initial-boundary value problem for a forward directed charged particle beam, simply because there are no particles with infinite energy. In the radio therapy applications, obviously the energy range is finite and we need to state the problem in a bounded domain, where the second derivative with respect to µ in (4.9) requires imposing an additional boundary condition at µ = 1, where we recall that I × I ⊥ := (0, L) × (0, 1) × (0, E 0 ) and We assume x i ∈ I x , i = 0, 1, . . . , N , to be distinct nodes partitioning I x into subintervals I n = (x n−1 , x n ), n = 1, 2, . . . , N . We divide each space "time" slabs S n = I n × I ⊥ into prisms I n × τ , where τ is an element in a triangulation of I ⊥ .
Our study concerns a Sd method based on using approximations consisting of continuous piecewise linear functions in x ⊥ and discontinuous piecewise constants in x, denoted by cG(1)dG(0). We define the trial space W k , with k being a common mesh parameter, to be the set of functions w(x, x ⊥ ) defined on I x ×I ⊥ such that the restriction w| Sn is continuous and piecewise linear in x ⊥ and constant in x. Thus, on each slab S n , we seek approximate solution w for (4.11), such that w| S n ∈ W n k , W n k is the streamline diffusion space and V n is the space of continuous piecewise linear functions vanishing on ∂I ⊥ . We denote by and h the mesh parameters in I x and I ⊥ , respectively, we choose k = max( , h/ β ) and use the usual notation w(x n + s). (4.14) In this setting, the method reads as: compute F ∈ W k such that for n = 1, . . . , N , For compatibility reasons, e.g.
The error estimate procedure and the numerical scheme are simplified using the transformation between Euler and Lagrange coordinates: (x, x ⊥ ) = (x,x ⊥ +xβ). Then forf (x,x ⊥ ) = f (x, x ⊥ ), by the chain rule and assuming µ = 0, we get Below we summarize a priori and a posteriori error estimates for the cG(1)dG(0) scheme derived for the local Lagrangian coordinates (x,x ⊥ ), in the presence of an artificial viscosity produced by replacingε, in (4.9)-(4.11), byε = max(1 − µ 2 , ε).
where [∂ S w] denotes the jump across the side S ⊂ ∂τ in the normal derivative of the function w ∈ V h and C i is an interpolation constant. Finally, the star indicates that the corresponding term is present only if V n−1 is not a subset of V n .
The convergence rates in this theorem hold in the Euler coordinates (x, x ⊥ ), provided that there exists an affine bijection G n : S n → S n defined by i.e., G n takes a non-oriented grid in (x,x ⊥ ) to an oriented one in (x, x ⊥ ). Let nowβ h n ∈ V n denote the nodal interpolant ofβ n =β(x n−1 , ·) and set ( where∇ denotes the Jacobian with respect tox ⊥ : where I is the identity operator. The proof of the above theorem is a lengthy modification of the a priori and a posteriori error estimates for heat equation based on a dual problem approach derived in [10]. The error estimates in Euler coordinates are more involved. Nevertheless, for (4.11) with diffusion coefficientε, a proof based on the Lagrange coordinate framework as outlined through (4.17)-(4.24) is rather technical and seemingly cumbersome. Below we shall focus on an approach considering, more closely, the study of a bilinear form of type (4.15).

Remark 2.
Here, we tackle the problem with its own small amount of diffusion (even lessε = min(1 − µ 2 , ε), but symmetric in µ and E), and without add of artificial diffusionε-term in (4.7), or as ofε-term in the case of Theorem 4.1.
The procedure (4.17)-(4.24) could be more simple if instead ofβ =β/µ in (4.17), we could divide the equation (4.9) by a µ. This, however, is inadequate in the presence of the diffusive term (1 − µ 2 )f µµ in (4.8) and the subsequent equations. We could circumvent this obstacle by using a weak form of (2.8). Such details are not in the scope of this study and therefore we skip them.

Streamline diffusion with discontinuity in x.
To study the distribution of the particle beams in a certain depth, e.g., x d , a reasonable initial guess would be obtained using the information in some previous distinct depths x = x i , i = 1, ..., n, x i < x i+1 . This corresponds considering discontinuities in x-direction. Now, for each n, we let W n h to be a finite element subspace of H 1 (S n ) based on the triangulation C h of the strip S n with the elements of size h > ε. (Because of the relation between the parameters h and n, in this section, we may only include the index h in the discrete function spaces and hence in our discrete functions). Let now T h be a triangulation of I ⊥ , with elements τ ∈ T h , and, for each n, define Here H 1 n are H 1 (S n ) functions satisfying the same boundary conditions as in (4.11), restricted to S n . If we now apply the streamline diffusion method successively on each strip S n for the broad beam problem (4.11) and impose the boundary conditions at the points x = x n−1 weakly, we obtain the following method: Find f h ∈ W n h such that for n = 1, 2, ..., N , , g + n−1 n = (Q, g + δB(g)) n , ∀g ∈ W n h , where we have used the notation: Then, the continuous bilinear form for this method reads as and in the discrete counterpart f is replaced by f h . The continuous linear form is Thus, we have We shall use a stability estimate for this Sd method in a norm ||| · ||| defined by , Lemma 4.2. Let g + = 0 in (4.11) and assume that σ > (1 − δ(2 + α)) for some α > 0, then Proof. We use the above definition for B and write Recalling the definition ofB(w) we have that where, since σ and δ are chosen to be independent of x, we have that Thus, rearranging the terms we may write (4.28) Further, using the equality: Hence, combining with the coefficients in the first term in B(w, w) above, we get Assuming (4.31), which is the same as in the assertion of lemma, we have that (4.32) Further to estimate the term involving ∆ ⊥ w we use the inverse estimate to obtain N n=1 Now choosingε ∼ h 3/2 (ε ∼ h for convex Ω) and δ ∼ h, for sufficiently small h, Finally, the ∂ n -term corresponds to an outflow at µ = 1 which vanishes by assumption that g + ≡ 0 in (4.11), and − w + , w + Γ − n ≥ 0. Summing up the estimates (4.26)-(4.34) and summing over n yield the desired result.
Now since ξ ∈ W h , by the Galerkin orthogonality B(e, ξ) = 0 and we can write where we used the above lemma. Using Green's formula yields (4.36) x = L} is the outflow boundary except the top surface for x = L. Thus using the same technique as in the proof of the above lemma, we may write an estimate of the form By a standard inverse estimate and recalling that δ ∼ h andε ∼ h 3/2 we have Similarly the following estimates hold true Combining with the interpolation estimate, see Ciarlet [9], and since |||η|||, the interpolation error, is of the same order as |||ξ|||, we get the desired result.
Remark 3. The main advantage of the Sd method is that: we have improved stability; and the rate of convergence for the SG method, in Section 3, both in the ||| · ||| and the L 2 -norm, by a factor of O( √ h), see the estimates (3.15), (3.21) and (3.22) which all are of order O(h r−1 ). Note that here k corresponds to r − 1 of Section 3. Moreover, using the Sd method we have been able to control the jump discontinuities in x. In the presence of such discontinuities, the proofs in Section 3 (Lemma 3.4, Proposition 3.5 and Theorem 3.6) are not extendable to I x × I ⊥ .
Below we focus on some numerical results comparing the Sd, bipartition and Monte-Carlo (MC) methods. To keep the computational costs in a realistic level, our simulations concern the energy ranges ≤ 30M eV , adequate for light ions and high-energy electron particles (for high energy ranges of ≈ 200M eV see [2]). For clinical applications, the accuracy of the simulations in a moderate energy range of, say, ∼ 50 − 100M eV are more desirable than the simulations for very high energies.
5. Numerical simulations. Given g := g(µ, E), we define the boundary conditions for broad beam model as: Here g is a smooth approximation of the source term Q = 1 2π δ(x)δ(1 − µ)δ(E 0 − E). We have considered the broad beam problem (4.8) associated with the boundary condition (5.1) and derived numerical algorithms for SG and Sd finite element methods, where the Sd method coincides with a characteristic Galerkin (CG)-method. We shall present only the results of the Sd part. For the sake of completeness, we briefly comment the results of our observations testing other methods than the Sd. We have carried out implementations to illustrate the applicability of the algorithms using different types of initial data approximating the Dirac δ function (and consequently the source Q). To begin with, semi-streamline diffusion (SSD) and characteristic streamline diffusion (CSD) methods are more stable and accurate than the SG and CG methods for all the canonical forms of the initial data (approximating a Dirac δ function beam source by a Maxwellian, cone, cylinder and a hyperbolic beam). As for the convergence: solutions with modified Dirac initial data are suited in CS and SSD. Maxwellian initial conditions produce accurate results in the CSD scheme, whereas the hyperbolic initial conditions produce more accurate results in the SSD scheme. The oscillatory behavior, appeared considering non smooth data, can be eliminated by modifying the L 2 -projection. The formation of layers can be avoided taking small steps in the penetration variable. However, a better approach to deal with this phenomenon is through adaptive refinements. 5.1. Results. We use the Sd-method to solve (4.8) with the boundary condition (5.1) and integrate f (x, µ, E) over µ and E to get the energy and angular distributions for different depths. We compare the results of dose delivery of intensity 10, 20, 30 M eV electron in water, with those of the bipartition model in [15] and [19], and MC method in [12] and [19]. Here, we have neglected the contributions from the secondary particles. In Fig. 2, 5 and 8 we can see that our results are very close to MC, only the positions of the maximum values are different. The reason is that the stopping power that we have used in our Sd model is somewhat different from that of MC code SHIELD-HIT+, see [12] for details. We use the same stopping power with the bipartition model, and find out that the positions of the maximum values are very close as seen in Fig. 1, 4 and 7. Bipartition model employs CSDA and singles out particles with larger angle and energy variations. Then, the energy distributions are very narrow and the maximal drops exponentially fast. Similar phenomena appears for the angular distributions in Fig. 3 and 6.

Conclusion.
We consider a broad beam model derived from the Fokker-Planck equation obtained, from the linear transport equation, by asymptotic expansion. We approximate the broad beam equation by a variety of finite element methods, where we (i) derive stability and convergence estimates for a semi-discrete standard Galerkin method, (ii) formulate a fully discrete scheme, (iii) study the streamline diffusion method and give sharp error bounds in local Lagrangian coordinates. Under certain assumption on the characteristics and by the inverse mapping theorem, these bounds are valid in Euler coordinates as well. (iv) We perform differentiation in µ and derive a non-degenerate, convection-dominated equation, introduce a diffusion correction (drop) within theε ≈ min(1 − µ 2 , ε)-term, and derive, without adding artificial viscosity, optimal convergence rate of O(h k+1/2 ) for the exact solution f ∈ H k+1 (Ω). (v) In implementations, we use the streamline diffusion method to calculate the energy and angular distributions for the electrons and light ions and compare the results with those obtained by bipartition and Monte Carlo simulations. In our knowledge this approach is not considered elsewhere.
We plan to study radiation beams with inhomogeneous data in anisotropic media. In a forthcoming work we shall study high energy ions, including secondary particles and perform implementations based on clinical data.