Zili He
email: zili.he@mines-albi.fr

Paule Lapeyre

Stéphane Blanco

Simon Eibner

Mouna El

Richard Fournier

Monte-Carlo estimation of geometric sensitivities in Solar Power Tower systems of flat mirrors

Optical optimizations for a Concentrated Solar Power (CSP) system are currently limited largely to gradient-free methods, since the gradient is hard to obtain by using the existing numerical optics tools available in the community. This article aims to build new algorithms of the Monte-Carlo type, which numerically estimate the gradient of power impacting the receiver with respect to the geometric parameters that characterize the geometric status of the heliostats in the heliostats field, for a Solar Power Tower (SPT) system. Physical models will be built for the specific intensity and also for its derivatives to the geometric parameters of the heliostats, also called geometric sensitivities of intensity. Similar to the intensity but with their own models, they are regarded as physical quantities emitted, absorbed and reflected in the system. They carry the perturbations of intensity as information, corresponding to the relations between the geometric parameters and the physical events in the SPT system: blocking, spillage, shadowing, etc. These relations will be distinguished and discussed.

Finally, not only the gradient but also the contributions of physical events (blocking, spillage, shadowing, etc.) to the gradient can be estimated for further sensitivity analysis and optimization processes.

Nomenclature

Concentrating solar power (CSP) technology is one of the promising options to replace fossil fuel resources. Solar Power Tower (SPT) is a mature technology, and commercial facilities that follow this principle are currently in operation [START_REF] Clifford | Advances in central receivers for concentrating solar applications[END_REF].

The optical system described by the heliostats field usually contributes to about 40% of the total cost of an SPT system [START_REF] Buck | Heliostat field layout improvement by nonrestricted refinement[END_REF]. The optical study of SPT systems has been an active research field [START_REF] Lorin | An educated ray trace approach to solar tower optics[END_REF][START_REF] Fw Lipps | User manual for the university of houston solar central receiver, cellwise performance model: Ns[END_REF][START_REF] Li | Optics of solar central receiver systems: a review[END_REF]. Evaluations and optimizations of optical performance for SPT systems are essential.

As regards the evaluation of optical performance, there are two common categories of optical modeling methods: Monte-Carlo ray-tracing methods [START_REF] Manuel | Preliminary validation of tonatiuh[END_REF][START_REF] Caliot | Validation of a monte carlo integral formulation applied to solar facility simulations and use of sensitivities[END_REF][START_REF] Leary | User's guide for mirval: a computer code for comparing designs of heliostat-receiver optics for central receiver solar power plants[END_REF][START_REF] Wendelin | Soltrace: a new optical modeling tool for concentrating solar optics[END_REF] and the cone optics convolution-based methods [START_REF] Schwarzbözl | Visual hflcal-a software tool for layout and optimisation of heliostat fields[END_REF][START_REF] Charles | User's guide to helios: a computer program for modeling the optical behavior of reflecting solar concentrators[END_REF][START_REF] Michael | Solarpilot: A power tower solar field layout and characterization tool[END_REF]. Compared to the latter kinds of method, the Monte-Carlo ray tracing methods are more versatile and can reproduce the real intersections of photons in complex geometries [START_REF] Garcia | Codes for solar flux calculation dedicated to central receiver system applications: A comparative review[END_REF]. Several Monte-Carlo-based tools are compared to the convolution methods in the study [START_REF] Wang | Verification of optical modelling of sunshape and surface slope error for concentrating solar power systems[END_REF], revealing that the Monte-Carlo methods have better performance in accuracy.

However, the optimizations of optical performance, which are usually related to geometric optimizations are limited largely to gradient-free methods [START_REF] Pitz-Paal | Heliostat field layout optimization for high-temperature solar thermochemical processing[END_REF][START_REF] Farges | Life-time integration using monte carlo methods when optimizing the design of concentrated solar power plants[END_REF]. An important reason for choosing gradient-free methods is usually that 'no gradient information is available from Monte-Carlo ray tracing methods' [START_REF] Diago | Net power maximization from a faceted beam-down solar concentrator[END_REF][START_REF] Zhou | Heliostat field layout design for solar tower power plant based on gpu[END_REF]. The gradient-based methods are usually used when the dimension of the parameter space of optimization is large [START_REF] Favennec | Spacedependent sobolev gradients as a regularization for inverse radiative transfer problems[END_REF], which is often the case in CSP optimization. Limited by the dimension of the parameter space, researchers can only perform optimization with small numbers of parameters. Diago et al. optimized only the height and the diameter of the tertiary concentrator in a Beam-Down system. In order to decrease the number of parameters to be optimized, Yu et al. optimized the aiming strategy by regrouping the heliostats, which point to some fixed points on the receiver. Furthermore, the gradient-free methods are usually treated as a "black box" tool, and the knowledge of the gradient helps to understand the parameters' influences on the optimization's objective.

Only a few studies have applied gradient-based methods to optimize the geometry of solar concentrators. Yang et al. [START_REF] Yang | Pattern-free heliostat field layout optimization using physics-based gradient[END_REF] suggested approximating it by using area ratios. However, this method introduces an empiric weight. In studies carried out by Marston et al., the gradient is approximated by a finite difference [START_REF] Aj Marston | Geometric optimization of concentrating solar collectors using monte carlo simulation[END_REF][START_REF] James | Geometric optimization of solar concentrating collectors using quasi-monte carlo simulation[END_REF], which is also called the resimulation method by Gobet [START_REF] Gobet | Stochastic differential equations and Feynman-Kac formulas[END_REF]. However, beyond the fact that these studies deal only with simple geometrical situations, these authors claim that associated uncertainties are difficult to control. Regarding complex geometries such as a heliostat field in an SPT system with a large number of parameters, they cannot be treated in such a way (using finite difference approximation) because of the huge computational time and the lack of accuracy. Besides this approach of approximation, in the framework of Integral Formulation Monte-Carlo (IFMC) [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF], the impacting power on the receiver in an SPT system is written as an integral formulation based on the physics of radiation. By differentiating and reformulating the corresponding integral formulation (for example, with respect to the angular spreading coefficient [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF]), the gradient is also formulated as an integral formulation and can be estimated through

H i R (a) Spillage effect. H i H i R (b) Blocking effect. H i H i R (c) Shadowing effect.
Figure 1: Spillage, blocking, and shadowing effects are common physical events in the SPT system. In this work, we further classify them into spillage, forwardblocking, backward-blocking, forward-shadowing, and backward-shadowing effects to estimate and analyze the sensitivity of geometric parameters. In the figure, H i and H i are the heliostats and R is the receiver.

the Monte-Carlo method [START_REF] De Lataillade | Monte carlo method and sensitivity estimations[END_REF]. However, when the differentiated parameter is a geometric parameter that interacts in the integral domain (for example, the size of heliostats), great formal development efforts are needed [START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF]. Integrating this approach to a complex system such as the heliostats field of an SPT system seems nowadays to be impracticable. A very recent work [START_REF] Lapeyre | Monte-carlo and sensitivity transport models for domain deformation[END_REF] proposes an alternative to the work performed by [START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF]. This new approach proposes that the intensity and the geometric sensitivity of intensity (the derivative of the intensity with respect to a geometric parameter) should be regarded as two physical quantities emitted, transported, and intersected in a radiative system [START_REF] Lapeyre | Monte-carlo and sensitivity transport models for domain deformation[END_REF]. Based on this point of view, we can benefit from years of research in the radiative transfer field, on modeling the intensity to model the geometric sensitivity in an SPT system, making estimating the gradient in such a complex system possible. In parallel to our works, the literature shows that shape sensitivity calculations using the radiative transfer equation and its derivative has been the subject of much work in the computer graphics community [START_REF] Zhang | Path-space differentiable rendering[END_REF][START_REF] Zhang | A differential theory of radiative transfer[END_REF][START_REF] Li | Differentiable monte carlo ray tracing through edge sampling[END_REF][START_REF] Zeltner | Monte carlo estimators for differential light transport[END_REF]. Their objective is to reconstitute shapes to create the desired image by inversion. For many applications, such as remote sensing or deep neural networks for artificial intelligence [START_REF] Kato | Differentiable rendering: A survey[END_REF] this "differentiable rendering" will become increasingly important in the coming years in this community. While we can take advantage of computer graphics tools to deal with the complex geometry of an SPT system, we still need to understand how the phenomena in SPT systems (blocking in Fig. 1b, shadowing in Fig. 1c, and spillage in Fig. 1a) influence the sensitivities, which is also the objective of this work. Therefore, we will build a model of sensitivity and analyze the sources of the sensitivities in this model, which is specific to our approach [START_REF] Lapeyre | Monte-carlo and sensitivity transport models for domain deformation[END_REF][START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]. The sensitivities are propagated instead of the intensities following the Radiative Transport Equation (RTE) and the sensitivity sources are identified by this approach on the borders or on the surfaces of the studied system when we change geometric parameters. Consequently, we will see in this work that sensitivities can be presented in an intuitive manner: displacements and rotations of heliostats create sensitivity sources on the edges and surfaces of the heliostats that move according to different parameters of displacement and rotation.

As stated above, this paper aims to develop a method to estimate the 1: Sensitivity and its contributions (estimated values and their standard deviations) for the translation X of the 676 th heliostat. For the 676 th heliostat in Fig. 2 (indexed as 676 in the table), the sensitivity of its translation X (indexed as 1 in the table) has two contributions. The translation will block its neighboring heliostat in its back (S blo 676,1 is negative). However, the translation helps it to be blocked less by its neighboring heliostat on its front (S tar 676,1 is positive). The sensitivity of translation X for this heliostat is then the sum of these two contributions (S 676,1 ). geometric gradient (sensitivities) of SPT systems and to indicate how the blocking, spillage, and shadowing effects contribute to the gradient (sensitivities). The implementation of this method leads to different results that will be of interest to the SPT community to design (Fig. 2a,2b,2c,2d) solar power plants, and also for their remote control (Fig. 2e,2f). In Fig. 2, x-y represents the ground plane, z is the height of the heliostat position, and each point is dedicated to a heliostat in the field. Furthermore, for each parameter of each heliostat, a detailed analysis of the contributions (spillage, blocking, shadowing) of sensitivity can be carried out (Table .

1).

Finally, the article is organized into three parts:

-Section.2 provides the models of transport for the intensity and its derivative towards each heliostat five degrees of freedom (three translations and two rotations) and its size.

-Section.3 provides the integral formulations and algorithms to estimate the power collected by the receiver and the different sensitivities towards these six geometric parameters for each heliostat.

-Section.4 is dedicated to validation purposes: the estimation of the different sensitivities is compared to the approximation by the finite difference method. Comparisons and discussions about these two methods will be undertaken.

General models

In this work, we address the question of estimating the power P impacting the receiver of an SPT system with flat squared heliostats (Fig. 3), as well as its derivatives S i,j ≡ ∂ πi,j P of geometric parameters π with i the index of a heliostat in the heliostats field and j the index of a geometric parameter (j ∈ {1, 2, 3, 4, 5, 6}, see Fig. 4). In order to numerically estimate P , we will build a model following the standard radiative transfer physics, using the intensity I as the descriptor. The transport of I in the medium is governed by a partial differential equation at any The corresponding color indicates its sensitivity with respect to its positions, rotations, and lengths. (10 6 Monte-Carlo realizations on each border of each heliostat and 10 8 Monte-Carlo realizations on each surface of heliostats).

H i • • • H i+1 H i+2 • • • R x r n r ω r
n i a i,φ a i,θ H i a i,θ
Figure 4: A heliostat (H i ) has 5 degrees of freedom, which includes the translation (π i,1 ,π i,2 and πi,3 ), following the three axes in the global coordinate system ( e 1 , e 2 and e 3 ) and the rotation (π i,4 and πi,5 ) corresponding to an altazimuth mount according to the two vectors of rotation ( a i,θ and a i,φ ). Also, the length of the size of the heliostat is included in this study (π i,6 ). Those six parameters are indexed by j. locations x inside the field G, following a vector of direction ω in a unit sphere S (RTE). On the boundaries ∂G, the transport of I is governed by the boundary conditions, where the locations are noted with index, see Section.2.4 for details.

Equivalent to P , in order to numerically estimate S i,j , we will build a radiative model in the medium and on the boundaries using the geometric sensitivity of intensity s ≡ ∂ πi,j I as the descriptor [START_REF] Lapeyre | Monte-carlo and sensitivity transport models for domain deformation[END_REF]. This descriptor carries the perturbation of I as information, corresponding to the relations between the geometric parameters πi,j and the physical events in the SPT system: blocking, spillage, shadowing, etc.

I and s i,j will be modeled locally in the medium and boundary conditions. After that, P and S i,j can be observed by formulating integral formulations of I and of s i,j .

Definitions

Geometric sensitivity of intensity

The specific intensity used to model the STP system is considered independent of the wavelength [START_REF] Wang | Verification of optical modelling of sunshape and surface slope error for concentrating solar power systems[END_REF][START_REF] Blanc | Direct normal irradiance related definitions and applications: The circumsolar issue[END_REF]. It is defined as I ≡ I( x, ω, π), where π is the matrix composed of πi,j (see Fig. 4):

π ≡ [π i,j ] . (1)
x, ω and π are independent variables and matrix s is defined as the derivative of I with respect to π [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]:

s( x, ω, π) ≡ [s i,j ( x, ω, π)] = ∂ π I( x, ω, π).
(

) 2
The components of s are:

s i,j ( x, ω, π) = ∂ πi,j I( x, ω, π). (3)

Geometric sensitivity of power

We implement this work within the framework of IFMC [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF]. P can be estimated by the product of spatial integral and angular integral of incoming I on the receiver:

P (π) = R d x r 2π | ω r • n r |d ω r I( x r , ω r , π) (4)
where x r is a position on the receiver, ω r the direction of observation, n r the normal of the receiver, and R the absorbing surface of the receiver (see Fig. 3).

Matrix S is defined as the derivative of P with respect to π, which can be developed in a straight forward way since x r and ω r are independent of π that characterizes only the heliostats:

S(π) ≡ [S i,j (π)] = R d x r 2π | ω r • n r |d ω r s( x r , ω r , π). (5)
The components of S are:

S i,j (π) = R d x r 2π | ω r • n r |d ω r s i,j ( x r , ω r , π). (6)
We aim to estimate P (π) and S i,j (π) by building a model of I and a model of s i,j in an SPT system.

Assumptions

The following assumptions are made for modeling I and s i,j in an SPT system:

1. The Pill-Box model is used to describe the sun shape (Eq.15).

2. The heliostats are considered perfectly specular, and beam pointing error is neglected.

3. The reflectivity of heliostats is considered homogeneous on the surface.

4. The medium between heliostats and the receiver is considered transparent.

It is noted that even though they are typical assumptions for modeling an SPT system, more realistic modelings can be realized based on the theoretical studies [START_REF] Lapeyre | Monte-carlo and sensitivity transport models for domain deformation[END_REF][START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF] without specific difficulties (more realistic sun shape model, spatially and angularly inhomogeneous reflectivity, imperfect reflection, semitransparent medium, etc.). In this work, we aim to build a basic model for I and s i,j in an SPT system. More realistic assumptions can be taken into account based on the developed basic model in this work. For example, if an angularly inhomogeneous reflectivity is considered in future work, it will have extra influences on the sensitivity with respect to the rotations of the heliostats. These extra influences can be considered as an additional source of sensitivity based on the developed model in this work.

Transport in the medium

Transport of intensity

The transport of intensity in the medium is modeled by RTE [START_REF] John R Howell | Thermal radiation heat transfer[END_REF]. In an inhomogeneous medium, the collisions of absorption and scattering are described by the collision operator C [START_REF] Zhang | A differential theory of radiative transfer[END_REF]:

∂I( x, ω, π) ∂ x = C[I( x, ω, π)] + κ a ( x)I b ( x) (7)
with

C[I( x, ω, π)] = -κ a ( x)I( x, ω, π) -κ s ( x)I( x, ω, π) + 1 4π ω =4π κ s ( x)I( x, ω , π)Φ( ω , ω)d ω (8)
where κ a is the absorption coefficient, κ s the scattering coefficient, x the position vector in the medium, ω the vector of propagation and Φ( ω , ω) is the phase-function. Under the assumption of local thermal equilibrium, the source κ a ( x)I b ( x) (in Eq.7) is isotropic and is a function of temperature, and I b is the specific intensity at equilibrium following the Planck's law.

In this work, we will consider a transparent and cold medium (where κ s = 0, κ a = 0, and I b ( x) = 0). Therefore, the equation of transport for I used in this work is:

∂I( x, ω, π) ∂ x = 0 (9)

Transport of geometric sensitivity of intensity

We differentiate Eq.7 with respect to the geometric parameter π on the both sides:

∂ 2 I( x, ω, π) ∂ x∂ π = ∂ C[I( x, ω, π)] ∂ π (10)
Referring to the definition of the sensitivity (Eq.2), the following equation transport is yielded:

∂s( x, ω, π) ∂ x = C[s( x, ω, π)] (11)
The development from Eq.7 to Eq.11 is relatively straightforward. First, it is because the medium's properties (κ a , κ s , Φ and I b ) are not dependent on π, and secondly, x, ω and π are three independent variables.

Referring to Eq.2 and Eq.3, the equation of transport of each of the component s i,j is [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]:

∂s i,j ( x, ω, π) ∂ x = C[s i,j ( x, ω, π)] (12)
with C the same collision operator as it is in Eq.8:

C[s i,j ( x, ω, π)] = -κ a ( x)s i,j ( x, ω, π) -κ s ( x)s i,j ( x, ω, π) + 1 4π ω =4π κ s ( x)s i,j ( x, ω , π)Φ( ω , ω)d ω (13)
Finally, the collision operator C is also applied to the geometric sensitivities when describing its transport in the medium. s i,j can be regarded as a physical quantity transported in the medium. Similar to the intensity I, the sensitivity of intensity s i,j is also absorbed and scattered. This new point of view allows us to study the geometric sensitivity of intensity practically, benefiting from years of research in the transport of intensity in the physics of radiative transfer.

In this work, we will consider a transparent and cold medium. The intensity and the geometric sensitivity will have the same model of transport in the medium:

∂s i,j ( x, ω, π) ∂ x = 0 (14)
In order to create a similar physical picture for the transport of I and s i,j in the medium, we can consider that the medium is also 'transparent' for the geometric sensitivities s i,j .

Boundary conditions

I and s i,j are transported following Eq.9 and Eq.14 in the medium, while on the boundaries, they can be emitted, reflected, or absorbed.

Seeing Fig. 5, positions and directions on the boundaries are noted with an index of s on O s , r on R, p on H + i , b on H - i and O l on O l . We will build the boundary conditions for I and then differentiate them with respect to πi,j to establish the boundary conditions for s i,j .

Upper outline O s

In an SPT system, the only source of intensity is the sun. Most of the time, we use the sun-shape models, which are usually angular functions, to describe the arriving intensity in the system [START_REF] Wang | Verification of optical modelling of sunshape and surface slope error for concentrating solar power systems[END_REF][START_REF] Buie | Sunshape distributions for terrestrial solar simulations[END_REF][START_REF] Charles | Six-gaussian representation of the angular-brightness distribution for solar radiation[END_REF]. To be noticed that this work is a stationary study, which means it is established for a fixed solar position.

The pill-box sun shape is used in this work:

H i • • • H i+1 H i+2 • • • R O s O l x p x b x r x x O l x s n r n i ω s ω r ω ω p ω b ω O l Figure 5:
The system consists of heliostats H i (with index of i) on the ground with H + i its reflecting surface and H - i its rearward surface, a receiver R, an upper outline O s and an lower outline O l . Solar intensity is reflected by H + i and collected on the R.

I( x s , ω s , π) ≡ I( ω s ) = I 0 H( ω s • ω c -cos (θ disk )) (15)
where ω c is a vector that characterizes the solar position (fixed in this work), I 0 the constant intensity coming from the sun, H the Heaviside function and θ isk = 0.0046 rad the radial angular radius of the solar disk. Moreover, this boundary condition is noted I because it is the source of I in the system.

The upper outline O s can be considered emitting intensity towards the system, but it is not reflecting the incoming intensity. Therefore, this boundary can be considered as a 'black body' for I.

The derivative of I on O s with respect to πi,j is therefore null, because πi,j characterizes the geometry of heliostats instead of O s :

s i,j ( x s , ω s , π) = 0. ( 16
)
Following Eq.16, no s i,j is emitted by the sun, and all s i,j transported to the upper outline O s are not reflected but absorbed. Therefore, in order to create a similar physical picture for the transport of I and s i,j on the boundaries, O s can be regarded as a 'cold black body' for s i,j .

Receiver R and lower outline O l

R and O l in our system (Fig. 5) are considered 'cold black body' for I. All the intensity is not hitting the heliostats, and the receiver will be lost. The radiative emission of the receiver is not included in the model since the heat transfer of the receiver is out of the scope of this work.

The boundary conditions of receiver R and lower outline O l are then established:

I( x r , -ω r , π) = 0 I( x O l , -ω O l , π) = 0 ( 17
)
whereω r andω O l are the outgoing direction of propagation of I from the receiver and the lower outline (Fig. 5).

Since πi,j only characterizes the heliostats, there are no s i,j going out from these boundaries:

s i,j ( x r , -ω r , π) = 0 s i,j ( x O l , -ω O l , π) = 0. ( 18
)
In order to create a similar physical picture for the transport of I and s i,j on the boundaries, R and O l can also be regarded as 'cold black body' for s i,j . All s i,j that are not hitting the heliostat or the receiver will be lost, and s i,j are not emitted from these boundaries.

Heliostat H i

Intensity:

We define n i as the normal of the reflecting surface of H i . The boundary conditions of I on a heliostat H i are modeled separately by two parts: one for the reflecting surface H + i (where ω p • n i > 0), the other for the rearward surface H - i (where ω b • n i < 0), see Fig. 6. Also, the heliostats are considered 'cold', which means they do not emit intensity. The transport of intensity on the reflecting surface of the heliostat can be described by:

I( x p , ω p , π) = ρI( x p , ω s , π) ( 19
)
where ρ is the reflectivity, I( x p , ω s , π) the incoming intensity and ω p the direction reflected by heliostat of ω s (Fig. 6a).

The transport of I on the rearward surface of H i can be described by:

I( x b , ω b , π) = 0. ( 20
)
The rearward surface of H i is also considered a 'cold black body'.

Geometric sensitivity of intensity:

The development of the boundary conditions of s i,j for H i is less straightforward. [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF] proposed a model of geometric sensitivity on reflecting surfaces. Based on this model, the boundary conditions of s i,j is yielded for reflecting surface (Eq.21, referring to Fig. 6a) and also for the rearward surface (Eq.22, referring to Fig. 6b). It is found that the boundary condition of the reflecting surface contains emission terms (the sources) and a reflection term, see details in Appendix.A.

s i,j ( x p , ω p , π) =s i,j ( x p , ω p , π) +ρs i,j ( x p , ω s , π) (21)
s i,j ( x b , ω b , π) = si,j ( x b , ω b , π) (22)
We note si,j as the sources of s i,j (as we noted I as the source of I).

The sources si,j can be physically regarded as the local perturbations of I with respect to πi,j on H i , while s i,j is such perturbations transported in the radiative system (see Fig. 7). Usually, the sources of sensitivity si,j are functions of I, where the two models of transport are coupled.

H i H i R I O s O l (a) Source of intensity I H i H i R πi,1 I O s O l H i si,1 s i,1 I (b)
Source of sensitivity si,1: An example.

Figure 7:

The source of intensity is on the boundary O s . In this example, The source of sensitivity is on the border of heliostat H i because the perturbation of πi,1 will create locally a perturbation of I on the border of H i , which causes a perturbation of P . Also, the source si,j is coupled with the arrived intensity I. Therefore, the model of s i,j is coupled with the model of I. The optical path of s i,j is in red, and the coupled optical path of I is in brown.

Moreover, different types of perturbations (sources) are distinguished, related to blocking, shadowing, and spillage effects. Therefore, we further classify the sources by their contributions to these different physical effects.

For the reflecting surface

H + i : si,j ( x p , ω p , π) = star i,j ( x p , ω p , π) + sblo i,j ( x p , ω p , π) + sshad-b i,j ( x p , ω p , π) (23)
and for the rearward surface

H - i : si,j ( x b , ω b , π) = sshad-f i,j ( x b , ω s , π). ( 24
)
star i,j is named the source of targeting. It is the perturbation caused by the targeting or missing of the reflected intensity on H + i , with respect to πi,j . It is separated into two parts:

star i,j = star,spatial i,j + star,angular i,j . (25)
star,spatial i,j

is originated from the change of position of H i with respect to πi,j , located on the border of H + i (see Fig. 8 for an example of πi,1 ). star,angular i,j

is originated from the change of normal, located on the surface of H + i (see Fig. 9 for an example of πi,4 , but to be noted that perturbation of πi,4 leads to the change of position and also the change of normal).

Concerning the change of position, the sources star,spatial i,j are located on the border of the heliostat (see formal developments in Appendix.A and Appendix.B). If we note the sum of star,spatial i,j captured on R to be S tar,spatial i,j , it will be 0 when star,spatial i,j are all captured or all missed, see Fig. 8. However, when the spillage effect (Fig. 8c), forward-blocking effect (Fig. 8d) or backward-shadowing effect (Fig. 8e) occurs, the perturbation of πi,j will perturb these effects and star,spatial i,j will be partially captured. Therefore, S tar,spatial i,j will be non-null . Concerning the change of normal, the sources star,angular i,j are located on the surface of the heliostat (see details in Appendix.A and Appendix.B). If we note the sum of star,angular i,j captured on R to be S tar,angular i,j , it will be non-null when star,angular i,j is all captured on R since the perturbation of πi,j perturbs the cosine effect, and it will be 0 when all is missing, see Fig. 9. Moreover, when the spillage effect (Fig. 9c), forward-blocking effect (Fig. 9d) or backward-shadowing effect (Fig. 9e) occurs, the perturbation of πi,j will perturb these effects and star,angular i,j will be partially captured.

We must distinguish the contributions from backward and forward shadowing in the model of s i,j . Backward-shadowing means that the perturbation of πi,j affects the shadows on the heliostat itself Fig. 10a. Forward-shadowing means that the perturbation of πi,j affects the shadows on a neighboring heliostat Fig. 11.

sshad-b i,j is named the source of backward-shadowing. This perturbation is observed on the surface of H + i , but originated from the border of a neighboring heliostat, see Fig. 10a and10b. It is also separated by two parts, corresponding to the change of position and direction of

H i : sshad-b i,j = sshad-b,spatial i,j + sshad-b,angular i,j (26)
We note S shad-b,spatial i,j as the sum of the source sshad-b,spatial i,j captured on R. When backward-shadowing occurs, the perturbation of πi,j will perturb the effect of backward-shadowing by the change of position (see an example of πi,1 in Fig. 10a) and S shad-b,spatial i,j = 0. We note S shad-b,angular i,j as the sum of the source sshad-b,angular i,j captured on R. When backward-shadowing occurs, the perturbation of πi,j will perturb

H i R πi,1 H i star,spatial i,1 (a)
All the sources are captured.

H i R πi,1 H i star,spatial i,1 (b)
All the sources are missing

H i H i R πi,1 star,spatial i,1
(c) Sources are partially captured: Spillage effect.

H i H i R πi,1 star,spatial i,1
(d) Sources are partially captured: Forward-blocking effect.

H i H i R πi,1 star,spatial i,1
(e) Sources are partially captured: Backward-shadowing effect. are all captured by R or they are all missing, the perturbation of πi,1 is not causing the perturbation of P and S tar,spatial i,1 = 0. When spillage effect, forward-blocking effect, or backward-shadowing effect occurs, star,spatial (e) Sources are partially captured: Backward-shadowing effect is on the border of the neighboring heliostat. The backward-shadowing is about the shadowing effect created by the neighboring heliostat. the effect of backward-shadowing by the change of normal (see an example of πi,4 in Fig. 10b) and S shad-b,spatial i,j = 0. sshad-f i,j is named the source of forward-shadowing. This perturbation of I is caused by the creation of shadows. With a perturbation of πi,j , shadows might be created on another heliostat in the field.

We note S shad-f i,j

as the sum of the source sshad-f i,j captured on R. When forward-shadowing occurs, a perturbation of πi,j will perturb the effect of forward-shadowing (see an example of πi,1 in Fig. 11) and S shad-f i,j = 0. sblo i,j is named the source of backward-blocking. This perturbation of I is caused by blocking. With a perturbation of πi,j , H i might block the I reflected on another heliostat in the field.

We note S blo i,j as the sum of the source sblo i,j captured on R. When backwardblocking occurs, a perturbation of πi,j will perturb the effect of backward-blocking (see an example of πi,1 in Fig. 12) and S blo i,j = 0. Finally, the geometric sensitivity of power S i,j is the sum of all captured sources:

S i,j = S tar i,j + S blo i,j + S shad-b i,j + S shad-f i,j (27)
with

S tar i,j = S tar,spatial i,j + S tar,angular i,j (28)
and

S shad-b i,j = S shad-b,spatial i,j + S shad-b,angular i,j . (29)
H i R πi, 1
H i sshad-f i,1
Figure 11: When the forward-shadowing effect occurs, with a perturbation of πi,1 , a perturbation of P is observed. The source of sensitivity sshad-f

i,1
is on the border of the heliostat itself. The forward-shadowing is about the shadowing effect created by the heliostat itself.

H i H i R πi,1 sblo i,1
Figure 12: When the backward-blocking effect occurs, a perturbation of πi,1 causes the perturbation of P . This is because π characterizes how much intensity reflected is blocked by H i . The source of sensitivity sblo i,j is on the border of

H + i .
Each contribution in Eq.27 is the sum of the sources si,j captured by the receiver respectively:

    S tar i,j (π) S blo i,j (π) S shad-b i,j (π) S shad-f i,j (π)     = R d x r 2π | ω r • n r |d ω r     star i,j ( x r , ω r , π) sblo i,j ( x r , ω r , π) sshad-b i,j ( x r , ω r , π) sshad-f i,j ( x r , ω r , π)     . (30)

Algorithms

In the previous section, general models for the intensity I and its derivatives s i,j are built for an SPT system. Based on the models, we aim to build algorithms to estimate the impacting power P and the matrix S by the Monte-Carlo ray tracing method.

Concerning the model of intensity I, the only source is the intensity from the sun I located on the surface that delimits the sky, O s . P will be estimated by summing the source I captured on the receiver R(see Eq.4).

Concerning the model of s i,j , there are several sources of sensitivities on the heliostat H i . S i,j will be estimated by summing all the sources si,j captured on the receiver R(see Eq.5).

The corresponding statistical ray tracing process can start from the receiver looking for the sources or, inversely, from the sources located on the heliostat, looking for the receiver. The strategy of proceeding statistical ray tracing process has enormous influences on the convergence rate.

In order to estimate the impacting power P , unlike the Collision-Based Monte-Carlo algorithms (CBMC) where the ray-tracing starts from a plane above reflectors (O s in our case) [START_REF] Caliot | Validation of a monte carlo integral formulation applied to solar facility simulations and use of sensitivities[END_REF], we integrate this work in the framework of IFMC, starting the ray-tracing process from the reflecting surfaces (H + in our case), which improves convergence rate [START_REF] Wang | Verification of optical modelling of sunshape and surface slope error for concentrating solar power systems[END_REF][START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF].

In order to estimate the matrix of sensitivity S, we start the statistical ray-tracing process by the sources, looking for the receiver R.

Estimation of impacting power

The algorithm to compute the impacting power by IFMC has been thoroughly described by [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF][START_REF] Piaud | Application of Monte-Carlo sensitivities estimation in Solfast-4D[END_REF]. Here is a quick reminder of the integral formulation for estimating P under the assumptions made in Section.2.2:

P = H + p X ( x p1 )d x p1 Ωs p Ωs ( ω s )d ω s {H( x p0 ∈ O s )H( x p2 ∈ R) ŵ} (31) with p X ( x p1 ) = 1 S H + , (32)
p Ωs ( ω s ) = 1

C d ω s = 1 2π(1 -cos θ disk ) (33) ŵ = DN Iρ| ω s • n i |S H + (34)
where S H + is the area of all reflected surfaces of the heliostats field. Equation 31 is the integral formulation of the model to compute the impacting power P and can also be seen as an expected value which is the core of an IFMC algorithm. The expected value P can be estimated by the mean operator (Algo.2 in Appendix.E):

P ≈ 1 N N q=1 [H( x p0,q ∈ O s )H( x p2,q ∈ R) ŵq ] (35)
H i R O s O l x p1 x p0
x p2 n i ω s ω p Figure 13: In order to estimate P , for each realization, the algorithm starts by uniformly sampling a random position x p1 (following Eq.32) on the reflected surface of the whole heliostat field H + . Then, a direction is sampled within the solar cone Ω s of angular radius θ disk following the density of probability p Ωs (Eq.33). In order to check if the shadowing effect occurs, x p0 is defined as the first intersection with a boundary of the system, starting at x p1 in the direction ω s , and it is tested if x p0 belongs to O s . If there is no shadowing, the reflected direction ω p is computed based on the law of reflection for specular surfaces from the sampled ω s and the local normal n i at position x p1 . In order to check if the spillage or the blocking effect occurs, x p2 is defined as the first intersection with a boundary of the system, starting at x p1 in the direction ω p , and it is tested if x p2 belongs to R. Finally, if no shadowing nor blocking occurs, the weight of Monte-Carlo for this realization ŵ is computed. where x p0,q and x p2,q are obtained through the described ray-tracing algorithm on Fig. 13 for the q th realization. ŵq is the value of ŵ for the q th realization.

The standard deviation σ(P ) is computed simultaneously with P [40].

Estimation of the geometric sensitivities of power

The whole matrix of sensitivity of impacting power S has a dimension of n H ×n π , where n H is the number of heliostats and n π = 6 is the number of geometric parameters of each heliostat, referring to Fig. 4. When the number of heliostats n H becomes large, it will be very costly in computing time if we estimate them one by one. Herein, we will build an efficient algorithm to estimate S. The matrix of S is composed of the vectors of sensitivities for each heliostat S ≡ S i (reminding the index i for i th Heliostat in the field):

S i (π) ≡         S i,1 (π) S i,2 (π) S i,3 (π) S i,4 (π) S i,5 (π) S i,6 (π)         T (36)
Following Eq.27, each S i,j has 4 contributions. We defined the vectors of contributions as well as the vectors of sources similarly, and finally, we have:

S i = S tar i + S blo i + S shad-b i + S shad-f i (37) and si = star i + sblo i + sshad-b i + sshad-f i . (38)
The contributions will then be the sum of captured sources respectively on the R. The corresponding integral formulations are detailed in Appendix.B. The algorithm to estimate S and its standard deviation σ(S) is Algo.1.

Algorithm 1 Estimate S and σ(S)

Input: Geometries of heliostats and receiver, ω c ω c characterizes the solar position (Eq.15)

Initialization: i ← 0 while i < n H do n H number
+ S shad-f i σ( S i ) ← σ( S mat i ) + σ( S blo i ) + σ( S shad-b i ) + σ( S shad-f i ) i ← i + 1 end while Output: S ← S i , σ(S) ← σ( S i )
In the following part of this section, we will discuss the estimation of each vector of contribution in Eq.37.

Contribution of the effect of targeting S tar i

S tar i,j is a component of S tar i , estimated by the sum of the sources of targeting star i,j (π) arriving to the receiver. The corresponding integral formulation is yielded in Appendix.B and is separated into two parts, corresponding to the change of position (Fig. 8) and the change of normal (Fig. 9), noted S tar,spatial i and S tar,angular i respectively. The corresponding sources of sensitivities are separated by star,spatial i,j and star,angular i,j in Eq.25. Those sources correspond to the spatial gradient of intensity I and the angular gradient of I [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]. It is noted that the πi,4 and the πi,5 characterize the change of normal and also the change of positions (because they are related to the rotations). In contrast, other parameters characterize only the change of position. x p1 is on the 4 borders of the heliostat ∂H i having k as the index of the border, x p0 the first intersection starting by x p1 following -ω s , x p2 the first intersection starting by x p1 following ω p , ω s the direction on the border of solar cone ∂Ω s , ω p the direction of reflection of ω s , x p0 the first intersection starting by x p1 following -ω s and x p2 the first intersection starting by x p1 following ω p . For each realization of estimating S tar,spatial i , we sample x p1 and ω s . If x p0 hits O s and x p2 hits R we count w 1 (Eq.42). For each realization of estimating S tar,angular i , we sample x p1 and ω s . If x p0 hits O s and x p2 hits R we count w 2 (Eq.43).

H i R O s O l x p1 x p0 x p2 n i ω s ω p x p1 x p0 x p2 ω s ω p
The integral formulation for estimating S tar i is detailed in Eq. 39 and Eq. 40. Corresponding probability density functions are given in Eq. 41. In Eq. 40, k refers to the four sides of the squared heliostat.

S tar i (π) = S tar,spatial i (π) + S tar,angular i (π) (39)
                     S tar,spatial i (π) = 4 k=1 l k p X ( x p1 )d x p Ωs p Ωs ( ω s )d ω s w 1 H( x p0 ∈ O s )H( x p2 ∈ R) S tar,angular i (π) = H + i p Xi ( x p1 )d x p ∂Ωs p Ω s ( ω s )d ω s { w 2 H( x p0 ∈ O s )H( x p2 ∈ R)} (40)
where

p Xi ( x p1 ) = 1 S H + i ; p X ( x p1 ) = 1 l k ; p Ω s ( ω s ) = 1 2π (41)
w 1 = star,spatial i ( x p , ω p , π) p X ( x p1 )p Ωs ( ω s ) ; ( 42
)
w 2 = star,angular i ( x p , ω p , π) p Xi ( x p1 )p Ω s ( ω s ) . ( 43
)
The notation and the Monte-Carlo ray-tracing algorithm are shown in Fig. 14. See also Algo.3 for the complete algorithm and Appendix.C for the explicit expression of w 1 and w 2 .

Contribution of backward-blocking effect

S blo i S blo
i,j is a component of S blo i , estimated by the sum of the sources of backwardblocking sblo i,j (π) arriving to the receiver. The corresponding integral formulation is yielded in Appendix.B.

Unlike the source of targeting, the evolution of πi,j has no effects on the change of normal of the reflected intensity (Fig. 15). Intensity is reflected following n i on H + i , while πi,j characterizes the heliostat H i . Therefore, sblo i,j is only about the change of position.

H i R H i O s O l n i x p2
x p1

x p0 x p3 ω s ω p x p1 is the first intersection starting by x p2 , following -ω p , x p0 the first intersection starting by x p1 following -ω s , x p3 the first intersection starting by x p1 following ω p . For each realization of estimating S blo i , we sample x p2 and ω s , and we calculate ω p . If x p0 hits O s , x p1 hits H i and x p2 hits R we count w 3 .

Seeing Fig. 15 for the ray-tracing algorithm, the integral formulation for estimating S blo i is as following:

S blo i (π) = n H ,i =i i =1 4 k=1 l k p X ( x p2 )d x p2 Ωs p Ωs ( ω s )d ω s w 3 H( x p0 ∈ O s )H( x p1 ∈ H + i )H( x p3 ∈ R) (44)
where

w 3 = sblo ( x p2 , ω p , π) p X ( x p2 ) p Ωs ( ω s ) ( 45
)
with n H the number of heliostats in the field. Moreover, the blocked heliostat H i can be all other heliostats except H i in the heliostats field. In order to capture all the sources sblo i,j (π i,j ), we need to proceed ray-tracing test for all other heliostats except H i in the field to detect the blocking effect. However, most of them will not be blocked by H i since the distance between the two heliostats is large. In order to accelerate and optimize the algorithm, we will proceed with the ray-tracing tests on the heliostats, of which the distance to H i is within a limited length. This length is defined as l d . See Algo.4 for the complete algorithm and Appendix.C for the explicit expression of w 3 .

Contribution of forward-shadowing

S shad-f i,j S shad-f i,j is a component of S shad-f i,j
, estimated by the sum of the sources of forward-shadowing sshad-f i,j (π) arriving to the receiver. The corresponding integral formulation is yielded in Appendix.B.

Seeing the ray-tracing algorithms on Fig. 16, S shad-f i can be estimated by the following integral formulation:

S shad-f i (π) = 4 k=1 l k p X ( x p1 )d x p1 Ωs p Ωs ( ω s )d ω s w 4 H( x p0 ∈ O s )H( x p2 ∈ H + )H( x p3 ∈ R) , ( 46
)
where

w 4 = sshad-f i ( x p1 , ω p , π) p X ( x p1 )p Ωs ( ω s ) . ( 47
)
See also Algo.5 for the complete algorithm and Appendix.C for the explicit expression of w 4 . (π) arriving to the receiver. The corresponding integral formulation is yielded in Appendix.B. Similar to the targeting effect, the evolution of πi,j might have two possible impacts on the heliostat H i : the change x p0 is the first intersection starting by x p1 , following -ω s , x p2 the first intersection starting by x p1 following ω s , x p3 the first intersection starting by x p2 following ω p . For each realization of estimating S blo i , we sample x p1 and ω s , and we calculate ω p . ω p is reflected following the normal n i which is the normal of H i . If x p0 hits O s , x p2 hits H + in Eq.26. They correspond to the spatial and the angular gradient of intensity I on the reflecting surface of the heliostat H + i [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]. Seeing Fig. 10a and Fig. 10b, these two gradients (spatial and angular) are both originated from the borders of another heliostat (from point x p1 in Fig. 17).

Contribution of backward-shadowing effect

H i R H i O s O l n i x p1 x p2 x p0 x p3 ω s ω p
Seeing the algorithm of ray-tracing in Fig17, S shad-b i can be estimated through the following integral formulation:

S shad-b i (π) = 4×n H k=1 l k p X ( x p1 )d x p1 Ωs p Ωs ( ω s )d ω s ( w 5 + w 6 )H( x p0 ∈ O s )H( x p2 ∈ H + i )H( x p3 ∈ R) , ( 48
)
where x p0 is the first intersection starting by x p1 , following -ω s , x p2 the first intersection starting by x p1 following ω s , x p3 the first intersection starting by x p2 following ω p , t i ,k the vector of circulation on the shadowing heliostat H i , l shad-b the length of the optical path between x p1 and x p2 . For each realization of estimating S blo i , we sample x p1 and ω s , and we calculate ω p . ω p is reflected following the normal n i which is the normal of H i . If x p0 hits O s , x p2 hits H + i and x p3 hits R we count w 5 + w 6 .

w 5 = sshad-b,spatial i ( x p1 , ω p , π) p X ( x p1 )p Ωs ( ω s ) ; ( 49
)
w 6 = sshad-b,angular i ( x p1 , ω p , π) p X ( x p1 )p Ωs ( ω s ) . ( 50
)
H i R H i O s O l n i x p1 x p2 x p0 x p3 ω s ω p
Similar to the backward-blocking effect, H i can be all other heliostats except H i itself in the heliostats field. In order to detect all backward-shadowing effects, we need to proceed with ray-tracing tests for all the borders of all heliostats in the field except the borders of H i (4 × (n H -1) in total where n H is the number of heliostats in the field). However, most of them will not create shadows on H i since they are too far from each other. In order to accelerate and optimize the algorithm, we will only proceed with the ray-tracing tests for the heliostats, of which the distance to H i is within l d . See also Algo.6 for the complete algorithm and Appendix.C for the explicit expression of w 5 and w 6 .

Results, validations, and comparisons

In this section, we will proceed with the following:

1. Validations of our method by finite difference method using four examples.

2. A comparison of the finite difference method with our method in calculating time.

Validations

We will here proceed with four simple examples (see Fig. 18) where we apply our method and finite difference method for validation purposes. The sensitivities of H 1 are calculated and results are shwon in Table .2, 3, 4 and 5. In case 1, only the spillage effect is detected. In case 2, the blocking effect is detected. In case 3, the backward-shadowing effect is detected. In case 4, the forward-shadowing effect is detected. The estimation results by our method are validated by the finite difference method, which will be discussed in the following subsection. Figure 18: The receiver R is a rectangular with lengths of 1.50 m and r 1 is its geometric center. The heliostats H 1 and H 2 are two rectangles with lengths of 1 m and h 1 , h 2 are their geometric centers respectively. In case 1, θ 1 is the horizontal angle of H 1 . In case 2, 3 and 4, the heliostat H 1 and H 2 point to r 1 according to the solar elevation angle.

e1 e2 e3 h 1 θ 1 H 1 r 1 R (a) Case 1: r1 = (0, 0, 0), h1 = (1, 5, 1) and θ1 = 30 • . Solar elevation angle = 60 • . e1 e2 e3 h 1 H 1 h 2 H 2 r 1 R (b) Case 2: r1 = (0, 0, 2), h1 = (0.25, 2, 0) and h2 = (0, 3, 0). Solar ele- vation angle = 60 • . e1 e2 e3 h 1 H 1 h 2 H 2 r 1

Approximation by finite difference method

The finite difference method is a standard method to approximate the sensitivities. When P is derivable around πi,j :

S i,j = lim →0 P (π i,j + ) -P (π i,j -) 2 . ( 51
)
The finite difference method approximates S i,j by giving a non-zero fixed value to . We note the result of approximation as Si,j : Table 3: Validation -Case 2: 10 6 realizations on each border and 10 9 realizations on the surface for an estimation of S 1,j and 10 9 realizations for each estimation of P when approximating S1,j by finite difference method.

S i,j ≈ P (π i,j + ∆π i,j ) -P (π i,j -∆π i,j ) 2∆π i,j ≡ Si,j . (52) unit W
Practically, the approximation of Si,j requires estimating P twice. Since the impacting power P is estimated by the Monte-Carlo method, its statistical standard deviation σ(P ) is also estimated at the same time [START_REF] William | Exploring monte carlo methods[END_REF]. When approximating Si,j , its standard deviation σ( Si,j ) can be obtained [START_REF] Gobet | Stochastic differential equations and Feynman-Kac formulas[END_REF]:

σ( Si,j ) ≈ σ(P (π i,j + ∆π i,j )) + σ(P (π i,j -∆π i,j )) 2∆π i,j . (53)
The major drawback of the finite difference method is that the related standard deviation σ( Si,j ) is hard to control [START_REF] Gobet | Stochastic differential equations and Feynman-Kac formulas[END_REF] and tends to infinity as ∆π i,j tends to zero.

Comparison with the finite difference method

In this section, we will apply and discuss the algorithms previously introduced in the context of a functioning solar tower station: Sierra SunTower [START_REF] Schell | Design and evaluation of esolar's heliostat fields[END_REF]. This case study has been chosen since our sensitivity model was developed for flat heliostats.

This solar tower station is located in Mojave Desert, California, at a latitude of 34.7 • . It consists of 2 solar towers, each with north and south heliostat subfields [START_REF] Schell | Design and evaluation of esolar's heliostat fields[END_REF].

The heliostat field and the two towers are symmetrically distributed. Therefore, we only focus on the field on the southwest side and its related tower, Zenith angle for the summer solstice at solar noon 79 • Table 6: Geometric configuration of the system [START_REF] Schell | Design and evaluation of esolar's heliostat fields[END_REF] The present work focuses on the optical performance (the impacting power P ). It is assumed that:

1. The sun's position is fixed (solar noon for the summer solstice) 2. Incoming Direct Normal Intensity (DNI) is considered homogeneous for the whole heliostat field (DNI = 1000W/m 2 ).

The incoming power, P (π) and its derivatives S(π) towards geometric parameters π ≡ [π i,j ] are estimated. It is reminded that i refers to the i th heliostat (i ∈ {1, 2 . . . , 6090}) and j to one of the six geometric parameters (j ∈ {1, 2, . . . , 6}) as shown on Fig. 4.

Results have already been shown at the beginning of this article (Fig. 2). It is noted that the accuracy of the simulation depends on the number of Monte Carlo realizations. Each Monte-Carlo realization consists of the corresponding sampling process and the ray-tracing tests described in Section.3.

When we focus on one parameter of one heliostat in the field (for example, the size of the 5297 th heliostat: π5297,6 ), it has a tinny impact on the total impacting power P . It causes a poor convergence performance when applying the finite difference method, as shown in Table .7. Compared to our method, the method of approximation by finite difference takes much more time to converge because the choice of ∆π enormously impacts the convergence rate and the precision for the approximation method of finite difference. In contrast, our method of estimation avoids this choice. These calculations are run parallelly in an ordinary laptop1 , and our method converges in 0.188s. According to [START_REF] William | Exploring monte carlo methods[END_REF], to reduce ten times the standard deviation of a Monte-Carlo estimation, approximately 100 times the calculation time is needed. In this application case, the approximation of the finite difference method will then spend approximately 2 × 10 11 more times of calculating time to attend the same order of standard deviation of our method.

Last but not least, it is found that in Table .2, Table .3, Table .4 and Table .5, the advantage of our method in convergence performance is not that significant compared to the finite difference method. We have the following comments:

1. The four cases are dedicated to validation purposes. Therefore, the geometries are relatively simple compared to a functioning SPT system. However, the finite difference method will become impracticable when the geometries become more complex.

2. Even with simple geometries, the standard derivation of the finite difference method (σ( Si,j )) is hard to be controlled. It is strongly dependent on the value of ∆π i,j chosen in Eq.52 and Eq.53.

3. The derivation related to the discretization can not be handled for the finite difference method. It is also largely related to the choice of ∆π in Eq.52 and Eq.53. However, our method can handle the statistical standard derivation. We avoid the derivation related to the discretization and the choice of ∆π i,j .

Finite difference method

Conclusion

We have presented general models for the intensity I and the sensitivity of intensity s i,j in a radiative system of an SPT system, as well as the complete algorithms for estimating the matrix of sensitivity of power S. These algorithms are validated by the finite difference method using four typical examples. Since the sensitivity of a heliostat of a geometric parameter S i,j is estimated by the accumulation of four contributions: S tar i,j , S blo i,j , S shad-b i,j and S shad-f i,j

, the contributions of different physical events to the sensitivity S i,j in an SPT system (blocking, shadowing, spillage) can be analyzed in detail. After that, an actual functioning SPT system was used as a test-case study. Firstly, a comparison of computing time was made for estimating sensitivity by our method and its approximation by the finite difference method. Secondly, the matrix of sensitivity of power S was calculated. Engineers can benefit from this information on sensitivities to optimize the optical design and the aiming strategy for the SPT system. Last but not least, this work makes the integration of gradient-based optimizations for the optical system in an SPT becomes possible. However, the sensitivity model developed in the present paper is only applied to a field of flat heliostats. Thus, a promising perspective would be to extend this model to curved heliostats.

This work received financial support from the from Region Occitanie and the SOLSTICE laboratory of Excellence (ANR-10-LABX-22-01).

A Boundary conditions of geometric sensitivities for a heliostat

In this work, a heliostat is modeled by two surfaces without thickness: the reflecting surface H + i and the rearward surface H - i . Based on the general model in [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], the boundary condition of geometric sensitivities for a cold, specular mirror with homogeneous reflectivity is developed, where I( x p , ω p , π) and I( x p , ω s , π) are spatially and angularly smooth:

s i,j ( x p , ω p , π) = -∂ 1, vi,j I( x p , ω p , π) + ρ∂ 1, vi,j I( x p , ω s , π) -ρ∂ a R I( x p , ω s , π) + ρ∂ a R I( x p , ω s , π) + ρs i,j ( x p , ω s , π) (54)
where ∂ 1, vi,j is the operator for spatial gradient following the vector v i,j which are the velocity of deformation (yielded in Appendix.D), ∂ a R the operator for angular gradient following the vector of rotation a, and a the reflected vector of a on H + i [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]. However, when we model the reflecting surface H + i , I( x p , ω p , π) is spatially discontinued because the heliostat has 4 edges and I( x p , ω s , π) can be angularly discontinued (when backward-shadowing occurs, see Fig. 10b).

In [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], the discontinuity of spatial gradient and angular gradient of intensity I are studied. The boundary condition of H + i can then be yielded. Furthermore, in this work, we classify the different sources of sensitivity based on the background of CSP applications:

s i,j ( x p , ω p , π) = star i,j ( x p , ω p , π) + sblo i,j ( x p , ω p , π) + sshad-b i,j ( x p , ω p , π) + ρs i,j ( x p , ω s , π). (
( x p , ω p , π) = -( ω p × v i,j ) • t i,k ρI( x p , ω s , π) (56)
and star,angular i,j

( x p , ω p , π) = -ρ∂ a R I( x p , ω s , π) + ρ∂ a R I( x p , ω s , π). ( 57
)
t i,k are the vectors of circulation of H + i , depending on x p , see Fig. 19. When the backward-shadowing does not occur, the angular gradient in Eq.57 will be the angular derivative of the boundary condition of O s (previously defined in Eq.15, which is the Pill-box sun shape model). The development is straightforward, and they are the sources on the border of the solar cone ∂Ω s :

∂ a R I( x p , ω s , π) = I 0 ( ω s × a) • ω c H(hit b ∈ O s ) (58) ∂ a R I( x p , ω s , π) = I 0 ( ω s × a ) • ω c H(hit b ∈ O s ) (59)
where H(hit b ∈ O s ) corresponds to a backward ray-tracing test which will be true if the direction ofω s starting from x p is not shadowed by other heliostats in the field. sblo i,j is the source on the border of heliostat ∂H + i :

sblo i,j ( x p , ω p , π) = ( ω p × v i,j ) • t i,k I( x p , ω p , π) (60)
When backward-shadowing effect occurs, H + i is shadowed by another heliostat H + i , a spatial gradient exists on the surface of H + i (Fig. 10a). An angular gradient exists within the solar cone Ω s (Fig. 10b). Furthermore, in our case, these sources on the surface of H + and within the cone Ω s all come from the border of the shadowing heliostat ∂H + i . Following the models of spatial gradient and angular gradient in [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], the source sshad-b i,j is on the borders of all heliostats in the field except the heliostat H i , which is ∂H + i , i = i. Similar to the source of targeting, we separate the source sshad-b i,j by sshad-b,spatial i,j and sshad-b,angular :

sshad-b i,j ( x p , ω p , π) = sshad-b,spatial i,j ( x p , ω p , π) + sshad-b,angular i,j ( x p , ω p , π) (61) where sshad-b,spatial i,j ( x p , ω p , π) = ρ v i,j • n i ω p • n i ( ω s × ω p ) • t i ,k I( x p , ω s , π) (62)
and sshad-b,angular i,j

( x p , ω p , π) = -ρl shad-b ( ω s × a × ω s ) • t i ,k I( x p , ω s , π) + ρl shad-b ( ω s × a × ω s ) • t i ,k I( x p , ω s , π). ( 63
)
t i ,k are the vectors of circulation of H + i [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], see Fig. 19 and l shad-b is the length of optical path between x p and x p (Fig. 17).

The boundary condition of the rearward surface H - i can be developed based on the general model of opaque-black surface in [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF]:

s( x b , ω s , π) = sshad-f i,j ( x b , ω s , π) (64)
and sshad-f i,j

is the source on the border of the heliostat

H - i : sshad-f i,j ( x b , ω b , π) = ( ω s × v i,j ) • t i,k I( x b , ω s , π) (65)

B Integral formulations

The vectors of contributions of S defined in Section.

where the Heaviside function H(hit f ∈ R) corresponds to a forward raytracing test which will be true if the source arrives at R (the ray staring at x p , following ω p impacts the receiver R).

S tar,angular i (π) is estimated by integrating star,angular i on H + i and on ∂Ω s :

S tar,angular i (π) = H + i d x p ∂Ωs d ω s H (hit f ∈ R) star,angular i ( x p , ω p , π) (67)
S blo i (π) is estimated by integrating sblo

Finally, the corresponding integral formulations of S tar (π), S tar (π), S blo (π), S shad-b (π) and S shad-f (π) are shown in Section.3.

C Explicit expression of Monte-Carlo weight C.1 Contribution of targeting

We substitute Eq.56 into Eq.42, as well as Eq.57 into Eq.43. The following explicit expressions of Monte-Carlo weight are then yielded:

w 1 = -l k DN I             ( ω p × v i,1 ) • t i,k ( ω p × v i,2 ) • t i,k ( ω p × v i,3 ) • t i,k ( ω p × v i,4 ) • t i,k ( ω p × v i,5 ) • t i,k ( ω p × v i,6 ) • t i,k             ;
(71)

w 2 = 2πI 0 S H + i | ω s • n i |          0 0 0 -2( ω s × a i,θ ) • ω c ω s × ( a i,φ -a i,φ ) • ω c 0          . ( 72
)
where S H + i the area of the reflecting surface of the heliostat H + i , l k the length of k th side of the heliostat H i and t i,k the vectors of circulation which is clockwise around the normal n i following the convention of [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], see Fig. 19. v i,j is the velocity of deformation linked to the point x p1 . Generally, it is a function of a point x on the heliostat H i : v i,j ≡ v i,j ( x).

( The vectors of v i,j are given by Table .8, where e = xx c,i and x c,i is the central position of the heliostat H i . The demonstrations are in Appendix.D and for example, the field of v i,4 ( x p ) and v i,5 ( x p ) are shown in Fig. 20 and Fig. 21.

v i,1 v i,2 v i,3 v i,4 v i,5 v i,6 [1, 0, 0] [0, 1, 0] [0, 0, 1] a i,θ × e a i,φ × e e l k
Table 8: Components of v i,j for a point x on H i In Eq.72, a i,θ is the axis of rotation related to πi,4 and a i,φ the axis of rotation related to πi,5 . a i,φ is the vector reflected by H + i from -a i,φ [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], see Fig. 19.

a i,θ v i,4 ( x p )

C.2 Contribution of backward-blocking

We substitute Eq.60 into Eq.45. The following explicit expression of Monte-Carlo weight is then yielded:

w 3 = l k DN I             ( ω p × v i,1 ) • t i,k ( ω p × v i,2 ) • t i,k ( ω p × v i,3 ) • t i,k ( ω p × v i,4 ) • t i,k ( ω p × v i,5 ) • t i,k ( ω p × v i,6 ) • t i,k             ; (74)

C.3 Contribution of forward-shadowing

We substitute Eq.65 into Eq.47. The following explicit expression of Monte-Carlo weight is then yielded:

w 4 = -l k DN I             ( ω s × v i,1 ) • t i,k ( ω s × v i,2 ) • t i,k ( ω s × v i,3 ) • t i,k ( ω s × v i,4 ) • t i,k ( ω s × v i,5 ) • t i,k ( ω s × v i,6 ) • t i,k             . (75)

C.4 Contribution of backward-shadowing

We substitute Eq.62 into Eq.49, as well as Eq.63 into Eq.50. The following explicit expressions of Monte-Carlo weight are then yielded: with

w 5 = ρl k DN I                        v i,1 • n i ω p • n i ( ω s × ω p ) • t i ,k v i,2 • n i ω p • n i ( ω s × ω p ) • t i ,k v i,3 • n i ω p • n i ( ω s × ω p ) • t i ,k v i,4 • n i ω p • n i ( ω s × ω p ) • t i ,k v i,5 • n i ω p • n i ( ω s × ω p ) • t i ,k v i,6 • n i ω p • n i ( ω s × ω p ) • t i ,k                        (76)
and

w 6 =ρl k DN Il shad-b         0 0 0 -2( ω s × a i,θ × ω s ) • t i ,k ( ω s × ( a i,φ -a i,φ ) × ω s ) • t i ,k 0         (77)

D Velocity of deformation

Following [START_REF] Lapeyre | A physical model and a monte carlo estimate for the specific intensity spatial derivative, angular derivative and geometric sensitivity[END_REF], we need to define a geometric space, and a material space for a heliostat H i , and special notations are needed to calculate the velocity of deformation. In geometric space, the positions y and directions ω are functions of the geometric parameter πi,j , while in material space, the positions are noted y and directions are noted ω. They are not depended on πi,j . The function Z and Ω link the two spaces so that the positions and directions in one space can refer to the other space: y = Z( y, πi,j ); ω = Ω( ω, π) (78)

We define the velocity of deformation v i,j as the derivative of Z with respect to πi,j : v i,j ≡ ∂ πi,j Z( y, πi,j )

For different πi,j , the velocities of deformations v i,j are summarized in Table .9, where a θ and a φ are two axes of rotation for the heliostat H i (in Fig. 4) and l k the original length of size of the heliostat.

  Sensitivity of azimuth angle for each heliostat.

Figure 2 :

 2 Figure 2: Sensitivity of the total impacting power P [w] for each heliostat in the field. Each point is dedicated to a heliostat, pointing to the center of the receiver at the moment of solar noon at the summer solstice. The layout corresponds to a real heliostat field: Sierra SunTower [34] (receiver at 50 m height at (0,0) and more details of the heliostats field are in Table.6 and Section.4.2). The corresponding color indicates its sensitivity with respect to its positions, rotations, and lengths. (10 6 Monte-Carlo realizations on each border of each heliostat and 10 8 Monte-Carlo realizations on each surface of heliostats).

  Figure 2: Sensitivity of the total impacting power P [w] for each heliostat in the field. Each point is dedicated to a heliostat, pointing to the center of the receiver at the moment of solar noon at the summer solstice. The layout corresponds to a real heliostat field: Sierra SunTower [34] (receiver at 50 m height at (0,0) and more details of the heliostats field are in Table.6 and Section.4.2). The corresponding color indicates its sensitivity with respect to its positions, rotations, and lengths. (10 6 Monte-Carlo realizations on each border of each heliostat and 10 8 Monte-Carlo realizations on each surface of heliostats).

Figure 3 :

 3 Figure 3: An SPT system consists of heliostats H i (with an index of i) on the ground and a receiver R at the top of a tower.

Figure 6 :

 6 Figure 6: Orientation convention of the heliostat H i

Figure 8 :

 8 Figure 8: The optical path of I is in yellow, and that of s i,j is in red. When all the sources of star,spatial i,1

i, 1

 1 is partially captured by R and S tar,spatial i,1 = 0.

Figure 9 : 1 = 0 .

 910 Figure9: With a perturbation of πi,4 , the perturbation of I is observed on all surfaces of H i and borders of Ω s (the solar cone). When all the sources of star,angular i,1 is all captured on R i,1 , the perturbation of π causes the perturbation of P because of the cosine effect and S tar,angular i,1 = 0. When all the sources are missing, S tar,angular i,1 = 0. Finally, when the spillage effect, forward-blocking effect or backward-shadowing effect occurs, star,angular i,1

  Perturbation due to the change of normal

Figure 10 :

 10 Figure 10: When the backward shadowing effect occurs, with a perturbation of πi,4 , P is perturbed because of its change of normal. The source of sensitivity sshad-b i,4is on the border of the neighboring heliostat. P is also perturbed because of its change of position. The source of sensitivity sshad-b i,1

Figure 14 :

 14 Figure 14: Algorithm to compute S tar i .x p1 is on the 4 borders of the heliostat ∂H i having k as the index of the border, x p0 the first intersection starting by x p1 following -ω s , x p2 the first intersection starting by x p1 following ω p , ω s the direction on the border of solar cone ∂Ω s , ω p the direction of reflection of ω s , x p0 the first intersection starting by x p1 following -ω s and x p2 the first intersection starting by x p1 following ω p . For each realization of estimating S tar,spatial

Figure 15 :

 15 Figure 15: Algorithm to compute S blo i .x p1 is the first intersection starting by x p2 , following -ω p , x p0 the first intersection starting by x p1 following -ω s , x p3 the first intersection starting by x p1 following ω p . For each realization of estimating S blo i , we sample x p2 and ω s , and we calculate ω p . If x p0 hits O s , x p1 hits H i and x p2 hits R we count w 3 .

  of S shad-b i , estimated by the sum of the sources of backward-shadowing sshad-b i,j

Figure 16 :

 16 Figure 16: Algorithm to compute S shad-f i .x p0 is the first intersection starting by x p1 , following -ω s , x p2 the first intersection starting by x p1 following ω s , x p3 the first intersection starting by x p2 following ω p . For each realization of estimating S blo i , we sample x p1 and ω s , and we calculate ω p . ω p is reflected following the normal n i which is the normal of H i . If x p0 hits O s , x p2 hits H +

i

  and x p3 hits R we count w 4 . of position and the change of normal, impacting the backward-shadowing effect. Therefore, the contribution S shad-b i can be divided by two parts: S shad-b,spatial i for the change of position and S shad-b,angular i for the change of normal. Also, the sources sshad-b,spatial i,j are divided by two parts respectively: sshad-b,spatial i,j and sshad-b,angular i,j

Figure 17 :

 17 Figure 17: Algorithm to compute S shad-b i .x p0 is the first intersection starting by x p1 , following -ω s , x p2 the first intersection starting by x p1 following ω s , x p3 the first intersection starting by x p2 following ω p , t i ,k the vector of circulation on the shadowing heliostat H i , l shad-b the length of the optical path between x p1 and x p2 . For each realization of estimating S blo i , we sample x p1 and ω s , and we calculate ω p . ω p is reflected following the normal n i which is the normal of H i . If x p0 hits O s , x p2 hits H + i and x p3 hits R we count w 5 + w 6 .

R 1 R

 1 (c) Case 3: r1 = (0, 0, 2), h1 = (0.25, 2, 0) and h2 = (0, 3, 0). Solar elevation angle = 30 • . (d) Case 4: r1 = (0, 0, 2), h2 = (0.25, 2, 0) and h1 = (0, 3, 0). Solar elevation angle = 30 • .

HHH

  (hit f ∈ R) sblo i ( x p , ω p , π)(68)S shad-b i (π) is estimated by integrating sshad-b i on all the borders of other heliostats in the field ∂H + i and on Ω s :(hit f ∈ R) sshad-b i ( x p , ω p , π) (hit f ∈ R) sshad-f i ( x b , ω b , π)

Figure 20 :Figure 21 :

 2021 Figure 20: Schema of the field of v i,4 ( x p )

  y + πi,4 ( a i,θ × y) a i,θ × y πi,5 y + πi,5 ( a i,φ × y) a i,φ × y

Table 2 :

 2 Validation -Case 1: 10 6 realizations on each border and 10 8 realizations on the surface for an estimation of S 1,j and 10 7 realizations for each estimation of P when approximating S1,j by finite difference method.

	unit	W/m	W/m	W/m	W/rad	W/rad	W/m
	j	1	2	3	4	5	6
	S tar,spatial 1,j S tar,angular 1,j S blo 1,j S shad-b,spatial 1,j S shad-b,angular 1,j S shad-f 1,j	-793.35 0.00 139.39 0.00 0.00 0.00	143.71 × 10 -3 0.00 -277.29 0.00 0.00 0.00	141.22 × 10 -3 456.62 0.00 -912.74 -415.94 -188.09 0.00 0.00 0.00 0.00 0.00 0.00	356.02 × 10 -3 -463.55 14.84 0.00 0.00 0.00	991.76 0.00 -234.30 0.00 0.00 0.00
	S1,j σ S1,j S1,j σ S1,j	-653.96 270.41 × 10 -3 -654.25 180.35 × 10 -3	-277.15 890.09 × 10 -3 -277.20 181.87 × 10 -3	-415.80 1.03 -415.90 181.38 × 10 -3	-644.22 13.81 -640.00 18.23	-448.36 2.08 -450.50 1.82	757.46 643.34 × 10 -3 757.50 1.82

Table 4 :

 4 Validation -Case 3: 10 6 realizations on each border and 10 9 realizations on the surface for an estimation of S 1,j and 10 9 realizations for each estimation of P when approximating S1,j by finite difference method.

	unit	W/m	W/m	W/m	W/rad	W/rad	W/m
	j	1	2	3	4	5	6
	S tar,spatial 1,j S tar,angular 1,j S blo 1,j S shad-b,spatial 1,j S shad-b,angular 1,j S shad-f 1,j S1,j σ S1,j S1,j σ S1,j	-1.29 × 10 -3 0.00 0.00 0.00 0.00 -68.35 -68.35 181.08 × 10 -3 -68.50 771.70 × 10 -3	-2.10 × 10 -3 0.00 0.00 0.00 0.00 375.08 375.08 220.47 × 10 -3 375.00 771.28 × 10 -3	-2.31 × 10 -3 0.00 0.00 0.00 0.00 -649.65 -649.66 382.90 × 10 -3 -649.50 770.40 × 10 -3	-849.47 1.69 × 10 3 0.00 0.00 0.00 -318.62 523.12 10.45 530.00 7.72	242.46 × 10 -3 -88.94 × 10 -3 0.00 0.00 0.00 17.21 17.37 973.89 × 10 -3 15.00 7.72	1.06 × 10 3 0.00 0.00 0.00 0.00 -232.07 823.17 206.71 × 10 -3 823.00 771.70 × 10 -3

Table 5 :

 5 Validation -Case 4: 10 6 realizations on each border and 10 9 realizations on the surface for an estimation of S 1,j and 10 9 realizations for each estimation of P when approximating S1,j by finite difference method. which includes 6090 flat heliostats. It is easy to extend the calculation for the whole field in the south (including 12180 heliostats). However, extra efforts are needed for the field on the north side since the towers' shadows will impact the sensitivities.Table.6 shows the geometric parameters of the radiative system, as well as the solar positions for the summer solstice.

	Characteristic	Value
	Overall north-south length	95 [m]
	Overall east-west length	175 [m]
	Distance between the tower and the first row of heliostats	12.5 [m]
	Number of heliostats per column	58
	Number of heliostats per row	105
	Size of heliostats Height of the centre of receiver	1.067 × 1.067 [m 2 ] 50 [m]
	Size of receiver	4 × 4 [m 2 ]

Table 7 :

 7 π5297,6 is originally equal to 1.067 (see Table.6), while herein, a reasonable value of ∆π 5297,6 is chosen: ∆π 5297,6 = 0.1[m]. Two estimations of impacting power (P (π 5297,6 -∆π 5297,6 ) and P (π 5297,6 + ∆π 5297,6 )) are realized following Algo.2 (10 9 Monte-Carlo realizations) . The sensitivity of power S5297,6 is approximated based on these two estimations. Also, the sensitivity of power S 5297,6 is estimated following Algo.1 (10 4 Monte-Carlo realizations for each border of the heliostats).

  3.2 are the sum of the vectors of sources captured by the receiver R respectively. S tar,spatial

	by integrating star,spatial i	on ∂H + i and on Ω s :	i	(π) is estimated
	S tar,spatial i	(π) =	i ∂H +	d x p	Ωs	d ω s

H (hit f ∈ R) star,spatial i ( x p , ω p , π)

  Figure 19: t i,k are the clockwise vectors of circulation around n i

		a i,φ t i,4	n i	a i,φ
		H i
	a i,θ a i,θ	t i,3		t i,1
			t i,2
				73)

Table 9 :

 9 The functions Z and v for the 6 geometric parameters of heliostat Geometry of heliostats and receiver, ω c , N N number of realizations Initialization: n ← 0, sum ← 0, sum2 ← 0 while n < N do Sample x p1 on S H + based on p X Sample ω s within Ω s based on p Ωs . Compute ω p by the law of specular reflection Get x p0 from x p1 followingω s Get x p2 from x p1 following ω p if x p0 ∈ O s & x p2 ∈ R then

	E Algorithms
	Algorithm 2 Estimate P
	Input: Compute	ŵ
	ŵ sum2 ← sum2 + ŵ2 sum ← sum + else
	sum ← sum + 0 sum2 ← sum2 + 0 end if
	n ← n + 1 end while
	Output: P ← sum N , σ(P ) ←	sum2 N -( sum N ) 2 N -1

4 CPUs of i5 Intel™, 8th generation

Algorithm 3 Estimate S tar i and σ( S tar i ) Input: Geometry of heliostats and receiver, ω c , N N number of realizations Initialization: S tar i ← 0, σ( S tar i ) ← 0 while k < 4 do H i has 4 borders Initialization: n ← 0, --→ sum ← 0, ---→ sum2 ← 0 while n < N do Sample x p1 on l k based on p X Sample ω s within Ω s based on p Ωs .

Compute ω p based on ω s and n i by the law of specular reflection Get

Ignore the too-far-away heliostats for backward-blocking effect Initialization:

Compute ω p based on ω s and n i by the law of specular reflection Get

Algorithm 5 Estimate S shad-f i and σ( S shad-f i )

Input: Geometry of heliostats and receiver, ω c , N N number of realizations Initialization:

Compute ω p based on ω s and n i by the law of specular reflection Get