High-Precision Angle Prediction of Sun for Space Remote Sensing Instrument

Abstract

In order to grasp the timing of sun calibration in advance, this paper introduces a high-precision method to predict the solar angle by using the current broadcast time and orbital instantaneous root of the satellite platform. By calculating the sun’s apparent right ascension and apparent declination, the conversion matrix from the geocentric inertial coordinate system to the orbital coordinate system, and the satellite attitude correction matrix, the sun vector in the satellite body coordinate system is obtained. This method is used to predict the sun angle of a sun synchronous orbit in the satellite coordinate system, and the prediction results are compared with the STK simulation results. The results show that the sun angle prediction error of this method is less than ±0.003°. It can meet the requirements of on-orbit solar calibration. The main error sources in the prediction method are analysed.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

With the further development of space remote sensing instruments and the continuous improvement of quantitative application requirements of remote sensing products, the high-precision calibration of space remote sensing instruments has become increasingly important and necessary. Among the proposed calibration methods [1], compared with laboratory calibration before launch and on-orbit alternative calibration, onboard calibration can directly reflect the actual situation of the remote sensing instrument after entering orbit. It is an important means to monitor the long-term performance changes of remote sensing instruments, and it is also the basis for quantitatively obtaining on orbit performance information of remote sensing instruments and correcting on orbit telemetric data [2].

In-orbit calibration system is usually divided into in-orbit internal calibration system and in-orbit external calibration system. The in-orbit calibration system can only monitor part of the optical system and increase the combined uncertainty of the final results [3, 4]. The sun is a stable radiation source. The method of using the sun as the standard light source for onboard calibration has high calibration accuracy and can make up for the shortcomings of the built-in standard light source which is easy to decay for a long time [5–8].

Because the incident angle of the sun changes with the seasons and the orbit parameters of the satellite, and the satellite attitude adjustment may also be needed during the calibration, it is necessary to grasp the timing of the sun calibration in advance, that is, the angle prediction of the sun. Usually, space remote sensing instruments use the time of the Julian day and orbit parameters to calculate the angle of the sun [9, 10]. Due to the limited computing resources of space remote sensing instruments, orbit parameters cannot be predicted or can only be predicted in a short time [11]. On the premise of considering both accuracy and simplification, this paper presents a high-precision solar angle prediction method for space remote sensing instruments.

2. Materials and Methods

The flow chart of the method for the angle prediction of the sun proposed in this paper is shown in Figure 1. The specific steps are as follows. (1) According to the current time broadcast by the current satellite platform ($t_0$) and the predicted time interval ($\Delta t$), the apparent right ascension and apparent declination of the sun at ($t_0+n\ast \Delta t$) are calculated. The time interval ($\Delta t$) represents the time difference between the previous and the next solar angle prediction, which can range from a few seconds to a few days according to the specific application. $n$ intervals mean that the current forecast time ($t_0+n\ast \Delta t$) has increased by $n$ time intervals compared with the start time ($t_0$). $n$ starts from 0. After the solar angle prediction at the current time, $n$ increases by 1 for the solar angle prediction at the next prediction time. (2) According to the instantaneous root of the satellite orbit at $t_0$ time broadcast by the satellite platform, the instantaneous root of the satellite orbit at $t_0+n\ast \Delta t$ is calculated, and then, the transformation matrix from the geocentric inertial coordinate system to the satellite local coordinate system is calculated. (3) According to the satellite attitude angle at $t_0$, the attitude correction matrix is calculated. (4) Using the results of the above steps, the solar angle at ($t_0+n\ast \Delta t$) time in the satellite body coordinate system is calculated until it reaches the end time of prediction ($t_1$).

[figure omitted; refer to PDF]

The vsop87 planetary theory consists of six tables of periodic coefficients. The vsop87D table can directly calculate the ecliptic longitude ($L$) and the ecliptic latitude ($B$) of the planets in the heliocentric coordinate system and the distance from the planet to the sun ($R$). The vsop87D table contains three parts of data, namely, the periodic term coefficient table for calculating the ecliptic longitude of planets in the heliocentric coordinate system ($L$ table), the periodic term coefficient table for calculating the ecliptic latitude of planets in the heliocentric coordinate system ($B$ table), and the periodic term coefficient table for calculating the distance between planets and the sun ($R$ table).

We use the date broadcast on the satellite platform to calculate the Julian day $J_D$. The formula for calculating Julian century $T$ is as follows. $$\begin{array}{}\text{(1)}& T=\frac{J_D-2451545.0}{365250}.\end{array}$$

Suppose that the three coefficients $A$, $B$, and $C$ are used to calculate the earth’s ecliptic longitude in the heliocentric coordinate system, then the expression of each periodic term $X$ is as follows. $$\begin{array}{}\text{(2)}& X=A\times \cos B+C\times T.\end{array}$$

$L_0$ is calculated by adding the items in L0 table, $L_1$ is calculated by adding the items in L1 table, and so on. The final ecliptic longitude in the heliocentric coordinate system is calculated according to Formula (3). $$\begin{array}{}\text{(3)}& L=\frac{L_0+L_1\times T+L_2\times T^2+L_3\times T^3+L_4\times T^4+L_5\times T^5}{{10}^8}.\end{array}$$

Using the same method to calculate the periodic coefficient table of the earth’s ecliptic latitude in the heliocentric coordinate system and the periodic coefficient table of the distance between the earth and the sun, the earth’s ecliptic latitude in the heliocentric coordinate system ($B$) and the distance between the sun and the earth ($R$) are obtained, and the results are corrected to the heliocentric ecliptic longitude $L$ and the heliocentric ecliptic latitude $B$ in the standard FK5 system.

The calculation method of ecliptic longitude in geocentric coordinate system ($\Theta$) is as follows. $$\begin{array}{}\text{(4)}& \Theta =L+{180}^{°}.\end{array}$$

The calculation method of ecliptic latitude in the geocentric coordinate system ($\beta$) is as follows. $$\begin{array}{}\text{(5)}& \beta =-B.\end{array}$$

Nutation correction is first used to calculate the apparent longitude and latitude. Nutation can be divided into nutation in longitude ($\Delta \psi$) and nutation in obliquity ($\Delta \epsilon$). Nutation can be described as the sum of some periodic terms. Each periodic term includes calculating the sine coefficient of the nutation in longitude, the cosine coefficient of nutation in obliquity, and the linear combination coefficient of the five basic angular distances ($M$, $M^{\prime }$, $D$, $F$, and $\Omega$) of the amplitude ($\theta$).

The formulas for calculating the five basic angular distances are as follows.

Mean elongation of the Moon from the Sun $$\begin{array}{}\text{(6)}& D=\frac{297.85036+455267.111480\times T-0.0019142\times T^2+T^3}{189474}.\end{array}$$

Mean anomaly of the Sun (Earth) $$\begin{array}{}\text{(7)}& M=\frac{357.52772+35999.050340\times T-0.0001603\times T^2-T^3}{300000}.\end{array}$$

Mean anomaly of the Moon $$\begin{array}{}\text{(8)}& M^{\prime }=\frac{134.96298+477198.867398\times T+0.0086972\times T^2+T^3}{56250}.\end{array}$$

Moon’s argument of latitude $$\begin{array}{}\text{(9)}& F=\frac{93.27191+483202.017538\times T-0.0036825\times T^2+T^3}{327270}.\end{array}$$

Longitude of the ascending node of the Moon’s mean orbit on the ecliptic, measured from the mean equinox of the date $$\begin{array}{}\text{(10)}& \Omega =\frac{125.04452-1934.136261\times T+0.0020708\times T^2+T^3}{450000}.\end{array}$$

The method of calculating the angle $\theta$ is to combine the values calculated from Formulas (6) to (10) with the corresponding five basic angle distance coefficients. According to the amplitude $\theta$, the formula for calculating the periodic terms of the nutation in longitude is as follows. $$\begin{array}{}\text{(11)}& S=S_1+S_2\times T\times \sin \theta .\end{array}$$

The nutation in longitude can be obtained by accumulating the value $S$ of each term, and the unit is ${0.0001}^{\prime }$.

The method of calculating the nutation in obliquity according to the amplitude $\theta$ is similar to that of the nutation in longitude. The formula of periodic term value is as follows. $$\begin{array}{}\text{(12)}& C=C_1+C_2\times T\times \cos \theta .\end{array}$$

In addition to nutation correction, aberration correction is also needed, and the formula is as follows. $$\begin{array}{}\text{(13)}& \text{AC}=\frac{-20^{\prime }4898}{R}.\end{array}$$

The apparent longitude $\lambda$ of the sun in the ecliptic coordinate system is as follows. $$\begin{array}{}\text{(14)}& \lambda =\Theta +\Delta \Psi +\text{AC}.\end{array}$$

The formula for calculating the apparent right ascension and apparent declination of the sun is as follows. $$\begin{array}{}\text{(15)}& \begin{array}{c}\tan \omega =\frac{\sin \lambda \times \cos \epsilon -\tan \beta \times \sin \epsilon }{\cos \lambda },\hfill \\ \sin \delta =\sin \beta \times \cos \epsilon +\cos \beta \times \sin \epsilon \times \sin \lambda .\hfill \end{array}\end{array}$$

In Formula (15), $\epsilon$ is the true obliquity of the ecliptic, which is the sum of the mean obliquity of the ecliptic ${\epsilon }_0$ and the nutation in obliquity $\Delta \epsilon$. The mean obliquity of the ecliptic ${\epsilon }_0$ is given by the following formula. $$\begin{array}{}\text{(16)}& {\epsilon }_0=23.43929111-0.013004167\times T-\frac{0.00059}{3600}\times T^2+\frac{0.00181}{3600}\times T^3.\end{array}$$

According to the application scenario of the method in this paper, the apparent right ascension of the sun is further corrected in order to improve the accuracy, and the correction method is as follows. $$\begin{array}{}\text{(17)}& \tan \alpha =\frac{\tan \omega -0.005}{1+0.005\times \tan \omega }.\end{array}$$

The calculated $\alpha$ and $\Theta$ should be kept in the same quadrant. According to Formulas (1)–(17), given any Julian day $J_D$, the apparent right ascension and apparent declination of the sun in the geocentric inertial coordinate system can be calculated. Then, the sun vector in the geocentric inertial coordinate system is obtained. $$\begin{array}{}\text{(18)}& S_1=\begin{array}{c}\cos \delta \times \cos \alpha \\ \cos \delta \times \sin \alpha \\ \hspace{1em}\sin \delta \end{array}.\end{array}$$

2.2. Instantaneous Root Prediction of Satellite Orbit and Coordinate System Transformation

At present, the prediction methods of instantaneous root of earth satellite orbit are mainly divided into numerical method and analytical method. These methods need too much resources for a single remote sensing instrument on the satellite. Therefore, on the premise of considering both accuracy and simplification, this paper proposes a simple orbital instantaneous root prediction method.

In the case of two-body problem, the satellite in orbit is only affected by the gravity of the central celestial body. In the orbital instantaneous root, except that the true anomaly $f$ changes with time, the other five parameters are semiaxis $a$, eccentricity $e$, inclination $I$, right ascension of ascending node $\Omega$, and the argument of perigee $w$ are constant. However, the satellite in the orbit is subject to various perturbations, and the instantaneous root of the orbit also changes. Two auxiliary quantities, the eccentric anomaly $E$ and the mean anomaly $M$, are introduced. The relationship between the eccentric anomaly $E$ and the true anomaly $f$ is shown in Figure 3.

[figure omitted; refer to PDF]

It can be seen from Figure 3 that the relationship between the eccentric anomaly $E$ and the true anomaly $f$ is as follows. $$\begin{array}{}\text{(19)}& \tan \frac{f}{2}=\sqrt{\frac{1+e}{1-e}}\times \tan \frac{E}{2}.\end{array}$$

The relationship between the mean anomaly $M$ and the eccentric anomaly $E$ is as follows. $$\begin{array}{}\text{(20)}& M=E-e\times \sin E.\end{array}$$

The relationship between the true anomaly $f$ and the mean anomaly $M$ is as follows. $$\begin{array}{}\text{(21)}& f=M+2\text{e}-\frac{e^3}{4}\text{sinM}+\frac{5}{4}{e}^2\sin 2M+\frac{13}{12}{e}^3\sin 3M+•••.\end{array}$$

Among the perturbations of LEO satellites, the earth perturbation $J_2$ term has a significant effect. From the Lagrange type perturbation equation of motion, the perturbation equation of motion of the satellite only with $J_2$ perturbation is obtained. $$\begin{array}{}\text{(22)}& \begin{array}{c}\frac{da}{dt}=\frac{de}{dt}=\frac{di}{dt}=0,\hfill \\ \frac{d\Omega }{dt}=-\frac{3}{2}{nJ}_2{\frac{R_e}{a}}^2\frac{\cos i}{{1-{\text{e}}^2}^2},\hfill \\ \frac{dw}{dt}=\frac{3}{2}{nJ}_2{\frac{R_e}{a}}^2\frac{2-2.5{\sin }^{2}i}{{1-{\text{e}}^2}^2},\hfill \\ \frac{dM}{dt}=n+\frac{3}{2}{nJ}_2{\frac{R_e}{a}}^2\frac{1-1.5{\sin }^{2}i}{{1-{\text{e}}^2}^{1.5}}.\hfill \end{array}\end{array}$$

In Formula (22), $J_2$ is the band harmonic coefficient, which is a constant. $R_e$ is the mean radius of the earth’s equator, which is a constant. In Formula (22), $n$ is the average angular velocity of the satellite. The calculation method of $n$ is as follows. $$\begin{array}{}\text{(23)}& n=\sqrt{\frac{{GM}_e}{a^3}.}\end{array}$$

In Formula (23), $G$ is the gravitational constant and $M_e$ is the mass of the earth.

Formulas (19)–(23) can be used to predict the instantaneous root of the orbit at any time according to the instantaneous root of the orbit at the starting time. Then, the transformation matrix from the geocentric inertial coordinate system to the orbital coordinate system at the prediction time is obtained as follows. $$\begin{array}{}\text{(24)}& T_1=\begin{array}{ccc}-\sin u\cos \Omega -\cos u\cos i\sin \Omega \hfill & -\sin u\sin \Omega +\cos u\cos i\cos \Omega \hfill & \cos u\sin i\hfill \\ -\sin i\sin \Omega \hfill & \sin i\cos \Omega \hfill & -\cos i\hfill \\ -\cos u\cos \Omega +\sin u\cos i\sin \Omega \hfill & -\cos u\sin \Omega -\sin u\cos i\cos \Omega \hfill & -\sin u\sin i\hfill \end{array}.\end{array}$$

In Formula (24), $u=f+w$ is the angle of orbit.

2.3. Calculation of Satellite Correction Matrix

The $X$ axis rotates the roll angle $\phi$, the $Y$ axis rotates the pitch angle $\theta$_, and the $Z$ axis rotates the yaw angle $\psi$ in the orbital coordinate system at time $t$, and the satellite attitude correction matrix can be obtained as follows. $$\begin{array}{}\text{(25)}& T_2=\begin{array}{ccc}\cos \psi \cos \theta -\sin \psi \sin \theta \sin \phi \hfill & \sin \psi \sin \theta +\cos \psi \sin \theta \sin \phi \hfill & -\sin \theta \cos \phi \hfill \\ -\sin \psi \cos \phi \hfill & \cos \psi \cos \phi \hfill & \sin \phi \hfill \\ \cos \psi \sin \theta +\sin \psi \cos \theta \sin \phi \hfill & \sin \psi \sin \theta -\cos \psi \cos \theta \sin \phi \hfill & \cos \theta \cos \phi \hfill \end{array}.\end{array}$$

When the satellite is in orbit, it faces the problems of space environment interference, flexible accessory jitter, and slow time-varying moment of inertia, so the attitude angle of the satellite is unpredictable. The existing methods of satellite attitude determination can accurately estimate the satellite direction relative to the reference system with the help of the measurement information of the attitude sensor. The methods of satellite attitude control, such as active three-axis stability, gravity gradient stability, and spin stability, can ensure high control accuracy and keep the satellite attitude in a stable state. For the application scenario of the method proposed in this paper, the attitude stability mode of the satellite is three-axis stability, and the three-axis pointing accuracy is not greater than 0.15° (3σ), inertial attitude three-axis measurement accuracy is not greater than 0.006° (3σ), and three-axis attitude stability is not greater than 0.006°/s (3σ). According to Formula (25), the error caused by attitude stability is the order of magnitude of 10^-4, so it can be considered that the satellite attitude will not change in the predicted time (unless the change command is injected). Therefore, the satellite attitude correction matrix at the current time is used as the satellite attitude correction matrix at the prediction time.

2.4. Angle Prediction of the Sun in the Coordinate System of the Instrument Body

From Sections 2.1 to 2.3, the sun vector $S_1$ in the geocentric inertial coordinate system, the transformation matrix $T_1$ from Earth inertial coordinate system to the satellite orbit coordinate system, satellite attitude correction matrix $T_2$, and the transformation matrix $T_3$ between the instrument body coordinate system and satellite body coordinate system obtained from the calibration of launch site at the forecast time are obtained. The solar vector $S_2$ in the coordinate system of the instrument body can be obtained by matrix calculation. $$\begin{array}{}\text{(26)}& S_2=T_3{T}_2{T}_1{S}_1.\end{array}$$

3. Results and Discussion

3.1. Results

In order to verify the accuracy of this method, the solar incident angles calculated by this method are compared with those calculated by STK (Satellite Tool Kit). The error of STK in short-time prediction is the order of magnitude of 10^-5, and its prediction error will increase with the increase of prediction time. However, as a mature orbit prediction tool, STK has been widely used, and the error of STK simulation is small, which has little impact on the error between the results of this method and STK simulation results. Therefore, STK simulation results are still used to verify the correctness of the method proposed in this paper.

The applicable scenario of the method proposed in this paper is the satellite carried by the remote sensing instrument under development, so the selected orbit is the orbit of the satellite. The sun synchronous orbit with 836 km orbit height and 9 : 30 : 00 descent intersection time is simulated. The starting time of solar angle forecast is summer solstice, that is, UTCG 2020/06/22 00 : 00 : 00, and the forecast duration is 125 days. Since the period of solar calibration is not fixed, when verifying the results of this method, try to increase the duration of prediction. Therefore, the prediction time is selected as 125 days, which includes the accuracy of different times of short period and long period, so as to verify the error of the method proposed in this paper in different durations.

Without considering the satellite attitude, the method described in Section 2 is used to calculate the $X$, $Y$, and $Z$ components of the sun vector in the satellite body coordinate system. The differences between the results and the STK simulation results are shown in Figures 4–6, respectively.

[figure omitted; refer to PDF]

As shown in Figures 4–6, the errors of $X$, $Y$, and $Z$ components of the solar vector calculated by this method are all within ±0.003°. It can meet the requirements of solar angle for on-orbit sun calibration. With the increase of prediction time, the error of solar angle predicted by the method proposed in this paper will become larger. Therefore, in order to ensure the accuracy of the prediction, it is suggested that the prediction duration should not be too long.

4. Discussion

The error of this method in calculating the solar right ascension and declination and variation of the satellite coordinate system is very small, and the main error comes from the prediction of the orbital instantaneous root. In this paper, only $J_2$ perturbation is considered, so the semiaxis, eccentricity, and inclination do not change with time. When a real satellite moves in orbit, it is affected not only by terrestrial gravity, but also by the earth’s nonspherical gravity, solar gravity, lunar gravity, solar light pressure, atmospheric resistance, and other perturbations, so the orbital instantaneous root changes with time. Taking STK simulation data as the true value, the error of right ascension of ascending node predicted by this method is shown in Figure 7, the error of argument of perigee is shown in Figure 8, and the error of true anomaly is shown in Figure 9. In order to better display the error law of true anomaly, Figure 9 shows the error in seconds within the time range of one day. It can be seen from Figure 7 that the error of the right ascension of ascending node increases with the increase of forecast time, and it can be seen from Figure 8 that the error of the argument of perigee increases with the increase of forecast time. As shown in Figure 9, the error of true anomaly shows an oscillatory change within -0.1°~0.2°, which does not become larger and larger with the increase of forecast time, and the error is always concentrated in a range. The law of true anomaly error at other times is the same.

[figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF]

5. Conclusions

In order to grasp the timing of solar calibration in advance, a simple method of high-precision solar angle prediction for space remote sensing instruments is proposed in this paper. This method considers both accuracy and simplification. This method is used to predict the solar angle of sun synchronous orbit in the satellite body coordinate system and compared with the STK simulation results. The results show that the prediction error of the solar vector is less than ±0.003° in 125 days. It can meet the requirements of solar angle accuracy for on-orbit sun calibration. At the same time, it provides a reference for other space remote sensing instruments on-orbit solar angle prediction.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (No. 62005268).

References

[1] G. Chander, T. J. Hewison, N. P. Fox, X. Wu, X. Xiong, W. J. Blackwell, "Overview of intercalibration of satellite instruments," IEEE Transactions on Geoscience & Remote Sensing, vol. 51 no. 3, pp. 1056-1080, DOI: 10.1109/TGRS.2012.2228654, 2013.

[2] C. G. CUI, L. LI, Y. F. LI, Y. X. GUO, Y. Q. LI, "External calibration methods for geosynchronous ultraviolet-visible hyperspectral instrument," Optics and Precision Engineering, vol. 29 no. 3, pp. 484-492, DOI: 10.37188/OPE.2020.0522, 2021.

[3] R. H. Evans, H. R. Gordon, "Coastal zone color scanner ‘system calibration’: a retrospective examination," Journal of Geophysical Research Oceans, vol. 99 no. C4, pp. 7293-7307, DOI: 10.1029/93JC02151, 1994.

[4] K. J. Thome, "Absolute radiometric calibration of Landsat 7 ETM+ using the reflectance-based method," Remote Sensing of Environment, vol. 78 no. 1-2, pp. 27-38, DOI: 10.1016/S0034-4257(01)00247-4, 2001.

[5] B. Schulz, S. Huth, R. J. Laureijs, J. A. Acosta-Pulido, M. Braun, H. O. Castañeda, M. Cohen, L. Cornwall, C. Gabriel, P. Hammersley, I. Heinrichsen, U. Klaas, D. Lemke, T. Müller, D. Osip, P. Román-Fernández, C. Telesco, "ISOPHOT-photometric calibration of point sources," Astronomy and Astrophysics, vol. 381 no. 3, pp. 1110-1130, DOI: 10.1051/0004-6361:20011538, 2002.

[6] M. A. Folkman, J. Pearlman, L. B. Liao, P. J. Jarecke, "EO1/hyperion hyperspectral image design, development, characterization, and calibration," Hyperspectral Remote Sensing of the Land & Atmosphere, vol. 4151, pp. 40-51, 2001.

[7] J. S. Pearlman, P. S. Barry, C. C. Segal, J. Shepanski, D. Beiso, S. L. Carman, "Hyperion, a space-based imaging spectrometer," IEEE Transaction on Geoscience & Remote Sensing, vol. 41 no. 6, pp. 1160-1173, DOI: 10.1109/TGRS.2003.815018, 2003.

[8] D. Roy, J. Ju, P. Lewis, C. Schaaf, F. Gao, M. Hansen, E. Lindquist, "Multi-temporal MODIS-Landsat data fusion for relative radiometric normalization, gap filling, and prediction of Landsat data," Remote Sensing of Environment, vol. 112 no. 6, pp. 3112-3130, DOI: 10.1016/j.rse.2008.03.009, 2008.

[9] Q. Zhang, X. Ye, W. Fang, "Calculation of sun’s location on the sun synchronous orbit satellite and its precision analysis," Journal of the Graduate School of the Chinese Academy of Sciences, vol. 28 no. 3, pp. 310-314, 2011.

[10] D. J. Yang, W. Fang, X. Ye, B. Q. Song, "Program-controlled sun-tracking precision of spaceborne solar irradiance monitor," Optics and Precision Engineering, vol. 23 no. 7, pp. 1813-1821, 2015.

[11] Z. F. Li, S. R. Wang, Y. Huang, G. Lin, Y. Shao, M. Yu, "High-Precision and Short-Time Solar Forecasting by Spaceborne Solar Irradiance Spectrometer," Acta Optica Sinica, vol. 39 no. 7, article 0712007, 2019.

[12] P. Bretagnon, G. Francou, "Planetary theories in rectangular and spherical variables: VSOP87 solutions," Astronomy & Astrophysics, vol. 202 no. 1/2, pp. 309-315, 1998.

[13] J. J. Michalsky, "The astronomical almanac 's algorithm for approximate solar position (1950–2050)," Solar Energy, vol. 40 no. 3, pp. 227-235, DOI: 10.1016/0038-092X(88)90045-X, 1988.

[14] J. Meeus, Astronomical algorithms, 1991.