Andrieu François
email: francois.andrieu@u-psud.fr

Douté Sylvain

Schmidt Frédéric

Schmitt Bernard

Radiative transfer model for contaminated rough slabs

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Hyperspectral imaging has become a major component in planetary surface observation since the past decades. Earth and other Solar system bodies are now observed in various spectral ranges at various resolutions and from various heights.

As a photon come across a surface, it interacts in two major ways. It can be either absorbed or deviated (scattering, diffraction, refraction). The objective of this radiative transfer model is to describe the interactions using a realistic surface description. In such descriptions, the reflectance of a surface is the result of multiple interactions, with multiple irregular interfaces of different materials. The exact resolution of the radiative transfer equations turns out to be a highly difficult and time consuming problem. This problem has been solved under certain hypothesis : if the characteristic size of the particle is much smaller than the wavelength, or if it is much bigger. In this study, we consider the geometrical optics domain, where the particles are much bigger than the wavelength. For example, in the visible and near infrared range (400 nm -5 µm), we suppose that the average particle size does not fall below 10 µm. This is in general valid for planetary surfaces [START_REF] Peltoniemi | Light scattering in planetary regoliths and cloudy atmospheres[END_REF]. Ray tracing algorithms [START_REF] Grynko | Scattering matrix calculated in geometric optics approximation for semitransparent particles faceted with various shapes[END_REF][START_REF] Chang | Ray tracing in absorbing media[END_REF][START_REF] Pilorget | A radiative transfer model to simulate light scattering in a compact granular medium using a monte carlo approach: Validation and first applications[END_REF][START_REF] Ben | Polarized radiative transfer in an arbitrary multilayer semitransparent medium[END_REF] that simulate the very complex paths of millions of photons through these surface can show very accurate results, but they depend on a huge number of parameters and are highly time consuming, se-riously limiting the the interpretation of extensive hyperspectral images. We aim at a radiative transfer model that is fast enough to be able to deal with a vast amount of data, such as planetary spectro-imaging databases. It is then necessary to make further simplifying assumptions that enable the formulation of approximate analytic or semi-analytic solutions to the radiative transfer problem. A possible simplification is to consider that the radiative properties inside a media can be described statistically only using local mean properties of scattering and absorption [START_REF] Kubelka | New contributions to the optics of intensely light-scattering materials. part i[END_REF][START_REF] Hapke | Bidirectional reflectance spectroscopy: 1. theory[END_REF][START_REF] Shkuratov | A model of spectral albedo of particulate surfaces: Implications for optical properties of the moon[END_REF]. The media is assumed to be homogenous at a mesoscopic scale. Another classical simplifying assumption considered in such problem is the two stream approximation [START_REF] Kubelka | New contributions to the optics of intensely light-scattering materials. part i[END_REF][START_REF] Hapke | Bidirectional reflectance spectroscopy: 1. theory[END_REF][START_REF] Kylling | A reliable and efficient twostream algorithm for spherical radiative transfer: Documentation of accuracy in realistic layered media[END_REF]. It has been shown that under certain conditions, it did not affect too much the solution compared to more accurate studies, but simplifies greatly the calculations [START_REF] Vargas | Applicability conditions of the kubelka munk theory[END_REF][START_REF] Vargas | Two-flux radiative transfer model under nonisotropic propagating diffuseradiation[END_REF]. To describe the reflectance of a surface, one also have to consider the geometry of illumination and observation. In our approach, these photometric effects are modeled by the properties of the interface between the media and the exterior. These properties of roughness can also be statistically described, using only one or a few parameters. Y. Shkuratov [START_REF] Shkuratov | A model of spectral albedo of particulate surfaces: Implications for optical properties of the moon[END_REF] and B. Hapke [START_REF] Hapke | Bidirectional reflectance spectroscopy: 1. theory[END_REF] developed analytical radiative transfer models for granular media, that are able to simulate the bidirectional reflectance of various granular surfaces.

If the media cannot be described as homogenous, it is possible to consider it pieciwise continuus, constituted of homogeneous strata. It is the case for example in the atmospheres, or stratified surfaces. A family of models describe the radiative transfer in stratified media, such as the DISORT algorithm [START_REF] Stamnes | Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media[END_REF]. In these discrete-ordinate modelisations, each layer is considered homogenous, and the total reflectance is calculated iteratively, by adding the contribution of each layer. This method has also been adapted to the study of the ocean-atmosphere coupled system with a rough surface [START_REF] Jin | Analytical solution of radiative transfer in the coupled atmosphere-ocean system with a rough surface[END_REF].

Starting from Hapke model, improving it, and combining it with a multi-layer method [START_REF] Stamnes | Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media[END_REF], S. Douté has also developed a model for stratified granular surfaces [START_REF] Douté | A multilayer bidirectional reflectance model for the analysis of planetary surface hyperspectral images at visible and near-infrared wavelengths[END_REF]. Using the same strategy, we developed a semi-analytical radiative transfer model for a compact layer (solid matrix containing inclusions) overlaying an optically thick granular layer. This two layers approach does not require an iterative DISORT-like method, but only adding coupling formulas. It is founded on three major assumptions : (i) the geometric optics conditions are observed, (ii) the media is piecewise continuous and (iii) the inclusions are close to spherical and homogeneously distributed in the matrix.

Model overview

We decompose the reflectance into two distinct contributions : specular and diffuse. We chose the Hapke [START_REF] Hapke | Bidirectional reflectance spectroscopy: 3. correction for macroscopic roughness[END_REF] probability density function of orientations, as it well describes the statistic distribution of slopes in the approximation of small angles. We consider a collimated incident radiation, at an incident angle i. We estimate the specular contribution, considering the geometry and the surface description. The specular reflection of rough surfaces have been studied in various cases [START_REF] Hapke | Bidirectional reflectance spectroscopy: 3. correction for macroscopic roughness[END_REF][START_REF] Cox | Measurement of the roughness of the sea surface from photographs of the sun's glitter[END_REF][START_REF] Muhleman | Symposium on Radar and Radiometric Observations of Venus during the 1962 Conjunction: Radar scattering from Venus and the Moon[END_REF][START_REF] Saunders | Shadowing on the ocean and the existence of the horizon[END_REF][START_REF] Van Ginneken | Diffuse and specular reflectance from rough surfaces[END_REF]. We use the same general idea of these methods, describing the rough surface as constituted of multiple unresolved facets. The specular contribution will result from the integration of the specular reflections on the facets, in the solid angles considered (i.e., the light source and the detector) as described in figure 1a.

Then we estimate the diffuse contribution. The total reflection coefficient at the first rough interface, that determines the amount of energy transmitted to the slab, is obtained by integrating specular contributions in every emergent direction, at a given incidence. We consider that the first transit through the slab is anisotropic (collimated), and that there is an isotropisation at the second rough interface (i.e. when the radiation reach the semi-infinite substrate). For the refraction and the internal reflection, every following transit is considered isotropic. The diffuse contribution is obtained using an analytical estimation of Fresnel coefficients [START_REF] Douté | A multilayer bidirectional reflectance model for the analysis of planetary surface hyperspectral images at visible and near-infrared wavelengths[END_REF][START_REF] Chandrasekhar | Radiative transfer[END_REF], and a simple statistical approach. The contribution of the semi infinite substrate is estimated using Hapke model [START_REF] Hapke | Theory of Reflectance and Emittance Spectroscopy[END_REF]. Finally, we consider that the slab is under a collimated radiation from the light source, and under a diffuse radiation from the granular substrate. We compute the resulting total bidirectional reflectance using adding doubling formulas (figure 1b).

SURFACE ROUGHNESS -FACETS DISTRIBUTION

The first step is to describe the roughness of the surface. We consider that it is composed of N facets that are not resolved, with N 1. These facets' orientations follow a probability density a (ϑ, ζ), where ϑ is the zenital angle between the normal to the facet and the local vertical direction, and ζ is the azimutal angle. To make our approach as general as possible, we chose to describe the surface as randomly rough. Such a roughness has already been widely studied (see for example [START_REF] Hapke | Bidirectional reflectance spectroscopy: 3. correction for macroscopic roughness[END_REF][START_REF] Cox | Measurement of the roughness of the sea surface from photographs of the sun's glitter[END_REF][START_REF] Muhleman | Symposium on Radar and Radiometric Observations of Venus during the 1962 Conjunction: Radar scattering from Venus and the Moon[END_REF][START_REF] Saunders | Shadowing on the ocean and the existence of the horizon[END_REF][START_REF] Van Ginneken | Diffuse and specular reflectance from rough surfaces[END_REF][START_REF] Lumme | Radiative transfer in the surfaces of atmosphereless bodies. i -theory. ii -interpretation of phase curves[END_REF] and the ... reference cited in these papers). These studies show that a slope distribution, with tan ϑ that is close to Gaussian is a good description of the surface. Such a description combines simplicity and efficiency reproducing the photometric variations. For the sake of simplicity and because it is widely used in the literature, we chose the following probability distribution function [START_REF] Hapke | Bidirectional reflectance spectroscopy: 3. correction for macroscopic roughness[END_REF] :

S i £ £ 1-S i 1-S i S i 1-S i R R spec spec R 0 ' T 0 ' F F (
a (ϑ, ζ) = 1 π 2 tan 2 θ exp - tan 2 ϑ π tan 2 θ sec 2 ϑ sin ϑ (1)
where

tan θ = 2 π π 2 0 a (ϑ) tan ϑ dϑ (2)
It is supposed that the azimutal distribution is uniform. The angle θ representing the mean slope angle completely characterizes the facets' orientations and the surface roughness. This slopes distribution considers the cases of small θ. Practically, the threshold of validity can be determined depending on the level of tolerance (see Figures 5 and6. The expression of a (ϑ, ζ) could be adapted in the future to extend the study to any type of terrain, as discussed in section A.

z y x I E W f # ³ i e à A f
Fig. 2. The local coordinate system ( x, y, z) is centered at the slab surface, with the z axis vertically upward and the x axis horizontally toward the sun. I and E are respectively the incident and the emergent directions. W f is the normal to the facet A f .

SPECULAR REFLECTANCE

In this part we determine the specular contribution to the reflectance at the detector. We first establish the relations between the orientation ϑ spec , ζ spec of a facet that is in specular conditions and the geometry of observation (i, e, ψ), where i is the in- cidence angle, e the emergence and ψ the azimuth (see Figure 2). The total specular contribution is obtained by integrating these relations, taking into account the geometry variations within one pixel and the statistics of the slopes a (ϑ, ζ).

A. Specular conditions for one facet

For one facet A f to satisfy the specular reflection conditions, then its normal W f must respect (see figure 2):

E = -I + 2 I| W f • W f (3)
where the operator | represents the scalar product. If we express Eq. 3 in the ( x, y, z) coordinates, :

       x e =
where x e = sin e cos ψ, y e = sin e sin ψ and z e = cos e are the coordinates of the emergent vector E in the ( x, y, z) frame. This leads to :

        

B. Expression of the specular reflectance for one pixel

We consider a pixel of area A formed of N facets of same area A f , orientated according the probability density a (ϑ, ζ), as de- tailed in section 2, with N 1. The number of facets satisfying the specular reflection conditions defined in section A will be H C N a ϑ spec , ζ spec d(ϑ, ζ), where H C is the set of values (ϑ, ζ) satisfying 5 within the range of observation geometries. Indeed, there is a range of different geometry of observation within one instruments' pixel. Let χ c be the emergence variations within a pixel. The facet orientation satisfying the specular conditions for this geometry is ϑ spec , ζ spec given by Eq. 5. At the incidence i, there is a range of orientations that satisfy specular condition within a pixel that is centered at ϑ spec , ζ spec and has the size δ (ϑ, ζ). δ (ϑ, ζ) is determined using a function

g i :   e ψ   →   g 1 (e, ψ) = ϑ g 2 (e, ψ) = ζ  
that transforms (e, ψ) into (ϑ, ζ), for the incident angle i. Not every facet that satisfy these conditions will send energy to the captor. Indeed, the roughness of the surface introduces a shadowing of the scene : some facets will not receive incident light, or will not be visible by the captor, or both. A shadowing factor S must be introduced at this point. Let N spec be the number of facets that satisfy both the geometrical condition defined in section A and the visibility condition :

N spec = H C N a ϑ spec , ζ spec S i, e, ψ, θ d(ϑ, ζ) (6)
with S i, e, ψ, θ a shadowing factor that depends on the geometry of observation and the roughness of the surface [START_REF] Hapke | Bidirectional reflectance spectroscopy: 3. correction for macroscopic roughness[END_REF]. Each one of these N spec facets receives an incident power P i and send back a reflected power P r :

P i = FA f cos α 2 (7)
P r = P i r f α 2 (8)
F being the incident power flux (e.g. : Solar flux) in the radiation direction, A f cos α 2 the projection of the facet in the plane orthogonal to the incident radiation, and r f α 2 the Fresnel reflection coefficient in energy at the phase angle α , r f = R 2 ⊥ (α) + R 2 (α). As α does not depend on the facet ori- entations, all these specular reflections will result in a specular power P spec = N spec P r , thus :

P spec i, e, ψ, θ = H C NFA f cos α 2 r f α 2 a ϑ spec , ζ spec S i, e, ψ, θ d(ϑ, ζ) (9)
The reflectance factor R is the ratio between the bidirectional reflectance r of the surface and the one r L of a perfectly lambertian surface , thus R = π r cos i . The bidirectional reflectance r is the ratio between the radiance L of the surface and the collimated incident power, perpendicularly to the incident direction. Thus r = L F , with L = P ΩA cos e , A being the illuminated surface and Ω c the solid angle subtended by a pixel. Finally,

R spec = π r spec cos i = π L spec F cos i =π P spec Ω c AF cos i cos e , thus :
R spec i, e, ψ, θ = H C π N A f cos α 2 Ω c A cos i cos e S i, e, ψ, θ r f α 2 a ϑ spec , ζ spec d(ϑ, ζ) (10)
where Ω c is the solid angle subtended by an instrument's pixel.

A is the sum of the horizontal projection of all the facets : A = N A f cos ϑ . The term cos ϑ is included in the shadowing function S i, e, ψ, θ described by B. Hapke [START_REF] Hapke | Bidirectional reflectance spectroscopy: 3. correction for macroscopic roughness[END_REF]. Thus we can simplify Eq. 10 as

R spec i, e, ψ, θ = H C π cos α 2 Ω C cos i cos e S i, e, ψ, θ r f α 2 a ϑ spec , ζ spec d(ϑ, ζ) (11) d(ϑ, ζ
) is derived from the integration angles e and ψ that are the emergence and azimuth angles. There is a bijection between Ω C and H C because 5 admits a unique solution for every (e, ψ). Considering that the incidence angle i is a constant, we can rigorously express R i spec as :

R i spec (i, e, ψ) = π Ω c Ω c cos α 2 cos i cos e S i, e, ψ, θ r f α 2 a ϑ spec , ζ spec det J g i (e, ψ) de dψ (12)
det J g e (e, ψ) is the Jacobian of the function g i . This expres- sion 12 assumes that the incidence i is a constant. In reality, the light source is almost never completely collimated, but ranges inside a solid angle (e.g. : the Solar disk). Let Ω S be the solid angle of the source. Then the total specular contribution within one pixel will be :

R spec (i, e, ψ) = 1 Ω S Ω S R i spec (i, e, ψ) det J g e (i, ψ) sin i di dψ (13
) where det J g e (i, ψ) is the Jacobian of the function g e :

  i ψ   →   g 1 (i, ψ) = ϑ g 2 (i, ψ) = ζ   that transforms (i, ψ) into (ϑ, ζ)
, for the given emergence angle e.

DIFFUSE REFLECTANCE

We consider a two layers model, with a slab overlaying a semiinfinite granular substrate. The collimated radiation from the sun is transmitted to the slab with an external reflection coefficient S e (the prime here represent the anisotropy). We suppose an isotropisation at the second interface. The slab is modeled as a compact isotropic and homogeneous matrix. It contains inclusions that are close to spherical and not identical to the matrix. The inclusions are the main contributors to the scattering of radiation in the layer. They are distributed homogeneously in the matrix. The determination of the Fresnel coefficients at the interface matrix/inclusion or inclusion/matrix is a key to estimate the transmission and reflection factors of the layer. An internal and external reflection coefficient S ik and S ek for each type of inclusion k must be defined.

In this part we describe the radiative transfer in the media. First we will characterize the transmission of light into the slab. By energy conservation this is equivalent to calculating the total reflected power which, normalized by the incident energy, stands for the reflexion coefficient (section A). Then we will describe the scattering of light by the inclusion during the transfer through the slab. This requires the calculation of the external and internal reflection coefficients of these inclusions (section B). Once the basic optical properties of the inclusions are known, we can consider fluxes of energy within the whole slab that will be governed by the radiative properties of the slab (section E). Solving this radiative transfer problem within the slab with an upper and lower optical interfaces will give the overall reflection and transmission factors of the slab (section F). Finally the radiative interaction of the two layers (substrate and slab) are considered and solved by adding doubling leading to the final result (section G).

A. Reflection coefficients for the slab

Anisotropic case

Let S e be the external reflection coefficient in a collimated case (interface atmosphere/ice matrix). It corresponds to the ratio between the incident power P i and the total reflected power, in every direction P tot r . The total reflected power can be estimated integrating the specular contributions for every emerging direction, at the given incidence angle i.

S e = 2π dP spec AF cos i ( 14
)
dP spec being the specular contribution in a given geometry. Using Eq. 9, the expression of S e becomes

S e = H cos α 2 r f α 2 a ϑ spec , ζ spec S i, e, ψ, θ cos i d(ϑ, ζ) ( 15
)
where H is the set of values taken by ϑ and ζ throughout the integration. Exactly like in section B, d(ϑ, ζ) is derived from the integration angles e and ψ that are the emergence and azimuth angles. There is now a bijection between B and H, B being the superior hemisphere that is the domain of variation of e and ψ. Considering that the incidence angle i is a constant, we can express S e as :

S e = π 2 0 2π 0 cos α 2 r f α 2 a ϑ spec , ζ spec S i, e, ψ, θ cos i × det J g i (e, ψ) de dψ (16)
where det J g (e, ψ) is the Jacobian of the function

g i :   e ψ   →   g 1 (e, ψ) = ϑ g 2 (e, ψ) = ζ  
that transforms (e, ψ) into (ϑ, ζ), for the incident angle i. The internal reflection coefficient S i in a collimated case at the interface ice matrix/atmosphere is not considered as we suppose an isotropisation of the radiation at the second interface (ice/granular regolith).

Isotropic case

In the isotropic case, the internal reflection coefficient S i is obtained integrating the Fresnel equations at the surface for the all geometries:

S i = π 2 0 R 2 ⊥ (α) + R 2 (α) cos α dα (17)
where R ⊥ (α) and R (α) are Fresnel reflectivites for perpendic- ular and parallel polarization with respects to the propagation plan, for an incidence angle α and will be detailed later.

The external reflection coefficient S e is estimated the same way :

S e = π 2 0 R 2 ⊥ (α) + R 2 (α) cos α dα (18)

B. Reflection coefficients for the inclusions

In the case of a spherical inclusion of the type k, the internal reflection coefficient S ik is obtained in the usual fashion integrating the Fresnel equations (see [START_REF] Hapke | Theory of Reflectance and Emittance Spectroscopy[END_REF] Hapke, 2012, sect. 5.4.4, pp.78-95)

S ik = π 2 0 R 2 ⊥ (α) + R 2 (α) cos α sin α dα (19)
For the estimation of the external reflection coefficient S ek , a differential absorption factor is taken into account. Indeed, as we deal with inclusions in an absorbing matrix, the parallel rays we consider in the integration touch the inclusion after different optic paths. For a ray that touches the inclusion with an incidence α, the differential path length in the matrix is ν = ρ k cos α, where ρ k is the radius of the spherical inclusion. Thus the differential absorption factor is e -a m ρ k (1-cos α) , where a m is the absorption coefficient of the matrix. Writing the matrix's optical index n m + i k m , the dispersion relation gives

a m = 4π λ k m .
Finally, the external reflection coefficient S ek at the interface matrix/inclusion is :

S ek = π 2 0 R 2 ⊥ (α) + R 2 (α) e -a m ρ(1-cos α) cos α sin α dα (20)
This differential absorption effect is already taken into account in the expression of R ⊥ (α) and R (α) in the case of the internal reflection at the interface inclusion/matrix.

C. Fresnel coefficients

The Fresnel reflectivities for perpendicular and parallel polarization with respects to the propagation plane, for an incidence angle α, R ⊥ (α) and R (α), are derived from Snell's law (see [21] Hapke, 2012, sect. 4.3, pp.46-60) :

R ⊥ (α) = (cos α -G 1 ) 2 + G 2 2 (cos α + G 1 ) 2 + G 2 2 (21) R (α) = n 2 -k 2 cos α -G 1 2 + [2nk cos α -G 2 ] 2 [(n 2 -k 2 ) cos α + G 1 ] 2 + [2nk cos α + G 2 ] 2 (22)
using n = n 1 n 2 +k 1 k 2 n 2 1 +k 2 1 and k = n 1 k 2 -n 2 k 1 n 2 1 +k 2 1
, with n 1 + ik 1 and n 2 + ik 2 the complex refractive indexes of the media considered.

G 2 1 = 1 2 n 2 -k 2 -sin 2 α + n 2 -k 2 -sin 2 α 2 + 4n 2 k 2 1 2 (23)
G 2 2 = 1 2 -n 2 -k 2 -sin 2 α + n 2 -k 2 -sin 2 α 2 + 4n 2 k 2 1 2 (24)
D. Integration S. Chandrasekhar showed (see Chandrasekhar, 1960, sect. 22, pp. 61-69 [START_REF] Chandrasekhar | Radiative transfer[END_REF]) that many radiative transfer integration can be approximated using the Gauss quadrature formulae. If f (µ) is a polynomial of order 2m -1, then

1 -1 f (µ) dµ = m ∑ j=1 c j f µ j (25)
where µ 1 , ..., µ m are the zeros of the Legendre polynomials P 1 , ..., P m of order 1, ..., m, and c 1 , ..., c j are the associated Christoffel numbers :

c j = 1 P m µ j 1 -1 P m (µ) µ -µ j dµ ( 26
)
Equation 25 is exact if f (µ) is a polynomial of order 2m -1.
When f (µ) is not a polynomial, then the quadrature formulae gives an approximation that converges to the exact value when m → ∞. The order m of the approximation directly governs its quality. We estimate analytically the internal and external reflection coefficients in the isotropic case S i an S e using the roots of the 32 th order Legendre's polynomial and the associated Christoffel's numbers as detailed in [START_REF] Davis | Abscissas and weights for gaussian quadratures of high order[END_REF]. We use a simple change of variable to transform the integration interval from [0, π 2 ] into [-1, 1]. All the integrations are performed using the Gauss quadrature formulae, except S e and R spec . In these cases, the integration being a double one, we cannot use the Gauss quadrature. We chose after numerical tests an adaptive grid and the rectangle method.

E. Radiative properties of a slab containing inclusions

We suppose an homogeneous distribution of isotropic inclusions inside the slab. The inclusions type is noted by k, with N i different types, defined by different geometrical and optical properties.

E.1. Proportions of inclusions

We define the slab compactness γ c as the volume of the ice matrix per unit of volume. We also define N the total number of inclusions per unit of volume, and N k the number of inclusions of the type k per unit of volume. The proportion of each type of inclusion is P k = N k N . Immediate geometrical calculations give :

N = 3 (γ c -1) 4π ∑ N i k=1 P k ρ k (27)

E.2. Cross sections

We suppose close to spherical inclusions. The scattering efficiency for a sphere has been described by B. Hapke ([21] Hapke, 2012, sect. 5.6, pp.95-99, eq 5.52a) in his equivalent slab model. For an inclusion of the type k :

Q sk = S ek + (1 -S ek ) (1 -S ik ) 1 -S ik Θ ik Θ ik ( 28
)
where S ik and S ek are respectively the internal and external reflection coefficients of an inclusion expressed in equations 19 and 20, and Θ ik is the internal transmission coefficient of and inclusion. In the two stream approximation, and assuming the isotropy of the phase function of the internal scatterers in an inclusion ( [START_REF] Hapke | Theory of Reflectance and Emittance Spectroscopy[END_REF] Hapke, 2012, sect. 6.5, pp.122-144, Eq. 6.26) the expression of Θ ik can be reduced simply to :

Θ ik = r ik + exp -ρ k a ik (a ik + s ik ) 1 + r ik exp -ρ k a ik (a ik + s ik ) (29)
a ik being a type k inclusion's absorption coefficient, s ik its scattering coefficient, and

r ik = 1 - a ik a ik +s ik 1 + a ik a ik +s ik (30)
The scattering cross section σ sk for one inclusion is

σ sk = σ k Q sk ( 31
)
where σ k is the geometrical cross section : σ k = πρ 2 k . Let σ s be the mean cross section of the inclusions :

σ s = N i ∑ k=1 P k σ sk (32)
We do the approximation of geometric optics, so the extinction cross section σ ek correspond to the geometrical cross section σ k .

E.3. Single scattering albedo and optical thickness

The single scattering albedo ω of an absorbing and scattering object is defined as the ratio of the total amount of power scattered to the total amount of power removed to the wave (absorbed or scattered). We propose a simple statistical approach to express the single scattering albedo of a unit of volume of slab containing inclusions. We use the same method as [START_REF] Hapke | Theory of Reflectance and Emittance Spectroscopy[END_REF] Hapke, 2012, sect. 7.4, pp.158-169, but we modify the medium description. After a travel of length dν, the probability p 1 for a photon to meet an inclusion and be scattered is:

p 1 = 1 -exp -N σ s ln γ c γ c -1 dν (33)
The probability p 2 that this photon has not been absorbed by the matrix before is:

p 2 = exp (-a m dν) (34)
Thus the probability p s for a photon to be only scattered per unit of length is:

p s = 1 dν exp (-a m dν) 1 -exp -N σ s ln γ c γ c -1 dν (35)
that becomes for an infinitesimal length dν:

p s = N σ s ln γ c γ c -1 + • (1) (36)
equally, the probability p 3 for a photon to be absorbed or scattered by an inclusion throughout dν is:

p 3 = exp (-a m dν) 1 -exp -N σ e ln γ c γ c -1 dν (37)
and the probability p 4 to be absorbed by the matrix during dν is:

p 4 = 1 -exp (-a m dν) (38)
so the probability of extinction p e per unit of length is:

p e = 1 dν 1 -exp -N σ e ln γ c γ c -1 + a m dν (39)
when dν is close to 0, it becomes:

p e = N σ e ln γ c γ c -1 + a m + • (1) (40)
Finally we obtain the single scattering albedo of a slab containing inclusions dividing p s by p e :

ω = N σ s N σ e + γ c -1 ln γ c a m (41)
Equation 39 gives the expression of the optical depth τ of a slab with inclusion :

τ = N σ e ln γ c γ c -1 + a m ν (42)

F. Diffuse reflectance and transmission factors of the contaminated slab

F.1. Diffuse reflectance of a slab under collimated illumination

In this section, we suppose that the slab is under a collimated radiation. As in section 2, we suppose that the surface is constituted of N unresolved facets that have a slope distribution given by the probability density function a (ϑ, ζ). Each facet will receive an illumination at an incidence i f depending on its orientation. We consider in that case that the first transit in the slab is collimated and will transmit rays of light into the slab at different inclinations. Our goal at this point is to determine the mean transmission path D trough a slab of a given roughness θ and a given thickness D. For a facet with the orientation defined by (ϑ, ζ) using Snell-Descartes law gives D = fact(ϑ, ζ)D, with:

fact (ϑ, ζ) = 1 -1 nm cos i + cos ϑ 1 nm cos i f -1 -1 n 2 m 1 -cos 2 i f (43)
i f being the incidence angle on the facet (i.e. the angle between the facet's normal and the incident radiation). Using basic trigonometric relations gives [START_REF] Hapke | Theory of Reflectance and Emittance Spectroscopy[END_REF]:

cos i f = sin i sin ϑ cos ζ + cos ϑ cos i (44)
We consider that only the first transit in the slab is anisotropic. The internal absorption factor for the first anisotropic transit Θ will depend on the mean length L of this transit, with D =factD , and

fact = 2π 0 π 2 0 fact (ϑ, ζ) a (ϑ, ζ) dϑ dζ ( 45
) thus fact = 1 π 2 tan θ × 2π 0 π 2 0 e -tan 2 ϑ π tan 2 θ sec 2 ϑ sin ϑ -1 nm cos i + cos ϑ 1 nm cos i f -1 -1 n 2 m 1 -cos 2 i f dϑ dζ (46)
and

Θ = r m + exp -D a m (a m + p s ) 1 + r m exp -D a m (a m + p s ) (47)
where p s is given by Eq. 36, r m = 1-

√ 1-ω 1+ √ 1-ω
, and ω is given by Eq. 41.

The internal absorption factor for an isotropic transit is ([21], Hapke, 2012 Eq. 6.26)

Θ = r m + exp -2D a m (a m + s) 1 + r m exp -2D a m (a m + s) (48)
Every following transit is considered isotropic. As illustrated on figure 1, we can express the reflectance of the slab under a collimated radiation R 0 as

R 0 = S e + 1 -S e Θ S i Θ (1 -S i ) 1 + ∞ ∑ n=1 (ΘS i ) 2 (49)
S i being the internal reflection coefficient of the slab. The term S e represents the integration over the sky of the specular reflectance, and the other represents the diffuse reflectance. Thus we can express the diffuse reflectance of the slab as

R 0 = (1 -S e ) Θ S i Θ (1 -S i ) 1 -(ΘS i ) 2 (50)
The diffuse transmission of the slab under a collimated radiation T 0 is obtained the same way :

T 0 = Θ (1 -S e ) (1 -S i ) 1 -(ΘS i ) 2 (51)

F.2. Diffuse reflectance of a slab under isotropic illumination

In this section we suppose that the slab is under an isotropic radiation. Indeed, at the lower interface, it is illuminated isotropically from below by the substratum. R 0 and T 0 have their usual expressions in this case :

R 0 = S e + (1 -S e ) S i Θ 2 (1 -S i ) 1 -(ΘS i ) 2 (52)
T 0 = Θ (1 -S e ) (1 -S i ) 1 -(ΘS i ) 2 (53)

G. Bidirectional reflectance of a contaminated slab overlaying a semi-infinite granular media

In realistic conditions, a slab will receive a collimated radiation from the solar disk, and a diffuse radiation from the granular medium underneath. There is a coupling between the two layers, illustrated on Figure 1. Using adding doubling formulas [START_REF] Douté | A multilayer bidirectional reflectance model for the analysis of planetary surface hyperspectral images at visible and near-infrared wavelengths[END_REF], we can express the total diffuse reflectance of the slab over a granular substrate as :

R Di f f = R 0 + T 0 T 0 r s ∞ ∑ n=0 (R 0 r s ) n = R 0 + T 0 T 0 r s 1 -R 0 r s ( 54
)
where

r s = 1- √ 1-ω s 1+ √ 1-ω s
is the lambertian reflectance of the substrate [START_REF] Douté | A multilayer bidirectional reflectance model for the analysis of planetary surface hyperspectral images at visible and near-infrared wavelengths[END_REF]. ω s is the single scattering albedo of the granular substrate. The last step is to simulate the diffuse contribution for one measurement. The total reflectance (BRDF) of the surface measured by the instrument is the sum of the specular and diffuse contributions :

R tot = R spec + R Di f f ( 55
)
where R spec is determined by Eq. 13 and R Di f f by Eq. 54.

DISCUSSION

A. Energy conservation At the first interface

We checked the conservation of the energy at different points in the model. We first checked it at the first interface, as it contains a complex numerical integration. To test the conservation of the energy at the first interface, we force the value of the Fresnel's reflection coefficient r f in Eq. 16 to one. Thus, all the energy is supposed to be sent back, and we have to obtain Q = 1 to have the energy conserved, where

Q = π 2 0 2π 0 cos α 2 a (ϑ s , ζ s ) S i, e, ψ, θ cos i det J g i (e, ψ) de dψ (56)
that is Eq. 16 where the Fresnel reflection coefficient is put to one. Figure 3 shows the value Q as a function of the incidence angle. Different roughness parameters θ were tested, ranging from θ = 0.01°to θ = 45°. Only values ranging from θ = 0.15°t o θ = 3.5°are displayed on Figure 3. This test illustrates the dependance of the validity of the model on both the incidence angle and the roughness parameter. To test the conservation of the energy for the whole model, we first have to set the complex value of the optical constant of the slab and the granular substrate to 0, to make the surface non absorbent. Then we integrate the energy sent back toward the sky. This energy must equal the incoming energy. To test this practically in the model, we set the sensor's angular aperture to a value that is equal to the integration step. where L is the surface radiance (W.m -2 .sr -1 ) A is the surface of a pixel , F is the incident flux in the incident direction (W.m -2 ) and i and e are the incidence and emergence angles. The relation R = π L F cos i between the reflectance factor and the radiance brings The symmetry of the model in azimuth leads to the quantity Q displayed in Figure 4. This figure shows that it is mostly the roughness parameter θ and the incidence angle i that control the validity of the model. Figure 5 shows the error in the energy conservation in percent, as a function of the roughness parameter θ and the incidence i. This gives the range of validity of the model according to a given tolerance. Roughness parameters larger than θ = 11°, always exhibiting error larger than 10 %, are not represented. For small slab real optical index (i.e. close to one), these errors decrease.

A.1. Slope distribution

As mentioned in section 2, this model is limited to the case of small θ. Figure 5 gives a quantification of that limitation. This is mostly due to the fact that the probability density function a (ϑ, ζ) that defines the repartition of slopes only makes sense if

(ϑ,ζ) a (ϑ, ζ) dϑ dζ = 1, which means that 1 2π π 2 0 a (ϑ) dϑ = 1 (59)
or that the value of

I = π 2 0 2 π tan 2 θ exp -tan 2 ϑ
π tan 2 θ sec 2 ϑ sin ϑ dϑ must be equal to 1. Figure 6

B. Behavior of the model

B.1. Specular reflection

This model is built to simulate observations. Thus, the specular spot characteristics will depend not only on the illumination divergence and the geometry, but also on the observation device. This model is designed to be adaptable to both conditions. Figure 7 shows a zoom on the specular spot for a water ice slab at 1 µm with a roughness parameter θ = 0.5°, illuminated at an incidence angle i = 50°illuminated with different light sources, and observed with two distinct detectors. In the first case (in Figure 7a), the surface is illuminated with a light source that has an aperture of 0.4°and observed with a sensor that has an aperture of 4.2°. It represents the conditions of a laboratory measurement with the instrument described in [START_REF] Brissaud | Spectrogonio radiometer for the study of the bidirectional reflectance and polarization functions of planetary surfaces. 1. design and tests[END_REF]. In the second case (in Figure 7b), the surface is illuminated with a light source that has an aperture of 0.2°and observed with a sensor that has an aperture of 6.92.10 -2 °. It represents the conditions of a measure with the OMEGA imaging spectrometer instrument orbiting the planet Mars [START_REF] Bibring | OMEGA: Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité[END_REF]. Both cases represent actual measurement situations.

As shown on Figure 7, both the amplitude and the shape of the specular spot depend on the characteristics of the illumination.

B.2. Influence of the parameters

To give a feeling on how the model behaves according to the different parameters, we chose a set of parameters, and plotted the dependence of the reflectance to the variation of one parameter around this first set. We chose arbitrary optical constants for the matrix and inclusions. We chose for the matrix n = 1.3 and k = 1.10 -3 , that are approximately the values for water ice at 270 K, and at the 1 µm wavelength. We selected as our standard set of parameter a 10 mm thick slab layer containing 1000 ppmv of 100 µm wide inclusions, overlaying a semi infinite granular layer of the same nature that the matrix. We tested the behavior of the model for various types of inclusions. We describe two types of behavior. The first type is when the absorption coefficient of the inclusions is smaller than the one of the matrix. This includes the particular case of a matrix contaminated with bubbles. It is illustrated in the figure 8, 10 and 9 with the green curves. In this case, the real part n of the optical index of the inclusions has very little influence. The second case is when the absorption in the inclusion is bigger than in the matrix (blue and red curves). In this case, the real optical index of the inclusions plays an important role.

Figure 8 shows the dependence of the reflectance on the thickness of the slab layer in different cases. In the case of an uncontaminated slab, the reflectance approaches 0 when the thickness increases. On the contrary, the reflectance of a slab containing inclusions will saturate at a value depending on the properties of the impurities. In the first case of a low absorption in the inclusions (green curve) the reflectance is higher than the reflectance of an uncontaminated layer whatever thickness the slab has. On the other case of high absorption in the inclusions, when the slab layer is thin, then the reflectance is lower than one of an uncontaminated slab, but as thickness increases, the value saturates, due to the scattering of light by the inclusions. This value gives an idea of the penetration depth of the light into a contaminated slab layer (i.e. the depth from which the layer becomes optically thick). Figure 9 show the dependence of the reflectance on the radius of inclusions in the slab layer. This illustrates the scattering properties of the inclusions. The reflectance factor of an uncontaminated at this geometry is approximately R = 0.372. In this figure, the volumetric proportion of inclusions in the matrix is held constant, so as the grain size increases, the number of inclusions per unit of volume decreases. Thus the scattering power decreases as well as a function of the grain size. In the case of inclusions with higher absorption than the matrix, at some point, the grain size reaches a value where the absorption in the inclusions becomes more efficient than the scattering effect, and the reflectance falls below the reference value of the uncontaminated slab (around 10 µm for the blue curve and 20 µm for the red one). Then, the decreasing probability of encountering an inclusion when grain sizes become too high make the reflectance approach the value of the uncontaminated slab. Indeed, when the grain size of the inclusion equals the thickness of the layer (at the extreme right of the plot), the probability of encountering one, knowing the the volumetric proportion is 1000 ppmv becomes very low and the influence of the inclusions become negligible.

Figure 10 shows the evolution of the reflectance in different cases where scattering or absorption dominates. On Figure 10a, the absorption is the dominant effect. In the case of inclusions with a higher absorption coefficient than the matrix, this makes (a) (b) Fig. 7. Zoom on the specular spot for a water ice slab at 1 µm with a roughness parameter θ = 0.5°, (a) illuminated at an incidence angle i = 50°with a light source that has an aperture of 0.4°, and observed with a detector that has an aperture of 4.2°, such as in the conditions of a laboratory measurement with the instrument described in [START_REF] Brissaud | Spectrogonio radiometer for the study of the bidirectional reflectance and polarization functions of planetary surfaces. 1. design and tests[END_REF], and (b) illuminated at an incidence angle i = 50°with a light source that has an aperture of 0.2°, and observed with a detector that has an aperture of 6.92.10 -2 °, such as in the conditions of a measure with the high resolution spectro imaging imaging instrument OMEGA orbiting Mars described in [START_REF] Bibring | OMEGA: Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité[END_REF]. The reflectance decrease as the roughness increase. This is due to the fact that a bigger roughness means more facets in specular conditions and thus less energy inserted into the system. The roughness parameter have a smaller impact on the diffuse reflectance than the other parameters of the model. the reflectance drop when the proportion of inclusions increases (blue and red curves). For the green curve, both absorption and scattering contribute in increasing the reflectance of the slab, thus it increases with the proportion of inclusions. Figure 10b shows more complexity. It is the same as Figure 10a except that the grain size of the inclusions is 10 µm instead of 100 µm. The green curve still represents the case of a lower absorption in the inclusion. Absorption and diffusion both tend to increase the reflectance, thus it increases with the proportion of inclusions. The red curve represents a case of higher absorption in the inclusions, when the scattering contribution is dominant. In this case, as diffusion limits the penetration of light into the layer, the reflectance increases with the proportion of impurities. The blue curve represents the limit case where diffusion and absorption contributions are of the same order of magnitude, leading to strong non-linear behavior.

Figure 11 shows the dependence of the reflectance of a 20 mm thick slab of water ice containing 1000 ppm of 100 µm wide inclusions on the roughness parameter θ. The optical indexes of the inclusions are n = 1.1 and k = 1.10 -9 . The diffuse reflectance of a slab decreases as the roughness increases. A bigger roughness means a bigger diversity in the slope distribution at the surface. This leads to an increased number of facets satisfying the specular reflection conditions defined in section A. Finally, less energy is inserted into the surface, and the diffuse reflectance is smaller. The dependance of the diffuse reflectance on the roughness if smaller compared to the dependance on the other parameters. This can be attributed to the relatively small range of values of θ. On the contrary the roughness parameter θ has a strong influence on the specular contribution.

CONCLUSIONS

We developed a radiative transfer model to simulate the bidirectional reflectance of a rough slab with inclusions. Typical calculation time is 1.10 -2 s per spectrum, considering 10000 wavelengths. Nevertheless, it can vary greatly depending on the sets of parameters desired.

Most of the constituting elements of this model have already been numerically [START_REF] Douté | A multilayer bidirectional reflectance model for the analysis of planetary surface hyperspectral images at visible and near-infrared wavelengths[END_REF] or experimentally validated [START_REF] Hapke | Theory of Reflectance and Emittance Spectroscopy[END_REF] but we adapted them to build a new model. In this study we tested numerically the conservation of energy and characterized the domain of validity of the model. We conducted sensibility studies in the case of a matrix containing only one type of inclusion. This shows the complexity and non linearity of the model with respect to its parameters. The sensibility study in the case of several types of inclusions were not conducted because of the great number of parameters. The experimental validation will be conducted in a following paper.

This model is designed to analyze massive hyperspectral data in the planetary science domain. In our favorite application, it calculates the radiative transfer in a contaminated ice slab overlaying an optically thick granular medium. The contamination in the slab can be of any type : ice, minerals or even bubbles. The matrix can be constituted as well of any ice. Thus our model can be applied on Earth with water ice, but also on Mars polar region covered with CO 2 ice [START_REF] Leighton | Behavior of carbon dioxide and other volatiles on mars[END_REF], on icy bodies, such as Jupiter's moon Europa (water ice), Neptune's moon Triton or dwarf planet Pluto (N 2 ice). Other applications in biology or industry are possible, as soon as the optical constants of each material are known.

We considered in all the calculations every wavelength independently. Thus, a spectrum in any spectral range can be built by computing every wavelength contribution at very high spectral resolution. The final objective is the comparison of the simulation to actual data, for analysis purposes. This makes this approach suitable for any spectroscopic measurement of slabs (made of ice or other material), overlaying optically thick material (granular or other material), from laboratory to spatial probe measurement. For the planetary science case, these results will be down-sampled at the instrument's wavelength resolution, using its PSFs.

One major hypothesis in this work is that we suppose an isotropic behavior of the inclusions. In the future, we plan to add the particles phase function to improve this point. We also plan to normalize the probability distribution function describing the roughness of the surface, to extend the applicability of the model.

NOTATIONS

Fig. 1 .

 1 Fig. 1. Illustration of the radiative transfer in a rough slab. (a) Radiative transfer for a slab ice layer only. Anisotropic transit are represented in red, and named with a prime. On the top left : illustration of the reflections and transmission at the first interface, used in the calculations of variables S e and fact, Θ . (b) Illustration of the adding coupling. The granular and slab layers are artificially separated in this figure to help the understanding of the coupling.

Fig. 3 .

 3 Fig.3. Q as a function on the incidence angle i, when forcing the value of the Fresnel's reflection coefficient r f in Eq. 16 to one. A value of one means that the conservation of energy is respected. This figure shows that in the cases of a roughness parameter ranging from θ = 0.15°to θ = 3.5°, the energy is fairly well conserved for incidences below 85°, and a roughness parameter below θ = 2.5°. Thus, this model will not be applicable to very high incidences values, and roughness over θ = 2.5°F

Figure 4 π 0 π 2 0

 42 Figure 4 show the value of Q = 2 π π 0 π 2 0 R tot cos e sin e de dψ. Practically, the energy is conserved if this integral equals 1. Indeed, the energy conservation gives

R

  tot cos e sin e de dψ = 1 (58)

Fig. 4 .

 4 Fig.4. Q as a function on the incidence angle i, where E r is the energy sent back from the surface, integrated over the hemisphere, and E i the incident energy.

Fig. 5 .

 5 Fig. 5. Error in the energy conservation as a function of θ and i.

  shows the value of I versus the roughness parameter θ. The function a (ϑ, ζ) only makes sense as a probability function if I = 1. For values of roughness larger than θ = 2°, I begins to fall to values below one. In a further development, we could extend the applicability of the model by defining a new probability density function a Norm (ϑ, ζ), that would be the normalization of the function a (ϑ, ζ) : a Norm (ϑ, ζ) = 1 I( θ) a θ,ϑ, ζ .

Fig. 6 .

 6 Fig. 6. I versus θ. For values of θ larger than θ = 2°, the probability density function for the slopes begins to drop. At θ = 10°, I = 0.957.

Fig. 8 .

 8 Fig. 8. Reflectance factor of a slab of water ice containing various types of inclusions, at a wavelength λ = 1 µm, as a function of the thickness of the slab layer, other parameters fixed. Black curve is the uncontaminated reference, and colored curves represent different optical indexes of inclusions.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Reflectance factor of a slab of water ice containing various types of inclusions, at a wavelength λ = 1 µm, as a function of the grain size of the inclusions, other parameters fixed. Colored curves represent different optical indexes of inclusions. When the absorption coefficient is higher in the inclusions than in the matrix (blue and red curves), there is a competition between scattering and absorption. The dominant effect depend on the grain size. The reflectance factor of an uncontaminated at this geometry is approximately R = 0.372.

Fig. 11 .

 11 Fig. 11. Reflectance factor of a 20 mm thick slab of water ice containing 1000 ppm of 100 µm wide inclusions as a function of the roughness parameter θ, other parameters fixed. The optical indexes of the inclusions are n = 1.1 and k = 1.10 -3 .The reflectance decrease as the roughness increase. This is due to the fact that a bigger roughness means more facets in specular conditions and thus less energy inserted into the system. The roughness parameter have a smaller impact on the diffuse reflectance than the other parameters of the model.

  slope angle of a facet in specular conditions(°) Θ ik transmission factor of a type k inclusion Θ transmission factor of the slab containing inclusion under isotropic illumination Θ transmission factor of the slab containing inclusion under collimated illumination ν optical path (m) ρ k radius of a type k inclusion (m) σ k geometrical cross section for type k inclusions (m 2 ) σ ek extinction cross section for type k inclusions (m 2 ) σ e mean extinction cross section (m 2 ) σ sk scattering cross section for type k inclusions (m 2 ) σ e mean scattering cross section (m 2 ) τ optical depth of the matrix containing inclusions ψ azimuth angle (°) ω single scattering albedo of the matrix containing inclusions ω s single scattering albedo of the granular substrate Ω C Solid angle subtended by the sensor (sr) Ω S Solid angle subtended by the source (sr) a m absorption coefficient of the matrix a ik absorption coefficient a type k inclusion a (ϑ, ζ) probability of occurrence for the slope(ϑ, ζ) A Surface of a pixel (m 2 ) A f Surface of a facet (m 2 ) D thickness of the slab layer (m) D apparent length of the first transit through the slab layer for one ray (m) D mean apparent length of the first transit through the slab layer for one ray (m) e emergence angle (°) F incident power flux in the radiation's direction (W.m -2 ) i incidence angle (°) k m Imaginary part of the optical index of the matrix k ik Imaginary part of the optical index of a type k inclusion L radiance (W.m -2 .sr -1 ) n m Real part of the optical index of the matrix n ik Real part of the optical index of a type k inclusion N number of facets within a pixel : N 1 N spec number of facets within a pixel satisfying specular conditions N total density of inclusions inside the matrix N k density of inclusions of type k inside the matrix p probability or probability per unit of length P power (W) Q sk scattering efficiency for type k inclusions r bidirectional reflectance (sr -1 ) r ik diffusive reflectance for a type k inclusion r m diffusive reflectance of the matrix r s diffusive reflectance of the granular substrate r f Fresnel reflection coefficient r f = R 2 ⊥ + R 2 R ⊥ Fresnel reflectivites for perpendicular polarization R Fresnel reflectivites for perpendicular and parallel polarization R 0 reflection factor of the slab containing inclusion under isotropic illumination R 0 reflection factor of the slab containing inclusion under collimated illumination R 0 R 0 -S e R Di f f diffuse reflectance factor of a slab containing inclusion over a granular substrate R spec specular reflectance factor of a slab containing inclusion over a granular substrate R tot reflectance factor of a slab containing inclusion over a granular substrate S shadowing function S e external reflection coefficient of the slab under collimated illumination S e external reflection coefficient of the slab under isotropic illumination S i internal reflection coefficient of the slab under isotropic illumination S ek external reflection coefficient of a type k inclusion S ik internal reflection coefficient of a type k inclusion T 0 transmission factor of the slab containing inclusion under isotropic illumination T 0 transmission factor of the slab containing inclusion under collimated illumination

  2 sin i sin ϑ spec cos ζ spec + cos ϑ spec cos i sin ϑ spec cos ζ specsin i y e = 2 sin i sin ϑ spec cos ζ spec + cos ϑ spec cos i sin ϑ spec sin ζ spec z e = 2 sin i sin ϑ spec cos ζ spec + cos ϑ spec cos i cos ϑ speccos i