- Email: [email protected]
An algorithm for the computation of the solar position
Available online at www.sciencedirect.com
Solar Energy 82 (2008) 462–470 www.elsevier.com/locate/solener
An algorithm for the computation of the solar position Roberto Grena
*
ENEA – Centro Ricerche Casaccia, via Anguillarese 301, Roma, Italy Received 7 June 2007; received in revised form 26 September 2007; accepted 11 October 2007 Available online 5 November 2007 Communicated by: Associate Editor Lorin Vant-Hull
Abstract A new algorithm for the determination of the solar position is proposed. The importance of tracking accurately the position of the sun is evident when considering high concentration thermal systems. In the literature are found many fast algorithms for the determination of the sun position within 0.01°, and complex astronomical algorithms that allow a precision of 0.0003° but at the price of a large amount of calculation. The algorithm proposed in this work has a precision that is half-way between the two cases (maximal error 0.0027°), that should be sufficient for all the applications in solar engineering, with a computational cost comparable to the cost of the fast algorithms. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Sun position; Solar vector; Right ascension; Declination; Zenith; Azimuth
1. Introduction In this work an algorithm for the computation of the solar position with high precision (maximal error 0.0027° over the period 2003–2023) and of small complexity is proposed. This precision should be enough for all the necessities of solar engineering. Many fast algorithms for the calculation of the solar position, used in engineering application, can be found in the literature. They require a small computational effort, and their maximal error is usually greater than 0.01°. The Spencer formula (Spencer, 1971) has a maximal error greater than 0.25°; Pitman and Vant-Hull algorithm (Pitman and Vant-Hull, 1978) reduced the error to 0.02°; the Walraven algorithm (Walraven, 1978), with subsequent corrections and improvements (Walraven, 1979; Archer, 1980; Wilkinson, 1981, 1983; Muir, 1983), has an error of 0.013°; the Michalsky algorithm (Michalsky, 1988), used for comparison in the present work, has a maximal error of 0.011°; the last algorithm proposed, *
Corresponding author. Tel.: +39 06 3048 6856; fax: +39 06 3048 6779. E-mail address: [email protected]
0038-092X/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.solener.2007.10.001
the PSA algorithm (Blanco-Muriel et al., 2001) has a maximal error of 0.008°. All these algorithms work correctly for limited periods of time: e.g. 1950–2050 for the Michalsky algorithm, and 1995–2015 for the PSA algorithm. There are also high-precision astronomical algorithms, such as the numerical algorithm proposed by Meeus (1988), which has been reviewed in a form suitable for solar application by Reda and Andreas (2004) and is known as SPA (solar position algorithm). This algorithm has a maximal error smaller than 0.0003° on a very long period of time (2000 b.C.–6000 a.C.), but it requires a large amount of calculation. Such a precision can be necessary in some applications, such as the calibration of pyranometers, as shown by Reda and Andreas. In the present article a new algorithm is proposed, which has a maximal error half-way between the first and the second category, for the period 2003–2022, with a computational complexity only slightly increased with respect to the complexity of the fast algorithms. In Section 2, the input and output data of the algorithm are shown, with some remarks on the precision needed for the input data. In Section 3, the algorithm is described step
R. Grena / Solar Energy 82 (2008) 462–470
by step. In Section 4, the errors of the algorithm and a comparison among the proposed algorithm and two other algorithms (Michalsky and PSA) are discussed. The maximal error on the solar vector is 0.0027°, and the standard deviation is 0.001°. This precision is sufficient for all the engineering applications, but the amount of calculations is far smaller than in SPA algorithm. 2. Input and output data The algorithm requires the following input data: – Fractional UT of the day: UT (Universal Time, or Greenwich civil time) is measured in hours from the Greenwich midnight; the minutes and seconds must be converted into fractions of an hour; – Date (day D, months M, and years Y); – Difference D between UT and TT. TT (Terrestrial Time) is the time scale for ephemeres, and is not dependent on the rotation of the Earth. The difference D is D ¼ TT UT:
ð1Þ
Since the terrestrial rotation is irregularly slowing down, D increases, but not in a regular way; D can only be determined through measurements and extrapolations. For this algorithm, D must be given in seconds. – Longitude h and Latitude u (in radians); – Pressure P (in atm) and Temperature T (in °C). The algorithm computes the global coordinates of the sun (Right ascension a and declination d), and the local coordinates (Hour angle h, Zenith z and Azimuth C) corrected taking into account the atmospherical refraction. Some remarks on the precision required on the input data may prove useful. Errors on the position on the earth (h and u) and on the time UT (from which the angle of rotation of the earth is obtained) do not affect the global position of the sun (a and d), but the local position will be affected by these errors: it is therefore advisable that such data are given with a precision higher than the one attained by the algorithm calculation (an error smaller than 0.001° on h and u, and smaller than 0.2 s on the time UT, is recommended). Since 0.001° of latitude corresponds to 110 m on the earth surface, plants distributed on large surfaces may need different calculations of the sun position for the different parts of the plant, if the maximum precision is desired. Pressure and Temperature determine the air refraction, and are not critical except for very low elevation angles. The difference D is not critical: an error within 5 s on D does not affect significantly the algorithm precision. A simple extrapolation from last years’ measurements of D can be used without loss of precision for the period 2003–2023. The algorithm calculates first the angular position of the earth with respect to the sun, in the ecliptic plane (the heliocentric longitude); from this angle, and from
463
the inclination of the earth rotation axis, the position of the sun in the geocentric coordinates (the right ascension and the declination) can be calculated. The topocentric coordinates differ from the geocentric because they have their origin on the earth surface, and not on the earth centre; the difference must be considered (it affects the sun positions by many seconds of arc). Then, calculating the hour angle of the sun, the topocentric right ascension and declination can be converted in local coordinates (Zenith and Azimuth); a refraction correction is introduced. The algorithm considers the main effects that can affect the sun position by more than half a second of arc (Moon perturbation, nutation, difference between topocentric and geocentric coordinates), but all the perturbations are adapted to the period 2003–2023, strongly reducing the amount of calculations needed, especially the number of trigonometric functions. Empirical corrections are also introduced in the calculation of the heliocentric longitude, to sum up all the other small perturbations too complex to be considered one by one. 3. Procedure The steps of the algorithm are listed below. All the angles are in radians. 3.1. Time scales The time scales tG and t used in the computation are the Julian Day and the Ephemeris Julian Day1 respectively, shifted to make them start at noon, 1st January 2003. They are defined as follows. In this formula, if M is 1 or 2, M must be increased by 12 and Y must be reduced by 1. tG ¼ INTð365:25ðY 2000ÞÞ þ INTð30:6001ðM þ 1ÞÞ UT 1158:5; 24 D t ¼ tG þ : 86400 þDþ
ð2Þ
The INT function rounds its argument to the nearest integer towards 0 (i.e. INT(7.8) = 7; INT(6.6) = 6). 3.2. Calculation of the heliocentric longitude of the earth This is the most critical point, and the main source of errors. For convenience, it will be split into Section 3.2.1,3.2.2,3.2.3,3.2.4:
1 The standard JD and JDE (starting at noon, 1st January 4712 b.C.) are:
JDE ¼ t þ 2452640; JD ¼ JDt þ 2452640;
464
R. Grena / Solar Energy 82 (2008) 462–470
3.2.1. Linear increasing with annual oscillation R ¼ 1:72019e 2t 0:0563;
ð3Þ
Ly ¼ 1:74094 þ 1:7202768683e 2t þ 3:34118e 2 sin R
ahr ¼ 12modða; 2pÞ=p;
þ 3:488e 4 sin 2R:
ð12Þ
and the fractional part can be converted in minutes and seconds. But this is not necessary when carrying out the calculation.
3.2.2. Moon perturbation Lm ¼ 3:13e 5 sinð0:2127730t 0:585Þ:
Usually the right ascension is measured in hours (1 h = 15°), from 0 to 24; the right ascension in hours can be easily obtained as
ð4Þ
3.5.3. Declination d ¼ asinðsin e sin cÞ:
ð13Þ
3.2.3. Harmonic correction 3.6. Local hour angle of the sun
Lh ¼ 1:26e 5 sinð4:243e 3t þ 1:46Þ þ 2:35e 5 sinð1:0727e 2t þ 0:72Þ þ 2:76e 5
h ¼ 6:30038809903tG þ 4:8824623 þ 0:9174Dc þ h a:
sinð1:5799e 2t þ 2:35Þ þ 2:75e 5
ð14Þ
sinð2:1551e 2t 1:98Þ þ 1:26e 5 sinð3:1490e 2t 0:80Þ:
ð5Þ
3.2.4. Polynomial correction t2 ¼ 0:001t;
Da ¼ 4:26e 5 cos u sin h:
ð15Þ
ð6Þ
Lp ¼ ððð2:30796e 7t2 þ 3:7976e 6Þt2 2:0458e 5Þt2 þ 3:976e
3.7. Parallax correction to right ascension
5Þt22 :
3.8. Topocentric sun coordinates 3.8.1. Topocentric right ascension
The time scale t2 is introduced in order to have more homogeneous quantities in the products within the polynomial, avoiding too rough rounding approximation. The longitude is the sum of these four terms:
at ¼ a þ Da:
L ¼ Ly þ Lm þ Lh þ Lp :
dt ¼ d 4:26e 5ðsin u d cos uÞ:
ð7Þ
3.3. Correction to geocentric longitude due to nutation Dc ¼ 8:33e 5 sinð9:252e 4t 1:173Þ:
ð16Þ
3.8.2. Topocentric declination ð17Þ
3.8.3. Topocentric hour angle ð8Þ
3.4. Earth axis inclination
ht ¼ h Da;
ð18Þ
cht ¼ cos h þ Da sin h ðapproximate cosine of ht Þ; sht ¼ sin h Da cos h ðapproximate sine of ht Þ:
ð19Þ ð20Þ
e ¼ 6:21e 9t þ 0:409086 þ 4:46e 5 sinð9:252e 4t þ 0:397Þ:
ð9Þ
3.5. Geocentric global solar coordinates
e0 ¼ asinðsin u sin dt þ cos u cos dt cht Þ:
3.5.1. Geocentric solar longitude c ¼ L þ p þ Dc 9:932e 5:
ð21Þ
3.10. Atmospheric refraction correction to the solar elevation ð10Þ De ¼ 0:084217P =½ð273 þ T Þ tanðe0 þ 0:0031376
3.5.2. Geocentric right ascension a ¼ atan2ðsin c cos e; cos cÞ:
3.9. Solar elevation angle, without refraction correction
=ðe0 þ 0:089186ÞÞ: ð11Þ
The function atan2 takes as input two quantities proportional to the sine and the cosine of the angle, and outputs a value in radians between –p and p.
ð22Þ
This formula is the same used by Reda and Andreas (2004), rewritten in radians. Since it is independent from the other steps of the algorithm, it can be replaced by any formula for the refraction correction that includes all the
R. Grena / Solar Energy 82 (2008) 462–470
465
Table 1 Errors on the main quantities calculated by the algorithm proposed in this work (ENEA), in seconds of arc (00 ), compared with the two most recent fast algorithms found in literature ENEA 00
2.6 0.900 2.000 3.700 3.600
a d z C Solar vector
PSA Range
r
00
Michalsky
r
Range
r
Range
6.400 10.600 9.400
[23.800 22.000 ] [92.000 86.600 ] [000 26.000 ]
11.200 13.000 14.100
[39.600 20.400 ] [114.200 80.600 ] [000 40.000 ]
00
[7.5 8.2 ] [2.700 3.100 ] [7.000 9.300 ] [29.100 24.100 ] [000 9.800 ]
Table 2 Computational cost required by the algorithms
PSA Michalsky ENEA SPA
Additions/subtractions
Products/divisions
Direct trigonometric functions
Inverse trigonometric functions
35 28 52 >1000
36 25 54 >1300
16 19 21 >300
4 4 4 7
parameters one chooses to consider, or even ignored if the refraction is considered negligible. 3.11. Local topocentric sun coordinates 3.11.1. Zenith z ¼ p=2 e0 De:
ð23Þ
3.11.2. Azimuth C ¼ atan2ðsht ; cht sin u tan dt cos uÞ:
ð24Þ
The azimuth obtained varies from –p to p, C = 0 towards the south direction, and the azimuth is positive in the western hemisphere. 4. Comparison with other algorithms The errors of this algorithm have been calculated for the period of validity (2003–2023), for a location at 0° longitude, 40° latitude; the errors have been computed as the difference between the outputs of the present algorithm and of the SPA algorithm, which qualifies as a good choice for a term of comparison, since its errors are far smaller than all the errors of the fast algorithms. In this calculation, the refraction formula has been applied when e0 > 0, otherwise the elevation considered is the uncorrected one. The results have been compared to the errors of the two most recent fast algorithms in the literature (Michalsky and PSA), reported in Table 72
2 From that table, the quadraticffi error can be calculated as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Average error2 þ Standard deviation2
in Blanco-Muriel et al. (2001), and obtained for the period 1999–2015 in a similar location (37°6 0 200 N, 2°21 0 3600 E). The comparison is shown in Table 1, where are reported the quadratic error (defined as the square root of the mean squared error) and the error range, in seconds of arc, for all the main quantities of interest. The most significant information is given by the error on the solar vector, which is the angular distance between the real and the calculated position. (Here the ‘‘real’’ position is the one calculated with SPA). Since the errors are small, the error on the solar vector is given by the formula qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi DV ¼ Dz2 þ ðDC sin zÞ2 ð25Þ where Dz is the error on the Zenith angle z, and DC is the error on the Azimuth. Considering both the quadratic error and the maximal error, the proposed algorithm reduces the error on the solar vector by 60% with respect to the PSA algorithm, and by 75% with respect to the Michalsky algorithm. The computational complexity of the algorithms has also been compared: from Table 2 one can see that the new algorithm only requires a slight increase in computational cost with respect to the fast algorithms, while being much faster than SPA algorithm. The time of calculation depends on the machine and on the programming language used; a rough comparison can be obtained by simply comparing the number of times each algorithm calls a trigonometric function (either direct and inverse), since such functions usually have a computational cost much bigger that of all the other operations. In this comparison, the PSA algorithm turns out to be a 25% faster than the proposed algorithm, but it does not calculate the refraction correction. The Michalsky algorithm is
466
R. Grena / Solar Energy 82 (2008) 462–470
Fig. 1. Error in right ascension.
Fig. 2. Error in declination.
Fig. 3. Error in Zenith distance.
Fig. 4. Error in Azimuth.
Fig. 5. Error in solar vector.
Fig. 6. Distribution of errors in right ascension.
R. Grena / Solar Energy 82 (2008) 462–470
Fig. 7. Distribution of errors in declination.
467
Fig. 10. Distribution of errors in solar vector.
a 10% faster than the proposed algorithm. All these three algorithms are more than 15 times faster than the SPA algorithm. Figs. 1–5 show the temporal variation of the errors, and Figs. 6–10 show their distribution.3
5. Conclusions
Fig. 8. Distribution of errors in Zenith distance.
The presented algorithm allows one to reach a very good precision, sufficient to all the needs of solar engineering. The sun position is obtained with an acceptable amount of calculations, comparable with the computational effort required by the algorithms usually adopted in the solar engineering, and much smaller than the one required by the astronomical algorithms. The period of validity is limited (20 years), but this is not a problem for engineering applications. The algorithm can be periodically modified so as to adapt it to future periods. The algorithm can be successfully employed in the control of high precision tracking systems, or for simulations that require a large number of sun position calculations.
Acknowledgement Thanks to Antonio De Luca and to all the SOLTERM Group of ENEA for useful discussions and support.
Fig. 9. Distribution of errors in Azimuth.
3 In the distribution of the declination errors, there is a small bias (0.22 s of arc); this is probably due to the fact that the sun is above the celestial equator for a longer time than below, and the error is reasonably the same if the declination changes sign. This bias could be easily eliminated, but the symmetry would be lost. Moreover, such a small bias does not affect significantly the precision of the algorithm.
468
R. Grena / Solar Energy 82 (2008) 462–470
Appendix A. C++ implementation of the algorithm
R. Grena / Solar Energy 82 (2008) 462–470
469
470
R. Grena / Solar Energy 82 (2008) 462–470
References Archer, C.B., 1980. Comments on ‘Calculating the position of the Sun’. Solar Energy 25, 91. Blanco-Muriel, M., Alarcon-Padilla, D.C., Lopea-Moratalla, T., Lara-Coira, M., 2001. Computing the solar vector. Solar Energy 70 (5), 431–441. Meeus, J., 1988. Astronomical Algorithms, second ed. Willmann-Bell Inc.. Michalsky, J.J., 1988. The astronomical Almanac’s algorithm for approximate solar position (1950–2050). Solar Energy 40 (3), 227–235. Muir, L.R., 1983. Comments on ‘The effect of atmospheric refraction on the solar azimuth’. Solar Energy 30, 295. Pitman, C.L., Vant-Hull, L.L., 1978. Errors in locating the Sun and their effect on solar intensity predictions. In: Meeting of the American
Section of the International Solar Energy Society, Denver, 28 August, pp. 701–706. Reda, I., Andreas, A., 2004. Solar position algorithm for solar radiation applications. Solar Energy 76 (5), 577–589. Spencer, J.W., 1971. Fourier series representation of the position of the Sun. Search 2 (5). Walraven, R., 1978. Calculating the position of the Sun. Solar Energy 20, 393–397. Walraven, R., 1979 Erratum. Solar Energy 22, 195. Wilkinson, B.J., 1981. An improved FORTRAN program for the rapid calculation of the solar position. Solar Energy 27, 67–68. Wilkinson, B.J., 1983. The effect of atmospheric refraction on the solar azimuth. Solar Energy 30, 295.
Solar Energy 82 (2008) 462–470 www.elsevier.com/locate/solener
An algorithm for the computation of the solar position Roberto Grena
*
ENEA – Centro Ricerche Casaccia, via Anguillarese 301, Roma, Italy Received 7 June 2007; received in revised form 26 September 2007; accepted 11 October 2007 Available online 5 November 2007 Communicated by: Associate Editor Lorin Vant-Hull
Abstract A new algorithm for the determination of the solar position is proposed. The importance of tracking accurately the position of the sun is evident when considering high concentration thermal systems. In the literature are found many fast algorithms for the determination of the sun position within 0.01°, and complex astronomical algorithms that allow a precision of 0.0003° but at the price of a large amount of calculation. The algorithm proposed in this work has a precision that is half-way between the two cases (maximal error 0.0027°), that should be sufficient for all the applications in solar engineering, with a computational cost comparable to the cost of the fast algorithms. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Sun position; Solar vector; Right ascension; Declination; Zenith; Azimuth
1. Introduction In this work an algorithm for the computation of the solar position with high precision (maximal error 0.0027° over the period 2003–2023) and of small complexity is proposed. This precision should be enough for all the necessities of solar engineering. Many fast algorithms for the calculation of the solar position, used in engineering application, can be found in the literature. They require a small computational effort, and their maximal error is usually greater than 0.01°. The Spencer formula (Spencer, 1971) has a maximal error greater than 0.25°; Pitman and Vant-Hull algorithm (Pitman and Vant-Hull, 1978) reduced the error to 0.02°; the Walraven algorithm (Walraven, 1978), with subsequent corrections and improvements (Walraven, 1979; Archer, 1980; Wilkinson, 1981, 1983; Muir, 1983), has an error of 0.013°; the Michalsky algorithm (Michalsky, 1988), used for comparison in the present work, has a maximal error of 0.011°; the last algorithm proposed, *
Corresponding author. Tel.: +39 06 3048 6856; fax: +39 06 3048 6779. E-mail address: [email protected]
0038-092X/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.solener.2007.10.001
the PSA algorithm (Blanco-Muriel et al., 2001) has a maximal error of 0.008°. All these algorithms work correctly for limited periods of time: e.g. 1950–2050 for the Michalsky algorithm, and 1995–2015 for the PSA algorithm. There are also high-precision astronomical algorithms, such as the numerical algorithm proposed by Meeus (1988), which has been reviewed in a form suitable for solar application by Reda and Andreas (2004) and is known as SPA (solar position algorithm). This algorithm has a maximal error smaller than 0.0003° on a very long period of time (2000 b.C.–6000 a.C.), but it requires a large amount of calculation. Such a precision can be necessary in some applications, such as the calibration of pyranometers, as shown by Reda and Andreas. In the present article a new algorithm is proposed, which has a maximal error half-way between the first and the second category, for the period 2003–2022, with a computational complexity only slightly increased with respect to the complexity of the fast algorithms. In Section 2, the input and output data of the algorithm are shown, with some remarks on the precision needed for the input data. In Section 3, the algorithm is described step
R. Grena / Solar Energy 82 (2008) 462–470
by step. In Section 4, the errors of the algorithm and a comparison among the proposed algorithm and two other algorithms (Michalsky and PSA) are discussed. The maximal error on the solar vector is 0.0027°, and the standard deviation is 0.001°. This precision is sufficient for all the engineering applications, but the amount of calculations is far smaller than in SPA algorithm. 2. Input and output data The algorithm requires the following input data: – Fractional UT of the day: UT (Universal Time, or Greenwich civil time) is measured in hours from the Greenwich midnight; the minutes and seconds must be converted into fractions of an hour; – Date (day D, months M, and years Y); – Difference D between UT and TT. TT (Terrestrial Time) is the time scale for ephemeres, and is not dependent on the rotation of the Earth. The difference D is D ¼ TT UT:
ð1Þ
Since the terrestrial rotation is irregularly slowing down, D increases, but not in a regular way; D can only be determined through measurements and extrapolations. For this algorithm, D must be given in seconds. – Longitude h and Latitude u (in radians); – Pressure P (in atm) and Temperature T (in °C). The algorithm computes the global coordinates of the sun (Right ascension a and declination d), and the local coordinates (Hour angle h, Zenith z and Azimuth C) corrected taking into account the atmospherical refraction. Some remarks on the precision required on the input data may prove useful. Errors on the position on the earth (h and u) and on the time UT (from which the angle of rotation of the earth is obtained) do not affect the global position of the sun (a and d), but the local position will be affected by these errors: it is therefore advisable that such data are given with a precision higher than the one attained by the algorithm calculation (an error smaller than 0.001° on h and u, and smaller than 0.2 s on the time UT, is recommended). Since 0.001° of latitude corresponds to 110 m on the earth surface, plants distributed on large surfaces may need different calculations of the sun position for the different parts of the plant, if the maximum precision is desired. Pressure and Temperature determine the air refraction, and are not critical except for very low elevation angles. The difference D is not critical: an error within 5 s on D does not affect significantly the algorithm precision. A simple extrapolation from last years’ measurements of D can be used without loss of precision for the period 2003–2023. The algorithm calculates first the angular position of the earth with respect to the sun, in the ecliptic plane (the heliocentric longitude); from this angle, and from
463
the inclination of the earth rotation axis, the position of the sun in the geocentric coordinates (the right ascension and the declination) can be calculated. The topocentric coordinates differ from the geocentric because they have their origin on the earth surface, and not on the earth centre; the difference must be considered (it affects the sun positions by many seconds of arc). Then, calculating the hour angle of the sun, the topocentric right ascension and declination can be converted in local coordinates (Zenith and Azimuth); a refraction correction is introduced. The algorithm considers the main effects that can affect the sun position by more than half a second of arc (Moon perturbation, nutation, difference between topocentric and geocentric coordinates), but all the perturbations are adapted to the period 2003–2023, strongly reducing the amount of calculations needed, especially the number of trigonometric functions. Empirical corrections are also introduced in the calculation of the heliocentric longitude, to sum up all the other small perturbations too complex to be considered one by one. 3. Procedure The steps of the algorithm are listed below. All the angles are in radians. 3.1. Time scales The time scales tG and t used in the computation are the Julian Day and the Ephemeris Julian Day1 respectively, shifted to make them start at noon, 1st January 2003. They are defined as follows. In this formula, if M is 1 or 2, M must be increased by 12 and Y must be reduced by 1. tG ¼ INTð365:25ðY 2000ÞÞ þ INTð30:6001ðM þ 1ÞÞ UT 1158:5; 24 D t ¼ tG þ : 86400 þDþ
ð2Þ
The INT function rounds its argument to the nearest integer towards 0 (i.e. INT(7.8) = 7; INT(6.6) = 6). 3.2. Calculation of the heliocentric longitude of the earth This is the most critical point, and the main source of errors. For convenience, it will be split into Section 3.2.1,3.2.2,3.2.3,3.2.4:
1 The standard JD and JDE (starting at noon, 1st January 4712 b.C.) are:
JDE ¼ t þ 2452640; JD ¼ JDt þ 2452640;
464
R. Grena / Solar Energy 82 (2008) 462–470
3.2.1. Linear increasing with annual oscillation R ¼ 1:72019e 2t 0:0563;
ð3Þ
Ly ¼ 1:74094 þ 1:7202768683e 2t þ 3:34118e 2 sin R
ahr ¼ 12modða; 2pÞ=p;
þ 3:488e 4 sin 2R:
ð12Þ
and the fractional part can be converted in minutes and seconds. But this is not necessary when carrying out the calculation.
3.2.2. Moon perturbation Lm ¼ 3:13e 5 sinð0:2127730t 0:585Þ:
Usually the right ascension is measured in hours (1 h = 15°), from 0 to 24; the right ascension in hours can be easily obtained as
ð4Þ
3.5.3. Declination d ¼ asinðsin e sin cÞ:
ð13Þ
3.2.3. Harmonic correction 3.6. Local hour angle of the sun
Lh ¼ 1:26e 5 sinð4:243e 3t þ 1:46Þ þ 2:35e 5 sinð1:0727e 2t þ 0:72Þ þ 2:76e 5
h ¼ 6:30038809903tG þ 4:8824623 þ 0:9174Dc þ h a:
sinð1:5799e 2t þ 2:35Þ þ 2:75e 5
ð14Þ
sinð2:1551e 2t 1:98Þ þ 1:26e 5 sinð3:1490e 2t 0:80Þ:
ð5Þ
3.2.4. Polynomial correction t2 ¼ 0:001t;
Da ¼ 4:26e 5 cos u sin h:
ð15Þ
ð6Þ
Lp ¼ ððð2:30796e 7t2 þ 3:7976e 6Þt2 2:0458e 5Þt2 þ 3:976e
3.7. Parallax correction to right ascension
5Þt22 :
3.8. Topocentric sun coordinates 3.8.1. Topocentric right ascension
The time scale t2 is introduced in order to have more homogeneous quantities in the products within the polynomial, avoiding too rough rounding approximation. The longitude is the sum of these four terms:
at ¼ a þ Da:
L ¼ Ly þ Lm þ Lh þ Lp :
dt ¼ d 4:26e 5ðsin u d cos uÞ:
ð7Þ
3.3. Correction to geocentric longitude due to nutation Dc ¼ 8:33e 5 sinð9:252e 4t 1:173Þ:
ð16Þ
3.8.2. Topocentric declination ð17Þ
3.8.3. Topocentric hour angle ð8Þ
3.4. Earth axis inclination
ht ¼ h Da;
ð18Þ
cht ¼ cos h þ Da sin h ðapproximate cosine of ht Þ; sht ¼ sin h Da cos h ðapproximate sine of ht Þ:
ð19Þ ð20Þ
e ¼ 6:21e 9t þ 0:409086 þ 4:46e 5 sinð9:252e 4t þ 0:397Þ:
ð9Þ
3.5. Geocentric global solar coordinates
e0 ¼ asinðsin u sin dt þ cos u cos dt cht Þ:
3.5.1. Geocentric solar longitude c ¼ L þ p þ Dc 9:932e 5:
ð21Þ
3.10. Atmospheric refraction correction to the solar elevation ð10Þ De ¼ 0:084217P =½ð273 þ T Þ tanðe0 þ 0:0031376
3.5.2. Geocentric right ascension a ¼ atan2ðsin c cos e; cos cÞ:
3.9. Solar elevation angle, without refraction correction
=ðe0 þ 0:089186ÞÞ: ð11Þ
The function atan2 takes as input two quantities proportional to the sine and the cosine of the angle, and outputs a value in radians between –p and p.
ð22Þ
This formula is the same used by Reda and Andreas (2004), rewritten in radians. Since it is independent from the other steps of the algorithm, it can be replaced by any formula for the refraction correction that includes all the
R. Grena / Solar Energy 82 (2008) 462–470
465
Table 1 Errors on the main quantities calculated by the algorithm proposed in this work (ENEA), in seconds of arc (00 ), compared with the two most recent fast algorithms found in literature ENEA 00
2.6 0.900 2.000 3.700 3.600
a d z C Solar vector
PSA Range
r
00
Michalsky
r
Range
r
Range
6.400 10.600 9.400
[23.800 22.000 ] [92.000 86.600 ] [000 26.000 ]
11.200 13.000 14.100
[39.600 20.400 ] [114.200 80.600 ] [000 40.000 ]
00
[7.5 8.2 ] [2.700 3.100 ] [7.000 9.300 ] [29.100 24.100 ] [000 9.800 ]
Table 2 Computational cost required by the algorithms
PSA Michalsky ENEA SPA
Additions/subtractions
Products/divisions
Direct trigonometric functions
Inverse trigonometric functions
35 28 52 >1000
36 25 54 >1300
16 19 21 >300
4 4 4 7
parameters one chooses to consider, or even ignored if the refraction is considered negligible. 3.11. Local topocentric sun coordinates 3.11.1. Zenith z ¼ p=2 e0 De:
ð23Þ
3.11.2. Azimuth C ¼ atan2ðsht ; cht sin u tan dt cos uÞ:
ð24Þ
The azimuth obtained varies from –p to p, C = 0 towards the south direction, and the azimuth is positive in the western hemisphere. 4. Comparison with other algorithms The errors of this algorithm have been calculated for the period of validity (2003–2023), for a location at 0° longitude, 40° latitude; the errors have been computed as the difference between the outputs of the present algorithm and of the SPA algorithm, which qualifies as a good choice for a term of comparison, since its errors are far smaller than all the errors of the fast algorithms. In this calculation, the refraction formula has been applied when e0 > 0, otherwise the elevation considered is the uncorrected one. The results have been compared to the errors of the two most recent fast algorithms in the literature (Michalsky and PSA), reported in Table 72
2 From that table, the quadraticffi error can be calculated as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Average error2 þ Standard deviation2
in Blanco-Muriel et al. (2001), and obtained for the period 1999–2015 in a similar location (37°6 0 200 N, 2°21 0 3600 E). The comparison is shown in Table 1, where are reported the quadratic error (defined as the square root of the mean squared error) and the error range, in seconds of arc, for all the main quantities of interest. The most significant information is given by the error on the solar vector, which is the angular distance between the real and the calculated position. (Here the ‘‘real’’ position is the one calculated with SPA). Since the errors are small, the error on the solar vector is given by the formula qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi DV ¼ Dz2 þ ðDC sin zÞ2 ð25Þ where Dz is the error on the Zenith angle z, and DC is the error on the Azimuth. Considering both the quadratic error and the maximal error, the proposed algorithm reduces the error on the solar vector by 60% with respect to the PSA algorithm, and by 75% with respect to the Michalsky algorithm. The computational complexity of the algorithms has also been compared: from Table 2 one can see that the new algorithm only requires a slight increase in computational cost with respect to the fast algorithms, while being much faster than SPA algorithm. The time of calculation depends on the machine and on the programming language used; a rough comparison can be obtained by simply comparing the number of times each algorithm calls a trigonometric function (either direct and inverse), since such functions usually have a computational cost much bigger that of all the other operations. In this comparison, the PSA algorithm turns out to be a 25% faster than the proposed algorithm, but it does not calculate the refraction correction. The Michalsky algorithm is
466
R. Grena / Solar Energy 82 (2008) 462–470
Fig. 1. Error in right ascension.
Fig. 2. Error in declination.
Fig. 3. Error in Zenith distance.
Fig. 4. Error in Azimuth.
Fig. 5. Error in solar vector.
Fig. 6. Distribution of errors in right ascension.
R. Grena / Solar Energy 82 (2008) 462–470
Fig. 7. Distribution of errors in declination.
467
Fig. 10. Distribution of errors in solar vector.
a 10% faster than the proposed algorithm. All these three algorithms are more than 15 times faster than the SPA algorithm. Figs. 1–5 show the temporal variation of the errors, and Figs. 6–10 show their distribution.3
5. Conclusions
Fig. 8. Distribution of errors in Zenith distance.
The presented algorithm allows one to reach a very good precision, sufficient to all the needs of solar engineering. The sun position is obtained with an acceptable amount of calculations, comparable with the computational effort required by the algorithms usually adopted in the solar engineering, and much smaller than the one required by the astronomical algorithms. The period of validity is limited (20 years), but this is not a problem for engineering applications. The algorithm can be periodically modified so as to adapt it to future periods. The algorithm can be successfully employed in the control of high precision tracking systems, or for simulations that require a large number of sun position calculations.
Acknowledgement Thanks to Antonio De Luca and to all the SOLTERM Group of ENEA for useful discussions and support.
Fig. 9. Distribution of errors in Azimuth.
3 In the distribution of the declination errors, there is a small bias (0.22 s of arc); this is probably due to the fact that the sun is above the celestial equator for a longer time than below, and the error is reasonably the same if the declination changes sign. This bias could be easily eliminated, but the symmetry would be lost. Moreover, such a small bias does not affect significantly the precision of the algorithm.
468
R. Grena / Solar Energy 82 (2008) 462–470
Appendix A. C++ implementation of the algorithm
R. Grena / Solar Energy 82 (2008) 462–470
469
470
R. Grena / Solar Energy 82 (2008) 462–470
References Archer, C.B., 1980. Comments on ‘Calculating the position of the Sun’. Solar Energy 25, 91. Blanco-Muriel, M., Alarcon-Padilla, D.C., Lopea-Moratalla, T., Lara-Coira, M., 2001. Computing the solar vector. Solar Energy 70 (5), 431–441. Meeus, J., 1988. Astronomical Algorithms, second ed. Willmann-Bell Inc.. Michalsky, J.J., 1988. The astronomical Almanac’s algorithm for approximate solar position (1950–2050). Solar Energy 40 (3), 227–235. Muir, L.R., 1983. Comments on ‘The effect of atmospheric refraction on the solar azimuth’. Solar Energy 30, 295. Pitman, C.L., Vant-Hull, L.L., 1978. Errors in locating the Sun and their effect on solar intensity predictions. In: Meeting of the American
Section of the International Solar Energy Society, Denver, 28 August, pp. 701–706. Reda, I., Andreas, A., 2004. Solar position algorithm for solar radiation applications. Solar Energy 76 (5), 577–589. Spencer, J.W., 1971. Fourier series representation of the position of the Sun. Search 2 (5). Walraven, R., 1978. Calculating the position of the Sun. Solar Energy 20, 393–397. Walraven, R., 1979 Erratum. Solar Energy 22, 195. Wilkinson, B.J., 1981. An improved FORTRAN program for the rapid calculation of the solar position. Solar Energy 27, 67–68. Wilkinson, B.J., 1983. The effect of atmospheric refraction on the solar azimuth. Solar Energy 30, 295.