© M. Weiss et al., Published by EDP Sciences, 2024

1 Introduction

Space agencies and private companies have several projects for lunar bases [1,2]. Structural developments of facilities, consideration of the Moon environment, work breakdown structures have already been studied to prepare future missions [1,2]. The space-like atmosphere of the Moon [3] and the necessary self-sufficiency due to the difficulty to supply materials from Earth led to consider photovoltaics (PV) as an energy source to electrically supply such lunar facilities. To determine the implementation sites and design the photovoltaic electricity production systems, it is necessary to predict the photovoltaic potential at any place on the Moon, for different technologies or installation modes. For that, researchers have been interested in solar illumination conditions [4–9], which represent the fraction of time when the Sun is visible at a defined place and for a given period. In the literature, Digital Elevation Models (DEMs) of the lunar topography are used to represent altitude variation of the surface and to anticipate shadows from the horizon topography [4–9]. DEMs are based on data collected from different sources such as Earth-based radars [4], lunar orbital satellite [5,7,9] or a combination of both [6]. Once the topography was determined, simulations of the Sun motion in the lunar sky were used to create illumination maps for a given period. Bryant assumed that, during an 18.6-year period, corresponding to the majority of the possible Sun-Moon relative positions during the period of the regression of the longitude of the ascending node of the Earth, the center of the Sun passes above all the longitudes and latitudes inside a rectangle of 180° in the East-West direction and only 1.54° in the North-South direction from the equator [4]. Thus, he did not compute the exact position of the Sun but used the area of the rectangle to compute illumination. His assumptions were verified thanks to the Jet Propulsion Laboratory HORIZONS tool, available online and able noticeably to compute the elevation and the azimuth of the Sun in the lunar sky from any location on the Moon [10]. Mazarico et al. used the SPICE toolkit to determine the position and the angular radius of the Sun at each step time in their computing algorithm [9,11]. The different methods listed above have determined points of interest where the illumination is predicted to be higher than 80% of the time for at least an 18.6-year cycle [4,7,10]. The points of interest can be different from one article to another but they are all located either near the North or the South pole of the Moon.

Solar irradiance and solar energy have been computed for the global lunar surface and for various specific locations on the Moon. Li et al. give the formula to compute the solar irradiance received by the Moon depending on the Moon-Sun distance [8]. Furthermore, by a simulation of an 18.6 yr cycle, they conclude that for a horizontal surface, the maximum monthly effective irradiance, computed using the angle between the surface and the Sun's rays, is higher near the lunar equator and decreased towards the poles. They also highlight that the lunar effective irradiance evolves in opposite direction between the northern and the southern hemisphere, due to the axial tilt of the Moon. Kaczmarzyk et al. proposed a model and estimated the solar energy received in one lunar day by different oriented plane surfaces at different latitudes below 80° [12]. They estimate that the maximum theorical energy a surface can receive tracking the Sun is 1740 MJ/m² per lunar day, which represents around 483 kWh/m². This value would be 57% higher than the energy received by a fixed surface with optimum angles, which would therefore be around 1100 MJ/m² or 305 kWh/m².

Other works led by Kaczmarzyk highlights that the irradiance received by a surface pointing towards the Sun using a two-axis or a one-axis tracker with an optimal tilt angle are very similar, no matter the latitude, because of a very small Moon axial tilt (about 1.54°) [13]. Furthermore, they assume that the Moon has no season, due to the relatively low tilt angle of the Moon. It results that a fixed solar array has almost the same production all along the terrestrial year.

Therefore, even if the photovoltaic potential on the moon has already been considered, there is no systematic method to compute the photovoltaic potential at any location on the Moon and to compare different modes of installation.

In this work, we predict numerically the solar energy received for any period of time by a surface located anywhere on the Moon. While power is an instantaneous data, we want to focus on a long-term study, like for several terrestrial years, this is why we choose energy instead of power. The computation is based on the position of the Sun in the lunar sky at any place and any time. We compute solar energy for 4 different types of Sun tracking modes: fixed, tracked on a horizontal or vertical axis and tracked on two axes. Highly illuminated site locations [4,5] are used as a case study. We have chosen high altitude locations, above 1000 m, to disregard lunar topography and minimize potential shadow effects.

2 Methods

The block diagram presented in Figure 1 describes the main steps of the algorithm we developed to obtain the irradiance on a flat surface at any instant t. Each block will be described in detail next.

First, the Earth-Moon-Sun system calculates for a given time t the selenographic longitude and latitude of the subsolar point defined as the point at the surface of the Moon where the Sun's rays are perpendicular to a horizontal surface. Longitude and latitude are counted positively towards East and North, respectively. The range of longitudes is [−180°, 180°] and [−90°, 90°] for latitudes. The Earth-Moon-Sun system also provides the lunar global irradiance. Then, the subsolar point coordinates are the inputs of the lunar sky block. Inputting selenographic coordinates of a point into this block allows to get the elevation and the azimuth of the Sun in the lunar sky for this point and for any instant.

The next step of the algorithm is to find the incidence angle between the Sun's rays and the flat surface at any time. The orientation of this surface is a controllable input composed of the azimuth and the tilt angles. The orientation can evolve at any time t.

Finally, the final irradiance is obtained by combining the lunar global irradiance, computed in the Sun-Moon system block, with the incidence angle previously determined. The elevation angle of the Sun is also used to introduce rules about the visibility or not of the Sun from the surface.

Fig. 1 Block diagram of the main algorithm. Green inputs are controllable and the red output is the final output of the algorithm.

2.1 The Earth-Moon-Sun system

2.1.1 The Earth-Sun and Earth-Moon systems sub-blocks

The Earth-Moon-Sun system is divided into four sub-blocks. The first two sub-blocks are called Earth-Sun and Earth-Moon systems. They compute in parallel the position of the Sun and the Moon, respectively, in the geocentric coordinates as well as their center-center distance to the Earth, for any instant t given as an input. The calculation uses the astronomical algorithms of Meeus [15,16]. The outputs of the Earth-Sun sub-block are the dimensionless Earth-to-Sun distance R_ES and the apparent geocentric longitude of the Sun λ_ES. The outputs of the Earth-Moon sub-block are the dimensionless Earth-to-Moon distance R_EM and the geocentric ecliptical longitude and latitude of the Moon λ_EM and ϕ_EM, respectively (cf. Fig. 2). Both dimensionless distances are in astronomical units (AU) and the coordinates are angles in degrees.

Furthermore, the Earth-Moon sub-block computes the physical librations of the Moon, thanks to the theory of Eckhardt and the Meeus astronomical algorithms [15–17]. There are three librations and they will be used later in the Moon-Sun system.

Fig. 2 Angles and distances between bodies in the Earth-Moon-Sun system. The primary direction is the direction where the longitude is 0°, corresponding to vernal equinox [14].

2.1.2 The Sun-Moon and Moon-Sun systems sub-blocks

The next two sub-blocks of the earth-moon-sun system are called Sun-Moon and Moon-Sun. The Sun-Moon system block has two functions in the algorithm. First, it computes the lunar solar irradiance. According to Li et al., the formula for the lunar solar irradiance (LSI, in W/m²) is [8]:

(1)

where S₀ is the solar constant (fixed to 1361.8 W/m²) [8] and R_SM is the dimensionless Sun-to-Moon distance relative to one astronomical unit (AU). The R_SM can be deduced from R_EM and ϕ_EM using:

(2)

The LSI can be computed for any period. As an example, Figure 3 presents the LSI at the moon equator for year 2022. We can see that the LSI fluctuates along the year with a small amplitude (+/− 0.3% in average) and a short period (∼14 days), and a large amplitude (+/− 3.5%) and a long period (∼1 terrestrial year). The longest period of oscillations corresponds to one lunar year while the shortest period corresponds to approximately half of a lunar sideral month, the time the Moon takes to complete one orbit around Earth.

Since the atmosphere of the moon is extremely thin, the light absorption and the scattering by particles in the lunar sky can be neglected [3]. In addition, the moon radius can be neglected compared to the Moon-Sun distance. Therefore, two surfaces facing the Sun at an instant t are considered to receive the same amount of solar power density, no matter their respective position on the Moon, as long as the Sun is visible. In the rest of the algorithm, we consider that the irradiance received by any illuminated point on the surface of the Moon is equal to LSI.

The second function of the Sun-Moon system block is to compute the heliocentric ecliptical longitude and latitude of the Moon λ_SM and ϕ_SM, respectively. The relations obtained from the Sun-Moon system block are the inputs of the computation and the formula is the following [18,19]:

(3)

(4)

The heliocentric coordinates of the Moon combined with the Moon's physical librations obtained in the Earth-Moon system block give, through the Moon-Sun system block, the optical and physical longitudes and latitudes of the Sun in the selenographic coordinate system. Then, summing the optical and the physical coordinates gives the final selenographic position of the subsolar point [18,19]. The outputs of the Moon-Sun system are then the subsolar point selenographic longitude and latitude λ_Msp and ϕ_Msp,respectively.

Fig. 3 Lunar global irradiance received in the lunar year 2022 (8520 h) by the Moon at the equator.

2.2 The lunar sky block

The Sun's position in the sky can be defined by the azimuth and the zenith angles (Fig. 4). In our computation, the North represents the azimuth zero and other azimuth angles are counted positively in the clockwise convention.

To get the azimuth of the Sun, the algorithm uses the components of a vector pointing the Sun via the formula of Zhang et al. [21]:

(5)

(6)

where S_x and S_y are the vector coordinated in the observer's referential and (λ₀, ϕ₀) are the observer's coordinates. The Sun azimuth angle γ_Sun is then calculated in the North clockwise convention as the arctan2 function of S_x and S_y:

(7)

To get the elevation angle of the Sun, corresponding to 90° minus the zenith angle, we first compute intermediate terms a₁ and a₂ from Li et al. [22]:

(8)

(9)

where R_Moon is the mean radius of the Moon considered equal to 1737.1 km [23].

The elevation angle β_Sun is finally obtained using:

(10)

Fig. 4 Angles used to compute the incidence angle on a flat surface (inspired by [20]).

2.3 The incidence angle block

The objective of this block is to determine the solar incidence angle (θ_S) on a flat surface. θ_S is defined as the angle between the Sun's rays and the normal to the surface (Fig. 4). This angle is calculated using γ_Sun and β_Sun, as well as the orientation of the surface, represented by the surface azimuth angle γ_Surface and the surface tilt angle β_Surface.

By combining angles of the Sun and angles of the surface, the incidence angle is then given by [20]:

(11)

2.4 The irradiance block

The last part of the main algorithm computes the irradiance received by a flat surface for an instant t, a location on the Moon and a surface orientation, the three controllable inputs of the algorithm. On Earth, the computation of the irradiance is a sum of three components, the direct, the diffuse and the reflected irradiances [24]. On the Moon, as mentioned before, there is negligible diffuse irradiance caused by particles or clouds occulting the Sun. Furthermore, due to a global surface albedo of 0.136 [25], the reflected irradiance can be neglected compared to Direct Normal Irradiance equal to the LSI previously computed in Section 1.2. Therefore, the global normal irradiance (GNI) is considered equal to LSI, providing the Sun is fully visible. To determine whether the Sun is fully visible or not, the altitude of the surface and the elevation of the Sun must also be considered.

First, the altitude of the surface adds a supplementary angle α_limit beyond horizon where the Sun remains visible. This angle α_limit can be computed thanks to R_Moon and the following formula [26]:

(12)

where h_point is the altitude of the point.

Furthermore, in the algorithm, the Sun is considered from the surface of the Moon as a sphere with an apparent diameter. We assumed that it is independent of the Moon-Sun distance and has a constant value of 0.53° (computed by deriving the formula from Weisstein [27]). Note that the moon topography is not considered here and that the Moon is assumed to be a perfect sphere.

Once the horizon limit angle and the apparent diameter are introduced in the algorithm, we can define the factor τ ∈ [0 ; 1] that indicates if the Sun is visible or not in the sky and following the rules:

• If β_Sun is at the limit angle  + the semi-apparent diameter and higher, it is fully visible. Therefore τ =1.

• If β_Sun is at the limit angle − the semi-apparent diameter or beyond, it is not visible. Therefore τ = 0.

• If the Sun is between the limit angle − the semi-apparent diameter and the limit angle + the semi-apparent diameter, the Sun is partially visible and therefore τ has a value between 0 and 1.

Therefore, the formula for τ is:

where δ_S= 0.265° is the semi apparent diameter of the Sun.

Then, we can compute the GNI thanks to τ using:

(13)

The evolution of the Sun's elevation and the GNI received by a flat surface facing the Sun are compared in Figure 5. The period was chosen to illustrate the rules listed before. This graph shows the correlation between the elevation angle and the GNI. The GNI and the Sun's elevation follow the same evolution during the period except for the moments when the rules impose that irradiance is null (τ = 0) or equal to the global lunar irradiance (τ = 1).

Finally, to calculate the plane of array (POA) irradiance POA_Irr, it is necessary to consider the incidence angle θ_S. When a flat surface is oriented towards the Sun (θ_S between 0° and 90°), POA_Irr is equal to GNI *cos(θ_S). When the Sun irradiates the backside of the plane (θ_S larger than 90°), POA_Irr is null (supposing that the surface is not a bifacial module). This leads to the last rule concerning the irradiance: if the absolute value of the incidence angle is superior to 90°, the irradiance received is null, POA_Irr = 0.

To summarize, the POA_Irr received by a flat surface at an instant t is:

Fig. 5 Evolution of the elevation angle of the Sun and the irradiance received when the Sun is visible at the point of interest with the latitude 89° from 22/01/2023 to 22/03/2023.

2.5 Solar energy calculations

The main algorithm can compute irradiance at any instant t. To compute the solar energy (in kWh) received by a flat surface for a given period, POA_Irr must be integrated over the period to obtain the irradiation (in kWh/m²) and then multiplied by the area of the surface. We used a time step of 1 h for the integration and a flat surface of 1 m² area.

We calculated the solar energy received considering four photovoltaic system installation modes. (1) A fixed mode meaning that the azimuth and tilt angles of the module are fixed. (2) A one-axis vertical tracker mode that considers the module elevation angle fixed and the module azimuth angle equal to the Sun azimuth angle. (3) A one-axis horizontal tracker mode where the module tilt follows the course of the Sun, while the azimuth is fixed. In this case, the controllable input of the one-axis horizontal tracker is the orientation of the horizontal axis. The tilt angle is represented by the absolute value of a special angle for this type of tracker called the “rotation angle” R, computed thanks to γ_Sun and β_Sun using [28]:

(14)

where γ_a the horizontal axis azimuth angle of the oriented surface.

In this tracker mode, γ_Surface takes only two values: γ_a ± 90° using:

(15)

(4) A two-axis tracker mode that assumes that both the azimuth and elevation angles of the surface are equal to the Sun's azimuth and 90° minus the Sun's elevation angles, respectively. In this case, POA_Irr is always considered equal to the GNI.

For each position and each tracker mode except the two-axis tracker, we computed the best angles to maximize the solar energy received. To find the best tilt and azimuth angles of a fixed surface, we first computed the solar energy every 2 degrees between −180° and 180° for azimuths and between 0° and 90° for elevations and then we computed the solar energy every 0.1 degrees around the maximum found before. For vertical and horizontal axis tracking systems, one of the two axes evolves along the day and the other is fixed. The optimal fixed angle is computed by the same method as for one angle of a fixed surface.

2.6 Localization of the points of interest for the computations

We consider 10 locations (see Tab. 1) on the Moon. The six first locations have a latitude ranging from 0° to 89°, while the longitude and altitude are fixed at 20° and 1,000 m, respectively. These locations will be used to compare the installation modes, depending on the latitude of the location. The last four locations are chosen as points of interest for future lunar installations [4,5,9]. The illumination fraction of these positions, which represents the fraction of time spent with a visible Sun over the total time spent for a defined period, has been estimated in the articles to be superior to 80% of the time on the period estimated, generally over an 18.6-year cycle. We verified it by computing the ratio of the number of hours when the elevation is larger than the limit angle minus the semi-apparent angle of the Sun on the number of hours in the lunar years between 2012 and 2031 for each location. The average illumination over the period and the associated standard deviation are presented in Table 1. Note that without considering the Moon's topography, the two points A1 and B1 are always exposed to Sun's rays with an illumination fraction of 100%.

At each location, we calculated the solar energy received by the 1 m² flat surface with the 4 installation modes for the lunar years between 2012 and 2031 included. Lunar years are shorter than terrestrial years for about 11 days and the start and the end are not fixed in the Gregorian calendar. Their duration is not fixed either. Between 2012 and 2031, the years lasted 8496, 8520, or 9216 h. The number of hours of the period has a non-negligible impact and the annual solar energies received are higher for the longest years. Consequently, we decided to take the average maximum solar energy values for one lunar year to compare the positions and the installation modes on their electrical potential on the Moon.

In the considered period, we covered an entire 18.6-year cycle. To test if the results can be translated to other cycles, we made the same computations for 20-year cycles between 1993 and 2031 and we found very similar results with on average less than 1% difference on solar energy received.

Finally, we also computed the limit angle α_limit defined in Section 2.4 and summarized the results in Table 1.

2.7 Implementation

The calculations and the code have been programmed in the Python language. The algorithm is available in open access via GitLab using the URL: https://gitlab.com/3it-cpv-public/pv-moon. The computation of the optimum angles of the four modes for all the positions and for all the lunar years between 2012 and 2031 took half a day with our laptop (Intel 10th Gen i7-1065G7, 1.30GHz clock frequency, 8 Go of random access memory). The code was also used to investigate tolerance angle around the best results as part of a master report chapter [29].

3 Validation of the algorithm

To validate the algorithm, we compared hour by hour for the lunar year 2022 the azimuths and the elevations angles calculated with our algorithm to the results given by the HORIZONS system [10], for the points of interest listed in Table 1 when the Sun's elevation is below 85°. Indeed, the notion of azimuth is vague when the Sun elevation is close to 90°. The comparison, summarized in Table A1 of Appendix A, shows a maximum difference of 0.031° for the elevation and 0.065° for the azimuth. This validates the Earth-Moon-Sun system and the lunar sky blocks.

4 Results for points with latitude between 0° and 89° and a longitude of 20°

Figure 6 presents the average annual solar energy received for the 4 installation modes (fixed, vertical or horizontal one-axis tracker, two-axis tracker) for the 6 locations on the Moon with the same longitude 20° and altitude 1,000 m. The only variable factor is the latitude, between 0° and 89°.

Fig. 6 Average maximum solar energy received in one lunar year by the flat surface considered between 2012 and 2031 for the 4 tracker modes and at the points of interest with longitude 20°.

4.1 Solar energy on a two-axis tracker (cross-shaped marker and dotted line)

Since a two-axis tracker always faces the Sun during the lunar day, it receives the highest solar energy compared to other installation modes. Therefore, this configuration is chosen as a reference. Figure 6 shows that the two-axis tracking mode leads to a solar energy of 6,176 kWh/yr on average for the location with latitude between 0 and 60°. When the latitude increases to 80° and 89° the solar energy increases to 6,758 and 10,988 kWh/yr, respectively. At first, the ∼1.8 times larger solar energy close to the pole can be attributed to the ∼1.8 times larger illumination fraction (see Tab. 1). We will see Section 5 that secondary order factors impact the solar energy received, especially near the poles.

4.2 Solar energy on a single vertical axis tracker (triangle-shaped marker and solid line)

Figure 6 shows also the solar energy received by a surface mounted on a vertical axis tracker. The difference with respect to the refence case (two-axis tracker) as well as the standard deviation are reported in Table 2.

We can see Figure 6 that the solar energy received by a vertical axis tracker increases with the latitude from 5,496 to 10,985 kWh/yr for a latitude from 0° to 89°, respectively. Reciprocally, the relative losses compared to a 2-axis tracker decrease from 10.4% at the equator to less than 0.1% close to the poles. This trend can be explained by a lower amplitude of the Sun's elevations when the latitude increases. Indeed, we report Figure 7 the range of elevation angles over an 18.6 yr cycle. In this figure, at the points with an absolute latitude close to 90°, including the point with a latitude of 89°, the Sun's elevation is equal to 0 +/− 3.5°. Therefore, the Sun's rays hit a vertical surface (tilt angle close to 90°) almost perpendicularly all along the period. The impact of elevation tracking is therefore negligible and the vertical-axis tracker behaves like the two-axis tracker.

To summarize, the use of a one-axis vertical tracker at the poles allows to obtain the same results as a two-axis tracker while being a simpler tracker, with only one mobile axis. At the equator and for latitudes under 40°, the losses are between 4 to 10% compared to a two-axis tracker, which is non-negligible.

Fig. 7 Elevation angle range of the Sun at the points of interest near the poles between the lunar years 2012 and 2031.

4.3 Solar energy on a single horizontal axis tracker (square-shaped marker and dashed line)

The solar energies calculated for a single axis horizontal tracker are presented Figure 6. The evolution of the solar energy for latitudes from 0° to 89° is inversed compared to the situation with a single vertical axis. When the latitude increases from 0 to 89°, the solar energy received decreases from 6,002 kWh/yr to 3,877 kWh/yr, respectively. It corresponds to 2.1 to 64.7% losses with respect to a two-axis tracker, respectively (see Tab. 2).

At the equator, due to the lunar equator's slight inclination of around 1.5° to the ecliptic (the orbital plane around the Sun), the Sun passes through all the elevations above a surface following the East-West direction and this path is similar every lunar day. Moving the tilt by the horizontal axe of rotation oriented in the North-South direction allows to follow the Sun continuously during all its path in the sky. That explains the similar solar energy received with a two-axis tracker at the equator.

To conclude, a horizontal axis tracker can yield a similar solar energy (within 98%) compared to a two-axis tracker at the equator and it is simpler with a single mobility axe.

4.4 Solar energy on a fixed surface with optimum tilt and azimuth angles (circle-shaped marker and dash-dot line)

Finally, we considered a fixed surface with optimal tilt and azimuth angles. For every latitudes considered from 0° to 89°, the total solar energy received by a fixed surface is equal to 3,822 +/− 0.5 kWh/yr. As the longitude of these points is the same, the Sun's path is similar with only an elevation difference. Even if the point with the latitude 89° has 95% of illumination fraction, as the surface is fixed and the Sun passes through all the azimuths from −180° to 180°, the Sun's rays hit the back of the surface as much as the surface itself. We computed the average incidence angle of the Sun on the surface, and we found that it is near 90° for all the points considered from latitude 0° to 89°. Thus, half the time, the Sun, whether visible or not, is behind the surface for every latitude considered. Considering monofacial surface for our study, the results for the point with the latitude 89° are the same than the other points.

Finally, if the solar energies received are globally similar for all the points of interest with the same longitude and altitude, the fixed surface has the highest losses of all the installation modes compared to the two-axis tracker. For points with a latitude below 80°, the losses are around 40% and for the point with the latitude of 89°, the losses are 65.2%.

5 Results for points of interests Shackleton crater, A1, B1 and North Pole

The four points Shackleton crater, A1, B1 and North Pole have maximum solar energies received for a two-axis tracker between 10,329 kWh/yr for the point North Pole and 11,889 kWh/yr for the point A1 (see Fig. 8).

At first, due to their latitudes near 89° and similar illumination fraction, a two-axis tracked surface at these points receives solar energy close to that received at the point with latitude 89°, see Section 4.1.

However, even if the two points A1 and B1 have both an illumination fraction of 100%, they do not receive the same exact solar energy. Furthermore, the solar energies considering only the illumination fraction (green curve) do not match with the solar energies considering our algorithm, see Figure 8. We can conclude that a secondary order factor impacts the solar energy received by a flat surface.

We attribute this second order factor to the rules introduced in Section 2.4 on Sun partial visibility. If the Sun is partially visible in the sky, the irradiance received is reduced which induces a lower solar energy for the period. The more the Sun is partially visible, the more the gap increases between the solar energy considering only the illumination fraction and the one of our algorithm.

Due to low elevations of the Sun near the poles (see Fig. 7), the Sun remains permanently around the horizon. Thus, it is more likely to be partially visible than at the other points, where the elevation fluctuates between greater limits, which amplifies the impact of this second order effect.

We reported Table 3 the average value of τ and the fraction of time the Sun is partially visible over the time it is visible. For the four points of interest, the Sun is partially visible 4.3% (point A1) to 7.8% (Shackleton crater) of the time with an average value of τ between 0.51% (Shackleton Crater and North Pole) and 0.81% (point A1). Point A1 has a shorter period of partial visibility of the Sun than the point B1 and the Sun is on average higher in the sky, represented through the average value of τ. The solar energy received at the point A1 is then higher than at the point B1.

To conclude, the combination of the partial visibility of the Sun and the average elevations of the Sun during that time explains why the points near the poles show lower results than would be expected if only the illumination fraction were considered.

To open the discussion, Figure 7 shows that the gap between the minimum and maximum elevation near the poles goes from (−1.97°,1.93°) to (−3.50°, 3.49°). This can be an indicator of the relevancy of the concentrator photovoltaics (CPV) with a single-axis tracker at the lunar poles, providing its acceptance angle covers those maximum values. Considering the almost absence of an atmosphere and therefore no scattering of solar rays and no diffuse light, a CPV system with an acceptance angle above +/−3.5° would be especially well suited for power generation on the Moon. The acceptance angle represents the maximum angle around the Sun's position at which a module produces 90% of its maximum energy.

Fig. 8 Comparison between the average maximum solar energy received by a surface following the Sun on two-axis for the points Shackleton crater, A1, B1 and North Pole and the average maximum solar energy received considering only the fraction illumination, for a lunar year.

6 Conclusion and discussion

This article compares the solar energy received by a flat surface using four types of tracking modes at different places on the Moon and for lunar years between 2012 and 2031, covering a cycle of most of the possible Sun-Moon relative positions.

First, we computed the elevations and azimuths of the Sun for every hour between 2012 and 2031. We used these angles to determine the incident angle between the Sun's rays and the surface and to calculate the solar energy received over a 20-year period. We performed this calculation for one-axis tracked surfaces (vertical or horizontal axis), for a two-axis tracked surface and for a fixed surface at optimal azimuth and tilt, for ten points of interest.

While the two-axis tracker presents the highest solar energies, similar solar energies are received near the poles by a vertical axis tracker, and near the equator with a horizontal axis tracker.

A fixed system would lose 37 to 64% of the solar energy compared to a two-axis tracker.

Finally, we discussed the phenomenon of partial visibility of the Sun which reduces the solar energy received, especially near the poles where the Sun remains around the horizon.

List of symbols

R_ES: Dimensionless Earth-to-Sun distance (in AU).

λ_ES: Apparent geocentric longitude of the Sun (in degrees).

R_EM: Dimensionless Earth-to-Moon distance (in AU).

λ_EM, ϕ_EM: Geocentric ecliptical longitude and latitude of the Moon, respectively (in degrees).

LSI: Lunar solar irradiance (in W/m²).

S₀: Solar constant, equals to 1367.8 W/m² in the article.

R_SM: Dimensionless Sun-to-Moon distance (in AU).

λ_SM, ϕ_SM: Heliocentric ecliptical longitude and latitude of the Moon, respectively (in degrees).

λ_Msp,ϕ_Msp: Selenographic longitude and latitude of the subsolar point, respectively (in degrees).

λ₀,ϕ₀: Selenographic longitude and latitude of the observer's position, respectively (in degrees).

γ_Sun,β_Sun: Azimuth and elevation of the Sun, respectively (in degrees).

R_Moon: Radius of the Moon (in km).

θ_S: Solar incidence angle (in degrees).

γ_Surface,β_Surface: Azimuth and tilt of the surface, respectively (in degrees).

α_limit: Limit angle of vision beyond horizon (in degrees).

h_point: Altitude of the point (in km).

τ: Factor that indicates if the Sun is visible or not in the sky (є [0;1]).

δ_S: Semi-apparent diameter of the Sun (in degrees).

GNI: Global normal irradiance (in W/m²).

POA_Irr: Plane of array irradiance (in W/m²).

R: Rotation angle (in degrees).

γ_a : Horizontal axis azimuth angle of the tracked surface (in degrees).

Acknowledgments

We wish to thank David Chuet for his support in the coding part. His advice has guided us to structure the code and make it clear.

We also wish to thank the Nanotechnologies Nanosystems Laboratory (LN2), which is a joint International Research Laboratory (IRL 3463) funded and co-operated in Canada by Université de Sherbrooke (UdS) and in France by CNRS as well as ECL, INSA Lyon, and Université Grenoble Alpes (UGA).

Funding

The research was funded by the following organizations:

• University of Sherbrooke, Canada.

• Laboratoire Nanotechnologies et Nanonsystèmes (LN2).

Conflicts of interest

We have a financial relationship with the organizations that sponsored the research, the University of Sherbrooke and the “Laboratoire Nanotechnologies et Nanonsystèmes” (LN2). We have full control of all primary data and we agree to allow the journal to review our data if requested.

Data availability statement

We have full control of all primary data and we agree to allow the journal to review our data if requested.

Author contribution statement

All authors have approved the submission of this manuscript. The results have not been previously published and are not being considered for publication in another journal. We confirm that this work is original and has not been published elsewhere, nor is it currently under consideration for publication elsewhere.

APPENDIX A:

Difference between the results from the algorithm and from the HORIZONS website

References

All Tables

All Figures

	Fig. 1 Block diagram of the main algorithm. Green inputs are controllable and the red output is the final output of the algorithm.
In the text

	Fig. 2 Angles and distances between bodies in the Earth-Moon-Sun system. The primary direction is the direction where the longitude is 0°, corresponding to vernal equinox [14].
In the text

	Fig. 3 Lunar global irradiance received in the lunar year 2022 (8520 h) by the Moon at the equator.
In the text

	Fig. 4 Angles used to compute the incidence angle on a flat surface (inspired by [20]).
In the text

	Fig. 5 Evolution of the elevation angle of the Sun and the irradiance received when the Sun is visible at the point of interest with the latitude 89° from 22/01/2023 to 22/03/2023.
In the text

	Fig. 6 Average maximum solar energy received in one lunar year by the flat surface considered between 2012 and 2031 for the 4 tracker modes and at the points of interest with longitude 20°.
In the text

	Fig. 7 Elevation angle range of the Sun at the points of interest near the poles between the lunar years 2012 and 2031.
In the text

	Fig. 8 Comparison between the average maximum solar energy received by a surface following the Sun on two-axis for the points Shackleton crater, A1, B1 and North Pole and the average maximum solar energy received considering only the fraction illumination, for a lunar year.
In the text