Pedro Granate

Maxime Pedro J Coelho
email: pedro.coelho@tecnico.ulisboa.pt

Maxime Roger

IDMEC Laeta

RADIATIVE HEAT TRANSFER IN STRONGLY FORWARD SCATTERING MEDIA USING THE DISCRETE ORDINATES METHOD

Keywords: Discrete ordinates method, Finite volume method, Spherical harmonics discrete ordinates method, Forward scattering media, Phase function normalization

Radiative heat transfer in strongly forward scattering media using the discrete ordinates method

INTRODUCTION

The discrete ordinates method (DOM) was initially developed by Chandrasekhar [START_REF] Chandrasekhar | Radiative transfer[END_REF] in his work about interstellar radiation. A few years later, Carlson and Lathrop [START_REF] Carlson | Transport theory -The method of discrete ordinates[END_REF] applied the method to the neutron transport theory. There were some early attempts to use the method to solve thermal radiation problems, but it did not become popular until the last quarter of the 20th century, with the pioneering work of Fiveland [START_REF] Fiveland | Discrete-ordinates solutions of the radiative transport equation for rectangular enclosures[END_REF]. Subsequent developments include improved spatial and angular discretization schemes, generalization to more complex grid structures and to complex enclosures, extensions to non-grey media, media with variable refractive index, and transient problems, as well as parallel implementation [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF]. Although the accuracy achieved using the DOM is often satisfactory and the computational requirements moderate, it is well known that errors due to ray effects and false scattering are generally present [57]. Another source of error arises in strongly anisotropic media, as a result of lack of conservation of scattered energy [89] and/or asymmetry factor of the phase function [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF]. This is the subject of the present work, which is concentrated on scalar radiative transfer.

Although the focus is on thermal radiation the methods discussed here may be applied to other fields, for example atmospheric radiation, astrophysics, neutron transport, etc.

When the DOM is used, the in-scattering term of the radiative transfer equation (RTE) is approximated by a quadrature over all discrete directions, and there is no guarantee that the integral of the scattering phase function over a spherical surface of unity radius yields 4π, as it should if there were no approximation. A simple normalization of the scattering phase function to enforce conservation of the scattering energy has been proposed a long time ago by Kim and Lee [START_REF] Kim | Effect of anisotropic scattering on radiative heat transfer in twodimensional rectangular enclosures[END_REF] and implemented by Liu et al. [START_REF] Liu | On the discrete ordinates method for radiative heat transfer in anisotropically scattering media[END_REF]. However, it was found that this method only yields accurate results for small or moderate asymmetry factors of the phase function [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF]. When the asymmetry factor increases, the normalization method changes the overall shape of the scattering phase function. As a consequence, even though the scattered energy is conserved, the asymmetry factor is not, and this may yield large errors in the case of highly anisotropic phase functions, which become larger in the case of optically thick media.

Boulet et al. [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF] found that this problem is largely overcome using the finite volume method (FVM) instead of the DOM. The effectiveness of the FVM to surmount this problem has also been found in problems with collimated irradiation [START_REF] Collin | Acute anisotropic scattering in a medium under collimated irradiation[END_REF].

A different normalization procedure, which conserves both the scattered energy and the asymmetry factor for the Henyey-Greenstein phase function, was proposed by Hunter and Guo [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF]. This method produced results that closely agree with FVM predictions in an axisymmetric cylindrical enclosure [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF], and Monte Carlo benchmark solutions in a cubic enclosure [START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part II: Benchmark comparisons[END_REF]. The application to a Legendre scattering phase function in ultrafast radiative transfer is reported in [START_REF] Hunter | Phase-function normalization for accurate analysis of ultrafast collimated radiative transfer[END_REF]. A drawback of this normalization procedure is the need to predetermine a normalization matrix. This may be avoided using a simpler normalization technique [START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF], which simultaneously conserves the scattered energy and the asymmetry factor of the phase function, and maintains most of the phase-function shape. The method was applied to the Henyey-Greenstein phase function, yielding results similar to those obtained using the previous normalization technique, but with a lower computational effort. This technique has recently been extended to collimated radiation for Henyey-Greenstein or Legendre scattering phase-functions [START_REF] Hunter | Applicability of phase-function normalization techniques for radiation transfer computation[END_REF].

It is also possible to approximate a strongly forward scattering phase function by another approximate phase function where the forward scattering peak is dealt with a Dirac delta function, and the remaining scattering contribution is accounted for using an isotropic or a moderate anisotropic phase function. In such a case, the lack of conservation of scattered energy or asymmetry factor is largely overcome.

In the present work, the normalization methods of Liu et al. [START_REF] Liu | On the discrete ordinates method for radiative heat transfer in anisotropically scattering media[END_REF] and Hunter and Guo [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF][START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF] are applied to three test problems along with two other methods that do not rely on the normalization of the phase function. One of these methods is the FVM, which is applied here only to the in-scattering term. This means that the DOM is still employed, except in the evaluation of that term, in contrast to previous works where the angular discretization of the RTE was fully carried out using the FVM. In this sense, the present approach is a DOM-FVM hybrid formulation. The other method is based on the spherical harmonics discrete ordinates method, which was proposed by Evans in the framework of atmospheric radiation [START_REF] Evans | The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer[END_REF]. We are not aware of previous application of this method to thermal radiation problems. The approximation of a strongly scattering forward phase function by a Dirac delta function plus a smooth scattering contribution will be also investigated.

MATHEMATICAL FORMULATION

The transient RTE for an emitting-absorbing-scattering grey medium may be written as follows [START_REF] Modest | Radiative heat transfer[END_REF]:

                           d , t , , I I t , , I t , , I t t , , I c s b s s s r r s r s r s s r      4 4 1 (1)
where   t I , , s r is the radiation intensity in direction s, r is the position vector, t is the time, c is the speed of light in vacuum, I b is the blackbody radiation intensity, ,  and  s are the absorption, extinction and scattering coefficients of the medium, respectively, and   s s ,   is the scattering phase function. he boundary condition for a grey surface that emits and reflects diffusely is given by [START_REF] Modest | Radiative heat transfer[END_REF]: are the radiation intensities at boundary point r w that leave the boundary along s direction and arrive along s  direction, respectively, I bw is the blackbody radiation intensity at the temperature of the boundary surface,  is the surface emissivity,  is the surface reflectivity, and n is the unit vector normal to the surface and pointing into the medium. When collimated radiation is present, the RTE is decomposed into a diffuse and a collimated component, as described in Modest [START_REF] Modest | Radiative heat transfer[END_REF].

                 d t
In the present work, the DOM was used to solve equation ( 1). The spatial discretization was performed using the CLAM scheme, and the angular discretization was carried out using the S N quadrature, except in the cases identified below. In the case of transient problems, the time discretization was achieved using a second-order Runge-Kutta method. Details on the discretization procedure may be found elsewhere [START_REF] Coelho | Discrete ordinates and finite volume methods: mathematical formulation[END_REF]. When the DOM is employed, the angular integrals are approximated by quadratures, e.g., the in-scattering term is evaluated as

        j N j j i j w t I d t I d           1 4 , , , 4 1 , , , 4 1 s s
s r s s s r    for i = 1, 2, …, N d (3)
where N d is the number of discrete directions and w j is the quadrature weight for the jth direction. The scattering phase function is defined such that the integral of the phase function is equal to 4, but the approximation introduced by the quadrature generally yields an error, i.e., the scattered energy is not exactly conserved. The methods used in the present work to address this problem are described below.

Normalization of Kim and Lee [9]

The technique proposed by Kim and Lee [START_REF] Kim | Effect of anisotropic scattering on radiative heat transfer in twodimensional rectangular enclosures[END_REF], which will be referred hereafter as KL, is straightforward, and consists in multiplying the phase function by a normalization coefficient, so that the new phase function,  ~, is given by:

      1 1 4 1                 j N j j i j i j i w , , , ~d s s s s s s  (4)
The scattering phase function  in the in-scattering term of equation ( 1) is replaced by the normalized one,  ~, and this guarantees that the discrete form of the integral on the left side of equation ( 3) is indeed equal to 4when the radiation intensity is equal to unity for all directions.

Normalization of Hunter and Guo[10]

Boulet et al. [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF] noticed that the KL normalization procedure may lead to strong modifications in the values of the discretized phase function when the asymmetry factor is large, and this compromises the accuracy of the results. Hunter and Guo [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF] suggested that this may be avoided by enforcing conservation of the asymmetry factor of the phase function, g, in addition to conservation of the scattered energy. This is achieved by using a normalized phase function,  ~, such that the following conditions hold:

  1 4 1 1     j N j j i w , d s s  for i = 1, 2, …, N d (5)     g w , ~j j i N j j i d      s s s s 1 4 1  for i = 1, 2, …, N d (6)
The normalized phase function is defined as follows:

        j i ij i j j i , A , , ~s s s s s s       1 (7)
There are 2N d equations (equations 5 and 6) to be satisfied and N d (N d +1)/2 unknows (equation 7), so that the system of equations has an infinite number of solutions for the normalization coefficients A ij . These coefficients were determined using the least norm solution via QR factorization [START_REF] Golub | Matrix computations[END_REF], where R is an upper triangular nonsingular

(2N d × 2N d ) matrix and Q is an orthogonal N d (N d +1)/2 × 2N d matrix, i.e., Q T Q is equal to the identity matrix of rank 2 N d .
Hereafter, this normalization procedure will be referred to as HG1.

Normalization of Hunter and Guo[16]

A simpler normalization procedure has been proposed by Hunter and Guo [START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF], which will be referred to as HG2. In this method, equations ( 5) and ( 6) are still enforced, but only the forward and the backward terms of the scattering phase function are normalized, while the others remain unchanged. Hence, equation ( 7) is only used when s i .s j = 1 or s i .s j = -1. This leads to a number of equations equal to the number of unknown A ij coefficients, which may be easily solved analytically after inserting equation (7) into equations ( 5) and ( 6), yielding explicit expressions for those coefficients (see [START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF] for details).

Finite volume method

An alternative to the previous normalization techniques of the scattering phase function is to replace the quadrature used in the DOM for the calculation of the integral on the left side of equation ( 3) by the integration procedure used in the finite volume method, which may be expressed as follows:

    j N j ij j d I d , t , , I           1 4 4 1 4 1    s s s r for i = 1, 2, …, N d (8a)
with ij  defined as follows:

  j i ij i j d d              s s , (8b)
The integrals on the right hand side may be evaluated analytically if the scattering phase function is simple enough. In the general case, however, a numerical integration is required, causing a numerical integration error, and preventing the scattering energy from being exactly conserved. This error may be greatly reduced by splitting the control angles  i and  j into smaller angles. In the present work, an adaptive integration method [START_REF] Patterson | The optimum addition of points to quadrature formulae[END_REF] was used for the quadruple integration in equation 8(b), with a prescribed relative error of 10 -4 . This guarantees that the integration error is also negligible for the vast majority of problems. In order to apply this method, the boundaries of the control angles must be defined, and therefore the S N quadrature cannot be used. Instead, the polar /azimuthal discretization commonly used in the FVM has been employed.

Spherical harmonics discrete ordinates method

The last method used in the present work to evaluate the in-scattering term of the RTE is based on the spherical harmonics discrete ordinates method introduced by Evans [START_REF] Evans | The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer[END_REF] in the field of atmospheric radiation. It will be referred to below as SHDOM. This method is rather different from the standard DOM, but the only feature of the SHDOM that we will retain here is in the evaluation of the in-scattering term, while the standard DOM is used elsewhere. This term is evaluated using spherical harmonics. Denoting the in-scattering term of the RTE by S SH (r,s), where the temporal dependence was omitted for simplicity, the following truncated expansion in spherical harmonics is taken:

          r s s s s r s r m l L l M M m m l s S Y d , , I , S              0 4 SH 4    (9) In this equation,   s m l Y are spherical harmonics,   r m l S
are coefficients that depend only on the spatial coordinates, and M ≤ L. Note that terms such that M > l are equal to zero. The real spherical harmonics are defined as follows:

                      m u P ! m l ! m l l m u Y m l mo m l m l       
        r s r m l L l M M m k j m l jk S m u , S        0 DOM   (11)
where  j =cos j and subscripts j and k identify one among the N d directions of the DOM angular discretization. The RTE is solved using the in-scattering term given by equation [START_REF] Liu | On the discrete ordinates method for radiative heat transfer in anisotropically scattering media[END_REF] to obtain the radiation intensity for every direction, namely I jk (r)=I(r,s jk ).

The radiation intensity at every grid node, evaluated from the DOM, is transformed to spherical harmonics space as follows:

      r s s r m l m l L l M M m I Y , I       0 SH (12)
The coefficients   r m l I are given by [START_REF] Evans | The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer[END_REF]:

                  jk k j m l k N j N k j m l m l m l , I m u d , I m u d , I Y I s r s r s r s r                        1 1 4 4 (13)
where N  and N  are the number of cosine polar angles and the number of azimuthal angles, respectively, and  j and  k are the quadrature weights. The values of N  and N  are related to L and M according to the following relations [START_REF] Evans | The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer[END_REF], which ensure that the spherical harmonics are orthogonal when their product is integrated:

1    N L (14) 1 2    N M (15)
In addition, since M must be lower or equal than L, the following relation also holds:

N  ≤ 2N  .
A Gauss-Legendre quadrature was used for the integration in the polar angle, i.e., the cosines of the polar angles of the discrete directions in the DOM,  j , correspond to the roots of the Legendre polynomial of order equal to that of the quadrature, i.e., N  . The azimuthal angles are evenly spaced, and the quadrature weights  k are equal to 2/N  .

The scattering phase function written in terms of Legendre polynomials may also be expressed in terms of truncated spherical harmonics:

        s s               m l l l m m l L l l L l l l Y Y l P 0 0 1 2 4 cos cos    (16)
where

s s     cos
and P l is the Legendre polynomial of order l. The Henyey-Greenstein phase function may be written in terms of Legendre polynomials by setting  l =(2l+1) g l . The in-scattering term is evaluated in spherical harmonics space as follows:

                  r s s s r s s s s r s r m l m l l L l M M m s j i L i i i j j i i m l m l L l M M m s s I Y l d Y Y i I ' Y d , , I , S 1 2 1 2 4 4 4 0 0 4 0 4 SH                                                           (17)
where the orthogonality property of the spherical harmonics has been used. It follows from equations ( 9) and ( 17) that the coefficients of the spherical harmonics expansion of the inscattering term are given by:

    r r m l l s m l I l S 1 2     (18)
The solution algorithm may be summarized as follows:

(i)
The in-scattering source term evaluated in the spherical harmonics space is transformed to discrete ordinates according to equation [START_REF] Liu | On the discrete ordinates method for radiative heat transfer in anisotropically scattering media[END_REF]. This step is ignored in the first iteration.

(ii)

The RTE is solved using the DOM, with the in-scattering term evaluated in step (i).

In the first iteration, the in-scattering term is calculated using a quadrature, as in the standard DOM.

(iii)

The radiation intensity field calculated in step (ii) is transformed to spherical harmonics according to equation [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF], with the coefficients given by equation [START_REF] Collin | Acute anisotropic scattering in a medium under collimated irradiation[END_REF]. Only these coefficients need to be calculated.

(iv)

The in-scattering term in spherical harmonics space is calculated from the radiation intensity field in spherical harmonics using equation [START_REF] Hunter | Applicability of phase-function normalization techniques for radiation transfer computation[END_REF]. Only the coefficients

  r m l S
, given by equation ( 18), of the spherical harmonics expansion of the in-scattering term, which is expressed by equation ( 9), need to be evaluated. This involves simple multiplications rather than summations that typically appear in quadrature rules.

Approximate scattering phase function

In addition to the methods described above, which are applied to the Henyey-Greenstein scattering phase function of the problem under consideration, we will further examine the approximation of the Henyey-Greenstein phase function using the transport, the delta-Eddington or the delta-M approximations. In all cases, the approximate scattering phase function is defined as follows:

                 cos 1 cos 1 2 cos * f f  (19)
In the transport approximation [START_REF] Davison | Neutron transport theory[END_REF], the phase function * is equal to unity and the parameter f is set equal to the asymmetry factor. If the transient term of RTE is neglected, then the transport approximation is equivalent to replace the anisotropic scattering problem by an isotropic one with the following modified scaled parameters:

  f f ˆ      1 1 (20a)     f ˆ  1 (20b)
In these equations,  and  stand for the albedo and optical thickness, respectively, and the circumflex denotes a scaled parameter. In the delta-Eddington approximation, the phase function * is defined as follows [START_REF] Joseph | The delta-Eddington approximation for radiative flux transfer[END_REF]:

        cos 3 1 cos g * ( 21
)
where parameters f and g´ are given by:

2 g f  (22a) f f g g     1 (22b)
In both cases, the DOM may be used without normalization of the approximated phase function, since the strong forward scattering peak is fully accounted for by the Dirac delta function, and the remaining term of the phase function is smooth, no matter the value of the asymmetry factor.

The delta-M phase function is defined as follows [START_REF] Wiscombe | The delta-M method. Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions[END_REF]:

               cos P f f g n n n M n * 1 1 2 cos 0 (23)
where parameter f is given by f = g 2M . In the case of M=1, the delta-Eddington approximation is recovered. When M increases, the delta-M phase function becomes closer to the Henyey-Greenstein phase function and the contributions of the first and second terms on the right side of equation ( 19) are smaller and larger, respectively. Therefore, the asymmetry factor of * increases. If M is small enough, the normalization of the phase function is not needed.

Otherwise, it becomes necessary. However, in such a case there is no advantage in the replacement of the Henyey-Greenstein by the delta-M phase function.

RESULTS AND DISCUSSION

Three test cases are presented in this section. The first two are steady-state radiative transfer problems, one in a 3D cubical domain and the other in a 2D axisymmetric enclosure. The last test case addresses a transient problem with collimated incident radiation. The methods outlined above are compared.

Their acronyms are summarized here to allow easier identification: MC  Monte Carlo method DOM  Discrete ordinates method KL  DOM using Kim and Lee [START_REF] Kim | Effect of anisotropic scattering on radiative heat transfer in twodimensional rectangular enclosures[END_REF] normalization method HG1  DOM using Hunter and Guo [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF] normalization method HG2  DOM using Hunter and Guo [START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF] normalization method FVM  DOM using the finite volume method to evaluate the in-scattering term of the RTE SHDOM -Spherical harmonics discrete ordinates method

Test case 1

The first test case consists of a cube of unit side length, L, containing a grey and cold (nonemitting) medium with an optical thickness =L=10. The walls are black and cold, except the bottom one (z=0), which has a unit emissive power. The scattering albedo is = s /=1, and the Henyey-Greenstein phase function was considered. The calculations were carried out using a uniform mesh with 25×25×25 grid nodes and an S 12 quadrature for the DOM and the KL, HG1 and HG2 normalization methods. A quadrature with 7 polar angles and 3 azimuthal angles per octant was used in the case of the FVM, yielding the same number of discrete directions as the S 12 . In the case of the SHDOM, N  and N  were set to 14 and 12, respectively, yielding again the same number of discrete directions. This test case was also used by Boulet et al. [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF] and Hunter and Guo [START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part II: Benchmark comparisons[END_REF][START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF]. Boulet et al.

[12] used the FVM for the spatial and angular discretization of the RTE, along with the STEP scheme, a slightly coarser mesh (21×21×21) and a little finer angular discretization with 4 polar angles and 6 azimuthal angles per octant. Hunter and Guo [START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part II: Benchmark comparisons[END_REF] employed the DOM with the same spatial and angular discretizations used in the present work, but employed the positive scheme instead of the CLAM scheme. Our results are very similar to those obtained in [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF] and [START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part II: Benchmark comparisons[END_REF]. The results presented below are compared with a Monte Carlo (MC) solution also reported in [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF].

Figure 1 shows the influence of the grid size and quadrature on the incident heat flux, normalized by the emissive power of the hot boundary, along a central line on the top boundary of the cubical enclosure. The plotted results were calculated using the DOM with the HG2 normalization method of the phase function, for an asymmetry factor equal to 0.93.

It can be seen that the grid size has a small influence on the results. The predictions obtained using the S 4 quadrature are rather different from all the others, but the S 12 and S [START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF] give results similar to each other, while the results of the S 8 quadrature depart a little from them.

Therefore, the mesh with 25×25×25 grid nodes and the S 12 quadrature are adequate for the present test case.

A comparison of the various methods used to handle the strong anisotropy of the phase function is displayed in figure 2. When the asymmetry factor is equal to 0.20, all methods yield similar results, which closely match the Monte Carlo reference solution. In fact, normalization of the scattering phase function is not needed in this case. If the asymmetry factor is equal to 0.80, the KL method overestimates by up to 7% the predictions of all other methods, which are in good agreement with each other and with the Monte Carlo reference solution. In the case of an asymmetry factor of 0.93, the KL results are rather poor. This is consistent with previous works [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF][START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part II: Benchmark comparisons[END_REF], and is explained by the change of shape of the scattering phase function, due to the normalization introduced to enforce conservation of the scattered energy. The asymmetry factor is modified due to this normalization, leading to the poor predictions observed. This problem is enhanced with the increase of the asymmetry factor. The remaining four methods employed here successfully overcome the limitations of the KL method, even though their predictions are not as close to the MC solution as for the lower values of the asymmetry factor. The HG1 and HG2 normalization methods yield very similar results, which underpredict by 8 to 10% the MC benchmark solution, in agreement with the results reported in the literature [START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part II: Benchmark comparisons[END_REF][START_REF] Hunter | A new and simple technique to normalize the HG phase function for conserving scattered energy and asymmetry factor[END_REF]. Both the SHDOM and the FVM perform very satisfactorily, and their predictions are even closer to the MC solution than those of the HG1 and HG2 methods. The very good performance of the FVM is consistent with the findings of Boulet et al. [START_REF] Boulet | On the finite volume method and the discrete ordinates method regarding radiative heat transfer in acute forward anisotropic scattering media[END_REF].

In general, the computational requirements of the SHDOM are the lowest, even though the need to switch between the discrete ordinates and the spherical harmonics spaces might suggest the opposite. Note, however, that the calculation of the in-scattering term in the spherical harmonics space only requires the evaluation of the coefficients of the spherical harmonics expansion by means of equation [START_REF] Evans | The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer[END_REF], which involves a straightforward multiplication. In the case of g=0.20, the CPU time for SHDOM was about 600 s in a personal computer with an Intel ® Core TM i7-4700HQ @ 2.4 GHz processor with 8Gb RAM, compared with about 950 s for the other methods. When g=0.80, the CPU time for SHDOM was about 300 s, compared with 370420 s for the other methods. Finally, in the case of g=0.93, the CPU time for the SHDOM was about 360 s, exceeding a little the computational requirements of the FVM (~300s), but remaining lower than those of KL (~540 s), HG1 (~840 s) and HG2 (~1030 s). The CPU time for the normalization procedure has not included, but it is marginal (less than 10s for HG1, and negligible for KL and HG2).

The decrease of CPU time when g increases from 0.20 to 0.80 is likely to be due to the fact that, in the former case, coupling between the different directions is stronger. When the asymmetry factor increases, only the directions closer to the direction of propagation play a role, while the remaining directions have a negligible contribution to the in-scattering term.

However, when the asymmetry factor approaches unity, the contribution to the in-scattering term of directions close to the direction of propagation becomes increasingly larger, and since directions different from that of propagation are treated implicitly, the number of iterations required to achieve convergence is likely to increase. This may explain the increase of CPU time observed for most methods when g increases from 0.80 to 0.93.

Test case 2

In the second test case, an axisymmetric cylindrical enclosure with radius R=1 m and height H=2 m is considered. The walls are black and cold, except the cylindrical surface, which has an emissive power of unity. The medium is cold, and scatters according to the Henyey-Greenstein phase-function. The calculations were performed using a uniform mesh with 25×25 control volumes. The S 16 quadrature was used for the DOM and the KL, HG1 and HG2 normalization methods, while 9 polar and 4 azimuthal angles per octant were used in the FVM. In the case of the SHDOM, N  and N  were set to 18 and 16, respectively. Hence, the same number of discrete directions is used for all methods. We first consider a purely scattering medium with an optical thickness, based on the radius of the cylinder, equal to 25.

This test case was proposed by Hunter and Guo [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF], who compared the FVM and DOM using the same grid size and angular discretization considered here. The results obtained using the STEP scheme, which are not shown here, closely match those reported in [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF] when the phase function is normalized according to HG1 method. When the CLAM scheme is employed, the predicted net heat fluxes on the hot boundary are lower than those obtained using the STEP scheme, particularly for the lowest value of the asymmetry factor (0.8).

Figure 3 shows a comparison between the net heat fluxes on the hot wall, normalized by its emissive power, obtained for the mesh and quadrature mentioned above, and for a much finer mesh with 100×100 control volumes and angular discretization, with 12 polar angle and 12 azimuthal angle subdivisions per octant. The DOM was used along with the CLAM scheme and the HG1 method. The figure reveals a very low sensitivity of the results obtained to the spatial and angular grid refinement. Hence, the results displayed in figure 4, where the different methods are compared, were obtained using the coarsest discretization.

Figure 4 shows that all methods perform similarly for g=0.80, as far as the accuracy is concerned. When g increases to 0.90, the solutions obtained using HG1, HG2, FVM and SHDOM are in good agreement with each other. The KL normalization method overestimates the net heat fluxes by about 20% at z=H in comparison with HG1. In the case g=0.95, the differences between the HG1, HG2, FVM and SHDOM remain small. The KL predictions are about 20% higher than those of the HG1 method in the middle of the cylinder, in agreement with previously reported results [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF].

Additional calculations were performed for an asymmetry factor of 0.95, an albedo of 0.90, and three different values of the optical thickness of the medium. Our predictions obtained using HG1 normalization method are in close agreement with those reported in [START_REF] Hunter | Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media[END_REF]. The influence of the grid size and angular discretization is relatively small, as evidenced in figure 5, obtained using the same normalization procedure, but the presence of ray effects is visible when the optical thickness of the medium is equal to unity, despite the fine discretizations employed. Figure 6 shows the comparison of the various methods employed for the coarsest discretization. The results obtained using KL, HG1 and HG2 are very close to each other for =1.0 and =5.0. The FVM and SHDOM yield smooth results for =1.0, i.e., ray effects are not visible, and the predictions are in good agreement with those obtained using the DOM for both =1.0 and =5.0. When the optical thickness increases to =25.0, the KL method overpredicts the DOM solutions for HG1 and HG2 normalization methods by about 5%, while the FVM underpredicts those DOM solutions by 1.5%. The SHDOM closely matches the DOM predictions.

The CPU times required for the solution of the RTE when the DOM along with HG1 and HG2 normalization methods are used are generally close to each other. The FVM and SHDOM require more CPU time than the DOM (along with HG1 or HG2 methods) for =25.0, no matter the asymmetry factor, but lower CPU time for =1.0 and =5.0, while an opposite trend was found for the DOM along with KL method. Regardless of the method used, the CPU time required to solve the RTE is typically lower than 25 s for g=0.95, and lower than 55 s for other asymmetry factors, but depends on the optical properties of the medium. These times are much smaller than in test case 1, because the present problem is two-dimensional.

The CPU times required for the KL and HG2 normalization of the scattering phase function are negligible, but those of the HG1 normalization procedure (about 110 s) and of the evaluation of the average scattering phase function in the FVM (about 170 s), largely exceed the CPU time needed to solve the RTE for the present test case.

Test case 3

In the last test case, a three-dimensional cubic domain of side L=1 m, containing a scattering and absorbing homogeneous medium, is studied. The medium is subjected to an incident short pulse laser source normal to the boundary x=0, while the remaining boundaries are nonreflecting and cold (non-emitting). The medium is also cold, with an optical thickness of 10 and a scattering albedo of 0.5. The scattering is described again by the Henyey-Greenstein phase function with an asymmetry factor of 0.90. The collimated intensity of the short pulse laser, denoted by I c , is Gaussian in time and uniform in space, according to:

            p p p o c t t , t t t I t , , z , y , x I 6 0 3 2 ln 4 exp 0 2        i s s  ( 24
)
where  is the Dirac delta function and I o is the maximum radiative intensity of the pulse, which occurs at t = 3t p . The dimensionless time parameter, defined as

p p t c t   * , with c
standing for the speed of light, is set equal to 0.5. After t = 6t p the medium is free from irradiation. In this test case we are concerned with the transmitted, T, and reflected, R, temporal signature of the incident radiation, which are given by:

       dz dy d t , , z , y , L x I t , , z , y , L x I L I t T d c L L o                 s i s s s  2 0 0 2 1 (25)             L L d o dz dy d t , , z , y , x I L I t R 0 0 2 2 0 1 s i s s  (26)
where I d is the diffuse component of the radiation intensity.

The temporal discretization was performed using a second-order Runge-Kutta scheme. The domain was discretized using a grid with 200×20×20 control volumes and an S 12 quadrature, and the time step t was defined as t=.x/c, where the stability parameter  was set equal to 0.5. A previously validated Monte Carlo algorithm [START_REF] Roger | A hybrid transport-diffusion model for radiative transfer in absorbing and scattering media[END_REF] was used to provide benchmark results.

The KL and HG1 normalization methods may be readily extended to account for the collimated component [START_REF] Hunter | Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer  Part I: Conservation of scattered energy and asymmetry factor[END_REF], as well as the FVM and SHDOM. However, the HG2 normalization procedure is not directly applicable, since the collimated component does not belong to the set of discrete directions, and therefore it is not feasible to normalize the forward and backward components of the phase function. However, the original idea underlying this method may be maintained by applying the normalization to the discrete directions with the scattering angle closest to 0º and 180º, as described in [START_REF] Hunter | Applicability of phase-function normalization techniques for radiation transfer computation[END_REF]. A drawback of this implementation is the possible occurrence of unphysical negative values of the normalized scattering phase function for a scattering angle close to 180º, as reported in [START_REF] Hunter | Applicability of phase-function normalization techniques for radiation transfer computation[END_REF] and confirmed in the present test case. Despite this, it was possible to obtain a converged solution, due to the minor role played by the backward-scattering terms in the case of a large and positive asymmetry factor.

The results obtained are displayed in figure 7. If the DOM is used without any normalization method for the phase function, the transmittance is overestimated by more than one order of magnitude, while the decay of the reflectance after the peak is also strongly overestimated.

Even though all other methods yield much better solutions, the discrepancies among them, and in comparison with the Monte Carlo reference solution, are larger than for the previous test cases. The KL normalization procedure provides a surprisingly good estimation of the transmittance, but does not perform well as far as the reflectance calculation is concerned.

The HG1 and HG2 methods underpredict the peak of the transmittance by about 15%, but match well the reference solution for shorter and longer times. The HG1 method further yields an accurate prediction of the reflectance, while the HG2 procedure does not, due to the problem with the normalization of the collimated radiation component mentioned above. The FVM predictions of the transmittance are similar to those of the DOM with HG1 or HG2 methods. The former method marginally overpredicts the peak of reflectance, in contrast with the latter, which slightly underestimates that peak, in comparison with the MC reference solution. The SHDOM provides the best estimate of the transmittance, and predicts well the reflectance.

Both the DOM with any of the tested normalization methods and the FVM require about 26 h of CPU time for this transient problem, while the SHDOM is far more efficient, and needs only about 11.5 h.

Approximation of the scattering phase function

We will now examine whether the approximation of the Henyey-Greenstein phase function by the phase functions mentioned in section 2.6 is satisfactory or not for the previous test cases.

Figures 8(a) and 8(b) show the predictions obtained using the DOM for test cases 1 and 2, respectively, and corresponding to the set of radiative properties considered in Figs. 2 and4, respectively. In the case of the Henyey-Greenstein phase function, the HG1 normalization method was used, while no normalization is required for the transport, delta-Eddington and delta-M approximations.

The transport and delta-Eddington approximations perform quite well for asymmetry factors up to 0.9 in the two test cases. Therefore, they can be used instead of the normalization procedures discussed above. The delta-Eddington approximation performs better than the transport approximation, as expected, since * is linearly anisotropic in the former approximation and isotropic in the latter. The delta-Eddington and delta-M phase functions yield similar results, which are graphically indistinguishable, and therefore only the results calculated using the former phase function are presented.

When the asymmetry factor is greater than 0.90, the differences between the results for the different phase functions become larger. Although the transport and delta-Eddington approximations remain an attractive option for test case 1, improved results are obtained using the delta-M with M=2, denoted as delta-M2 in the figure. The results for delta-M approximation are not shown for test case 2, since they are almost coincident with those of the delta-Eddington approximation.

In test case 3, however, the results obtained using the transport and delta-Eddington approximations are rather different from those calculated using the Henyey-Greenstein phase function, as illustrated in Figs. 8(c) and 8(d). The transmittance peaks are overestimated by a factor of two or more. Moreover, the maximum transmittance occurs earlier, and the subsequent decay is faster when those approximate phase functions are used. The reflectance peak calculated using the transport approximation is marginally overestimated, and that computed using the delta-Eddington phase function is significantly underestimated. Accordingly, these approximate phase functions should not be used in the case of transient radiative transfer with collimated irradiation. However, the delta-M function provides much better results, particularly for M=6, even though the transmittance peak occurs slightly earlier in comparison with the solution computed using the Henyey-Greenstein phase function.

Overall, and since the Henyey-Greenstein phase function is also an approximation of an actual phase function, the delta-M may be used instead of the Henyey-Greenstein phase function, allowing the DOM to be applied to strongly forward scattering media without any normalization technique.

However, the value of M in the delta-M approximation needs to be carefully chosen. If it is too small, the approximation may be poor, while if it is too large, the weight of the Dirac delta function in Eq. ( 19) is small, and the results may be worst. In fact, even though the delta-M function approximates better the Henyey-Greenstein phase function when M increases, the asymmetry factor of * increases as well. When this asymmetry factor is large enough and the DOM is used, the error in the calculation of the scattered energy and asymmetry factor tends to increase. In such a case, the phase function needs to be normalized to avoid poor predictions, but then there would be no advantage in replacing the Henyey-Greenstein by the delta-M phase function, as pointed out in section 2.6. In test case 3, the results obtained using M=8 (not show in the figure) are worst then those computed using M=6, because the asymmetry factor of * for M=8 is large (g=0.87).

CONCLUSIONS

An investigation of different strategies to deal with highly anisotropic scattering media in the framework of the DOM was carried out. Different proposals for the normalization of the phase function previously reported in the literature were compared with two other approaches that do not rely on the phase function normalization. One of them is based on the finite volume method for the evaluation of the in-scattering term of the RTE, and the other is based on the spherical harmonics discrete ordinates method, formerly used in the field of atmospheric radiation. The use of approximate phase functions given by the sum of a Dirac delta function and a smoother phase function is also investigated as an alternative to the normalization of the strongly forward scattering phase function.

The KL normalization procedure does not perform well for highly anisotropic scattering, in agreement with past work, since it does not preserve the asymmetry factor of the phase function. The HG1 method was highly effective for all studied cases, but the computational requirements of the normalization procedure may exceed those needed for the solution of the RTE. The HG2 technique overcomes that drawback, and provides solutions of similar accuracy for diffuse radiation. In the case of collimated radiation, however, this technique may yield negative normalization coefficients and physically unrealistic values of the normalized scattering phase function. These may lead to negative values of the radiation intensity. Even though this did not prevent convergence for the test case selected in this work, the predictions of the reflectance for the transient test case were not satisfactory. Therefore, this normalization procedure should be used with caution for collimated radiation.

The FVM is an interesting alternative to the normalization procedures, since it does not modify the phase function. However, an accurate integration is needed to obtain the average phase function for every pair of discrete directions, before and after scattering, and this may be time-consuming in the case of acute forward scattering. In the present work, an adaptive integration method was employed, which yields satisfactory results for an asymmetry factor up to 0.95.

The SHDOM performed quite well for all test cases, and the computational cost was the lowest, except in the two-dimensional axi-symmetric problem and for optically thick media. It proved to be a very interesting alternative to the DOM anisotropically scattering media, and deserves further investigation.

The approximation of the Henyey-Greenstein phase function using the transport or the delta-Eddington approximation is very satisfactory when the transient term of the RTE is negligible and the asymmetry factor does not exceed 0.90, and still performs reasonably for higher asymmetry factors, avoiding the normalization of the phase function. In the case of transient problems with collimated irradiation, those approximations do not perform well, but the more accurate delta-M phase function yielded very satisfactory results. Therefore, it constitutes an interesting alternative that relaxes the need to normalize the phase function when the DOM is applied to strongly forward scattering media. However, the value of M needs to be carefully chosen, and should be as large as possible, provided that the scattered energy and the asymmetry factor are accurately evaluated.

FIGURE CAPTIONS

  polynomials, and u(m) = cos(m) if m≥0 and u(m) = sin(m) if m<0. The in-scattering term is transformed to discrete ordinates as follows:

Figure 1 .

 1 Figure 1. Influence of the grid size (a) and quadrature (b) on the incident heat flux, normalized by the

Figure 2 .

 2 Figure 2. Incident heat flux, normalized by the emissive power of the hot boundary, along a central

Figure 3 .

 3 Figure 3. Influence of the spatial and angular discretizations on the net heat flux on the hot

Figure 4 .

 4 Figure 4. Net heat flux on the hot boundary, normalized by its emissive power, along the axial

Figure 5 .

 5 Figure 5. Influence of the spatial and angular discretizations on the net heat flux on the hot

Figure 6 .

 6 Figure 6. Net heat flux on the hot boundary, normalized by its emissive power, along the

Figure 7 .

 7 Figure 7. Time history of the dimensionless transmittance and reflectance for test case 3.

Figure 8 .

 8 Figure 8. Influence of approxiamted scattering phase functions.

Figure 1 .

 1 Figure 1. Influence of the grid size (a) and quadrature (b) on the incident heat flux,

Figure 2 .

 2 Figure 2. Incident heat flux, normalized by the emissive power of the hot boundary, along a

Figure 3 .

 3 Figure 3. Influence of the spatial and angular discretizations on the net heat flux on the hot

Figure 4 .

 4 Figure 4. Net heat flux on the hot boundary, normalized by its emissive power, along the axial

Figure 5 .

 5 Figure 5. Influence of the spatial and angular discretizations on the net heat flux on the hot

Figure 6 .

 6 Figure 6. Net heat flux on the hot boundary, normalized by its emissive power, along the axial

Figure 7 .

 7 Figure 7. Time history of the dimensionless transmittance and reflectance for test case 3.

Figure 8 .

 8 Figure 8. Influence of approxiamted scattering phase functions.

ACKNOWLEDGEMENTS

This work was financially supported by national funds by FCT-Fundação para a Ciência e Tecnologia under the project PTDC/EMS-ENE/1028/2012, and through IDMEC, under LAETA, project UID/EMS/50022/2013.


Amelia Carolina Sparavigna

Mauro Maria Baldi

Symmetry and Golden Ratio in the Analysis of Regular Pentagon

Keywords: Geometry, Symmetry, Mathematics, Golden Ratio, Mathematics Education

Regular pentagon is a geometric figure which had a relevant symbolic meaning in Pythagorean and Platonic philosophies and a subsequent important role in the Western thought, appearing also in arts and architecture. A property of the pentagon, which was probably found by Pythagoreans, is that the ratio between diagonal and side is the golden ratio. Here, we will study some links between regular pentagon and golden ratio. We will focus first on the group of five-fold rotation symmetry, to find the position of the vertices of this geometric figure in the complex plane; then, we will give an analytic method to solve the same problem in the Cartesian plane. In the analytic approach, we find the golden ratio without any specific geometric consideration. This study can be interesting also for a general education in mathematics, because it can convey and link several concepts, requiring only a general pre-college education.

Introduction.

In geometry, pentagons are five-sided polygons, which can be regular when the interior angles are equal to 5 / 2π π -, that is, 108°. With regular pentagons we can construct one of the Platonic solids, the dodecahedron. Plato wrote about the solids having congruent regular polygons as faces in his Timaeus, c.360 B.C., associating them to the four classical elements (earth, air, water, and fire) [START_REF] Kotrč | The Dodecahedron in Plato's "Timaeus[END_REF]. The dodecahedron is the fifth solid, of which Plato told simply: "There was yet a fifth combination which God used in the delineation of the universe" [START_REF]Plato's Timaeus[END_REF]. Besides this role in Plato's cosmology, the regular pentagon had a general importance in geometry, art and philosophy, probably because of the fact that its diagonals are in the golden ratio to its sides. We find the pentagon, for instance, in Agrippa's De Philosophia Occulta, to illustrate the proportion, measure and harmony of human body [START_REF] Cornelius Agrippa Von Nettesheim | 1486 -1535) was a German Polymath, physician, legal scholar, soldier, theologian and an occult writer. He used a regular pentagram for the proportion of the man's body[END_REF].

The property that the ratio between diagonal and side is the golden ratio ϕ, that is 2 / ) 5 1 ( + = φ , was first reported in a written text at Euclid's time. Euclid gives two constructions of regular pentagon in Book IV, as Propositions 11 & 12. According to Sir Thomas Little Heath, who translated the Euclid's works, these methods had been probably previously developed by the Pythagoreans [START_REF] Heath | The thirteen books of Euclid's Elements[END_REF]. Besides this link to the golden ratio, another important feature related to the Pythagorean philosophy is observed in a regular pentagon: diagonals are creating another geometrical figure, the pentagram, which is the symbol of Pythagoreanism [START_REF] Ghyka | The geometry of art and life[END_REF]. As told in [START_REF] Fossa | On the pentagram as a Pythagorean emblem[END_REF], a striking visual feature of the regular pentagon was that its diagonals create the pentagram and, in it, a new pentagon inside the original pentagon. This self-reproductive property of the pentagon was probably evocative to Pythagoreans of the transmigration of souls [START_REF] Fossa | On the pentagram as a Pythagorean emblem[END_REF].

With its five-fold symmetry, the pentagon attracted an artist, Albrecht Dürer, and the astronomer Johannes Kepler, for the problem of polygonal tiling of the plane using pentagons [START_REF] Caspar | Five-fold symmetry in crystalline quasicrystal lattices[END_REF]. In his Harmonices Mundi, when discussing on the congruence of harmonic figures, Kepler investigated the symmetry of polyhedra and explored orderly arrangements of plane pentagons. According to [START_REF] Caspar | Five-fold symmetry in crystalline quasicrystal lattices[END_REF], this was an anticipation of Penrose aperiodic tiling models, that served for the geometry of quasicrystal structures.

Since regular pentagons and pentagrams have a large impact on arts and philosophy, we should not be surprised if the pentagon is also a symbol of the Sacred Geometry, a discipline ascribing sacred meanings to geometric shapes and proportions [START_REF] Skinner | Sacred Geometry: Deciphering the Code[END_REF][START_REF] Lawlor | Sacred Geometry: Philosophy and Practice[END_REF]. In a previous paper, we have already discussed one of these symbols, the Vesica Piscis, a shape coming from the intersection of two circles having the same radius, arranged so that the center of each circle lies on the perimeter of the other [START_REF] Sparavigna | A Mathematical Study of a Symbol: the Vesica Piscis of Sacred Geometry[END_REF].

Here, we will discuss some properties of the regular pentagon; as we did in [START_REF] Sparavigna | A Mathematical Study of a Symbol: the Vesica Piscis of Sacred Geometry[END_REF], we will focus our analysis on the algebra that we can obtain from the symbol. We will see that the study of such symbol can be interesting also for a general education in mathematics, because it can convey several concepts of mathematics to people having just a pre-college education. In particular, we will link the symmetry of the pentagon to an analytic method to solve it, demonstrating some properties linking the pentagon to the golden ratio. Let us start our discussion from the symmetry of pentagon (with its five-fold rotation group), a symmetry which is also observed in nature (Fig. 1). The five-fold symmetry. The regular pentagon possesses a five-fold symmetry, that is, a rotational symmetry of order 5, given by rotations in the plane of the polygon about its center of angles 2π/5, 4π/5, 6π/5 and 8π/5. The rotation of 10π/5 = 2π is the identity. It has also five lines of reflection symmetry.

A symmetry is coming from a transformation of a planar figure, which is a rigid motion that leaves the figure unchanged. New transformations can be obtained by following one symmetry after another. In this way, an algebra of symmetries is obtained, where the multiplication is given by the composition of symmetries. This algebra contains the identity element and the inverse elements. In other words, we have a symmetry group for each geometric figure, because every figure has at least the symmetry group consisting only of the identity. Let us consider an example; a rectangle, which is not a square, has four elements in its symmetry group, namely rotations about an axis perpendicular to the figure and passing through its center, of angles 0 and π, and a horizontal and vertical flip about the midlines. This is the Klein 4-group. A square, which is more symmetric, has 8 elements in its symmetry group.

Let us remember that a regular polygon has all its sides having the same length and all its angles the same value. Then, we can have an equilateral triangle, which is a regular 3-gon, a square which is a regular 4-gon, a pentagon a regular 5-gon and so on. The symmetry group of a regular n-gon has n rotational and n reflection symmetries, for a total of 2n symmetries. It is called the dihedral group of the n -gon and denoted by D n . The pentagon has a dihedral group of 10 elements. There are also two subgroups, consisting of the five rotation symmetries and the subgroup of the five mirror symmetries.

Here, we consider a representation of the rotation group of the pentagon in a complex plane. This plane is a geometric representation of the complex numbers established by the real axis and the orthogonal imaginary axis. The position of the vertices of a pentagon can be described by complex numbers, which can be assumed in a rectangular form, ) ( ib a + , or in a polar form or in the exponential form. This last form is giving the position of a point Z in the complex plane by the Euler's identity:

) sin (cos θ θ θ i A Ae Z i + = = .
Using the pentagon in Figure 2(a), inscribed in a circle of unit radius, positions of vertices O, A, B, C and D are given by the following numbers: . These numbers are also the representation of the elements of rotation group, that, with the multiplication, give us an algebra.

Here two examples of this algebra: let us consider two elements of the group, for instance Let us consider diagonal AD and side BC. Using Euler's identity, the ratio between their lengths in the complex plane has the value:

= =       - = = = BC AD (1)
It means that, following Euclid who told that the ratio between diagonal and side is the golden ratio, we have that φ is two times the cosine of 5 / π . However, from Figure 2(a), we have also immediately that the apothem of the regular pentagon, inscribed in the goniometric circumference, is equal to the cosine of

+ + - = + + - = + + + - = + + + - = + = i i i i i e i π π π (2)
We can repeat the same calculation for the other four vertices.

However, pentagon can appear as in the Figure 2(b), rotated by 90 degrees. We have that the position of A is now given by: . So, we have that triangle ABC is isosceles. Triangles OAB and ABC have congruent angles, therefore they are similar. We have: OA / AB = AB / BC. Since AB = AC = OC = x and OA = OB = 1, we have:

) 1 ( 1 x x x - = = φ (4)
Eq.4 is the equation that defines the golden ratio ϕ. This means that AB = OC is the golden ratio of the radius. From (4), we have

x 2 = (1 ─ x) and 2 / ) 5 1 ( + - = x
. Then: . From this result, (2) had been evaluated.

Using the rotation group and complex numbers, and also a little bit of geometry concerning the golden ratio, we have seen how to find the vertices of regular pentagon. From them, all the other properties can be obtained. ( )

2 1 2 2 2 2 1 2 2 2 2 1 1 1 1               - - + =               - + - + - β β β α β α
If we square both sides of this equation, we have to guarantee that one of the two radicands must be non-negative, in order to avoid the introduction infeasibility due to the squaring operations. In other words, we have to solve the following system:

( ) ( )        ≥       - + - + -       - - + =       - + - + - ) ( 0 1 1 ) ( 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 b a β α β α β β β α β α
(b) is always satisfied, because it is the sum of squared quantities. Thus we have to consider (a), that is:

( ) 2 2 2 2 2 2 2 1 1 1 1       - - + =       - + - + - β β β α β α (5)
Let us also consider that the regular pentagon is made of five congruent triangles. Triangles are congruent when they have exactly the same sides and angles. In this manner, a triangle can become another triangle using turns and flips. In the case of Figure 4, we can start from triangle AOE and, after a rotation, having triangle AOB. Then, in the same manner, we can have BOC. In fact, as we have discussed before, the regular pentagon has an order of rotation symmetry equal to 5. Since AOE is isomeric to BOC, the areas of the two triangles are equal. Thus:

( ) 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 α α β β α α - ⋅ = ⇒ ⋅ ⋅ = - ⋅ ⋅ ⇒ ⋅ = ⋅ MB OC HO AE (6)
We have two equations ( 5) and ( 6) and two unknown variables β α , .

( )

       - ⋅ =       - - + =       - + - + - 2 2 2 2 2 2 2 2 1 2 1 1 1 1 α α β β β β α β α (7)
After simple passages:

( )( )      - ⋅ = - - = - - 2 2 2 2 1 2 1 1 1 α α β β αβ β α (8)
A shortcut in the resolution of this equation is the squaring of both sides. The squaring operation gives no problems; however, it must be performed with care. Since 0 > α , we have that:

α β α 2 1 2 = -
If we square both sides of this equation, we also have to guarantee that β/2α ≥ 0; this condition is already satisfied because, by construction, we have that α, β > 0.

Therefore, from the second equation in [START_REF] Skinner | Sacred Geometry: Deciphering the Code[END_REF], we have:

( ) ( ) 2 2 4 2 2 4 2 2 2 2 2 1 4 4 1 1 4 4 1 4 α α α β α α α α β - = + - = - ⇒ - = - = (9)
Using [START_REF] Lawlor | Sacred Geometry: Philosophy and Practice[END_REF] in the first equation of ( 8), with the substitution , we have a noteworthy simplification. However, radicals in (8) must be handled with caution. In fact we have that:

     - ⋅ = - - - ⋅ = - - 2 2 2 2 2 2 1 2 2 1 1 2 1 2 1 α α β α α α α α (10)
It is therefore necessary to distinguish two possible cases:

0 2 1 2 < -α or 0 2 1 2 > -α
. Let us consider the second case which is that fitting our assumptions (the first case is discussed in the Appendix). ( 10) becomes:

( ) ( )      - ⋅ = - - - ⋅ = - - 2 2 2 2 2 2 1 2 2 1 1 2 1 2 1 α α β α α α α α (11)
In the first equation of (11), we can use the substitution 2 α = t

. The system becomes:

:

( )        ≥ - - = - = - - (**) 0 2 1 ) 1 ( 2 (*) 2 1 1 1 4 t t t t t t β
Squaring (*) requires the following condition:

) 0 2 1 0 1 4 ( ) 0 2 1 0 1 4 ( < - ∧ < - ∨ ≥ - ∧ ≥ - t t t t
. The second part of the condition cannot be satisfied because it violates assumption (**). Therefore, squaring both sides of (*) implies the resolution of the following system:

( ) ( )        ≤ ≤ - = - = - - 2 1 4 1 ) 1 ( 2 2 1 ) 1 ( 1 4 2 2 t t t t t t β (12)
After some passages we have:

( ) 0 5 20 16 2 = + -t t t
(13) Equation ( 13) has three solutions. One solution (t = 0) implies α = 0, which is in contrast with assumption 0 < α < 1. The other two solutions are: After some simple passages, we have the desired coordinates of the five vertices of pentagon: Note the remarkable power of the analytic method; we have obtained the positions in the Cartesian reference, without needing sine and cosine of angles or the golden ratio.

( )         + - - - =       - - -         - + - =       - - =         - + =       -         + - - =       - -
The apothem and the golden ratio. In the previous discussion, the golden ratio does not explicitly appear. However, it must be obtained also from the analytic calculation.

And in fact, let us evaluate the ratio between diagonal d and side l of the pentagon (see Figure 4 (16) shows that the ratio between diagonal and side is the golden ratio.

Moreover, in the analysis of the previous section, we have found a relevant square root, which is also the length of OH, the apothem of the regular pentagon (see Figure 4):

+ = + = - = = α a OH (17)
It is easy to verify that the following identity exists:

2 2 2 2 B A A B A A B A - - ± - + = ± (18)
Squaring the members of (18), we have:

4 2 2 2 2 2 2 2 2 2       - -       - + ± - - + - + = ± B A A B A A B A A B A A B A
If we use the identity (18): . Therefore, we have:

Conclusion.

In this paper we have discussed some properties of the regular pentagon, focusing first on its group of rotation symmetry, to find the position of the vertices in the complex plane, then on an analytic method to solve the same problem in the Cartesian plane. In the solution of the analytic method, we do not need the use of trigonometric functions or geometry for the golden ratio φ , if this ratio is assumed as solution of equation . We have also shown, from the analytic method, that the ratio between the apothem of a regular pentagon and the radius of its circumscribed circumference is one-half the golden ratio.

Appendix. For Eq. [START_REF] Sparavigna | A Mathematical Study of a Symbol: the Vesica Piscis of Sacred Geometry[END_REF], it is necessary to distinguish two possible cases:

Figure 1 :

 1 Figure 1: The five-fold symmetry (pentamerism) is common in nature. Many flowers possess it as shown by the picture on the left. Note the pentagon in the central part of the flower. On the right, pentamerism in an apple (Courtesy: Rasbak, Wikipedia).

Figure 2 :

 2 Figure 2: Positions of the vertices in the complex plane for a pentagon (a), and the same figure rotated of 90 degrees (b).

  5 /πand then it is one-half the golden ratio.Let us conclude considering the rectangular representation of the point A with 5 2π i e . It is given by:

  so on, for all the vertices.Golden ratio. To have equation (2), we need to use a property of the remarkable angle γ = π/5 = 2π/10 = 36°. Let us define η=2γ=72°, and consider Figure3.

Figure 3 :

 3 Figure 3: AB = OC = x is the golden ratio of radius OB.

Figure 4 :

 4 Figure 4: Coordinates of vertices in the Cartesian plane, for a pentagon inscribed in a circumference of radius 1. Note that β α 2 , 2 = = BD AE .

  that the golden ratio φ is the positive root of equation

.

  second case was considered in the paper. Here we discuss the first case;[START_REF] Sparavigna | A Mathematical Study of a Symbol: the Vesica Piscis of Sacred Geometry[END_REF] becomes: Zero and the negative solution cannot be acceptable, because we are requiring 0 > α . If we use the positive solution, we find after some calculations that2 / 3 ± = β. The negative root is not acceptable. But also the positive, because we should have 2