On the Insolation at the Top of the Earth’s Atmosphere and Its Variation under Smooth Changes of Astronomical Elements during the Past Four Centuries

Gerhard Kramm¹, Nicole Mölders²

¹Engineering Meteorology Consulting, Fairbanks, AK, USA.
²Department of Atmospheric Sciences and Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK, USA.

DOI: 10.4236/ns.2025.176009                

Abstract

Our objective is to (a) investigate the variation of insolation at the top of the atmosphere (TOA) under secular changes of astronomical elements (e.g., Earth’s heliocentric distance, Sun’s apparent geocentric declination) during the past four centuries, and (b) assess the accuracy of historical work on insolation in this historic context. Therefore, we predict the insolation at the TOA for timescales from the diurnal course to the annual course and arbitrary periods of days using a 10-minute timestep and an increment of 5˚. Our predications are based on the ICRS geocentric-rectangular coordinates and their rates provided by the JPL planetary and lunar ephemeris DE440 and the adjustment of them to the equator and the equinox-of-date using the subroutines for aberration, frame-bias, precession, and nutation of the Naval Observatory Vector Astrometry Software, Version F3.1. At 1AU, we apply the solar constant of the respective year, taken from the reconstructed total solar irradiance (TSI) based on the “Community-Consensus TSI Composite” and SATIRE-T model. We quantify the variation as the differences in the daily mean solar irradiance over the annual course compared to 2010. Although the reconstructed TSI varies only by about 2.2 W∙m⁻², seasonal daily mean solar radiation differs notably from that of 2010, particularly for the polar regions of the southern hemisphere (SH) and the northern hemisphere (NH). There, the differences have mainly opposite signs, resulting in a butterfly-like distribution across latitudes and seasons. The 1910 daily mean solar irradiance differs the largest from that of 2010 (SH: −7.3 W∙m⁻² to 7.4 W∙m⁻²; NH: −7.4 W∙m⁻² to 7.3 W∙m⁻²) followed by 1810 (SH: −5.0 W∙m⁻² to 5.5 W∙m⁻²; NH: −5.8 W∙m⁻² to 5.2 W∙m⁻²). For 1910, 1810, 1950, and 1710, the daily mean solar irradiance is most sensitive to smooth changes in the astronomical elements from the polar circles poleward and least sensitive around the equator. The smallest differences compared to 2010 occur for 1750 (SH: −1.9 W∙m⁻² to 1.1 W∙m⁻²; NH: −2.0 W∙m⁻² to 1.2 W∙m⁻²), when the largest positive differences occur around the equator. Our results reveal a noteworthy nonlinear relationship between the annual course of daily mean solar irradiance and the combined effects of Earth’s heliocentric distance and the Sun’s apparent geocentric declination. Over the annual course, the Sun’s apparent geocentric declination dominates the deviations in daily mean solar irradiance. Given the magnitude of these differences in solar irradiance, the smooth changes of the astronomical elements must be considered on the multi-decadal to multi-centennial scales.

Keywords

Solar Irradiance, Top of the Atmosphere, Heliocentric Distance, Geocentric Declination, Hour Angle, JPL Planetary and Lunar Ephemeris DE440, Naval Observatory Vector Astrometry Software (NOVAS) F3.1, International Celestial Reference System (ICRS), Geocentric Celestial Reference System (GCRS), Frame Bias, Precession of the Equator, Precession of the Ecliptic, Nutation, Earth Rotation Angle (ERA), Celestial Intermediate Origin (CIO), Celestial Intermediate Pole (CIP), Terrestrial Intermediate Origin (TIO), Total Solar Irradiance (TSI), Community-Consensus TSI Composite, SATIRE-T Model

1. Introduction

The Earth-atmosphere system (EAS) is continuously rotating in the Sun’s radiation field from which the entire EAS receives solar energy during each daily rotation. The input of solar energy into the EAS at a given location on the surface spanned by the top of the atmosphere (TOA) depends on the Earth’ heliocentric distance, $r$ , the geocentric declination, ${\delta }_S$ , and the local hour angle, $h$ , of the Sun (e.g., [1-12]). The TOA may be interpreted as the height of the intervening atmospheric layer. Above this height, neither solar radiation nor infrared radiation is notably affected by gaseous or particulate matter. Generally, $r$ , ${\delta }_S$ and $h$ depend on time, but on different time scales. While $h$ serves to describe the daily course of insolation at a given location between sunrise and sunset, $r$ and ${\delta }_S$ are necessary to capture, for instance, the seasonal variation and annual course of insolation. Since $r$ and ${\delta }_S$ are affected by the variation of the Earth’s orbit around the Sun due to smooth changes in the astronomical conditions caused by precession of both the equator and the ecliptic, nutation, changes in the obliquity (and long-term changes in eccentricity), it is indispensable to predict the astronomical elements $r$ and ${\delta }_S$ using trustworthy astronomical calculation methods.

Since the findings of Meech [1] and Wiener [2, 3] that the global average of the daily mean solar irradiance at the TOA is equal to one quarter of the solar constant, later confirmed by Milankovitch [4, 5], List [6], Fortak [13], Raschke et al. [14], Liou [9], Fu [10], Berger and Yin [15] and many others, has recently been disputed by experts of the German organization European Institute for Climate and Energy (EIKE), the goals of our study are: 1) to assess the prediction of insolation at the TOA for different time scales ranging from the diurnal course to the annual course and arbitrary periods of days in between, 2) to review the historical work on solar irradiance performed by Meech and Wiener and assess the accuracy of their findings using our results obtained for the year 1853 and the tropical years 1874/1875 and 2009/2010, and 3) to quantify the impact of smooth changes in the astronomical elements $r$ and ${\delta }_S$ on the daily mean irradiance over the annual course during the past four centuries by predicting the daily mean solar irradiance for the years 1610, 1650, 1710, 1750, 1810, 1850, 1910 and 1950 and comparing the results with those obtained for 2010. We also confirm the results of Kopp [16] for 2023.

Our predictions are based on the ICRS geocentric rectangular coordinates and their rates provided by the planetary and lunar ephemeris DE440 of the Jet propulsion Laboratory (JPL), California Institute of Technology (Park et al. [17]), hereafter called the JPL ephemeris DE440, and the adjustment of these data to the equator and the equinox of date (see Section 3) using the subroutines for aberration, frame bias, precession and nutation of the Naval Observatory Vector Astrometry Software (NOVAS), Version F3.1 (Kaplan et al. [18, 19]), hereafter called the NOVAS F3.1 subroutines. Our predictions were performed for 90˚S to 90˚S N latitudes using an increment of 5˚, where a 10-minute time step was used. For 2009 and 2010, the value of the total solar irradiance (TSI) at 1AU, the so-called solar constant, of $S=1361\text{W}\cdot {\text{m}}^{-2}$ was chosen, in accord with Laue and Drummond [20], Raschke et al. [14], Kopp and Lean [21], and Kopp et al. [22]). For all other years considered in our study, we took the data from the reconstruction of the TSI for the past four centuries that are based on the “Community-Consensus TSI Composite” [23] and SATIRE-T model [24] (with modifications to fix spurious values prior to 1650).

Our paper is organized as follows: In Section 2, we give an overview regarding the historical work on solar irradiation. In Section 3, we discuss the astronomical aspects, including the state-of-the-art coordinate frame. The reconstruction of TSI for the past four centuries is briefly discussed in Section 4. In Sections 5 to 8, we discuss the solar irradiance from the diurnal course to the annual course and an arbitrary period of days in between, as well as the interannual, and seasonal variations of the solar irradiance during the past four centuries, all based on our simulations. We conclude on the impact of astronomical conditions on daily mean solar irradiance in the annual course in Section 9.

2. historical overview

The solar irradiance reaching the TOA at a certain location is given by (e.g., [1, 3-5, 8, 9, 25])

$F_{s\downarrow }\left(\varphi ,\lambda \right)=F\cos {\Theta }_0$ ,(1)

where $\varphi$ is the latitude, $\lambda$ is the longitude, and $F$ is the TSI at the Earth’s actual heliocentric distance. The cosine of ${\Theta }_0$ occurs because only the component of the solar radiation flux density, $F$ , that is perpendicular to the (infinitesimal) surface element is considered. This component is obtained using the scalar product between the vector $-F$ and the normal vector $n$ of this surface element (positively counted if $n$ it is directed outwards), i.e., $-F\cdot n=-\left|F\right|\left|n\right|\cos \left(\pi -{\Theta }_0\right)=\left|F\right|\cos {\Theta }_0$ , where $F=\left|F\right|$ denotes the magnitude of $F$ , and $\left|n\right|=1$ is the magnitude of $n$ . If the surface element is horizontally aligned as always considered here, ${\Theta }_0$ is the local zenith distance of the Sun’s center.

The TSI at the Earth’s heliocentric distance can be expressed by

$F={\left(\frac{r_S}{r}\right)}^2{E}_S$ .(2)

Here, $r_S\cong 696300\text{km}$ is the radius [26] and $E_S$ the emittance of the Sun. This formula implies that the radiant power ($=4\pi r_S^2{E}_S$ ) of the Sun is constant when the solar radiation is propagating through the space because of energy conservation principles in the absence of an intervening medium (e.g., [5, 9, 11, 27, 28]). In accord with Equation (2), we obtain for the Earth’s mean heliocentric distance, $r_0$ , for which the so-called solar constant, $S$ , is defined (e.g., [29, 30])

$S={\left(\frac{r_S}{r_0}\right)}^2{E}_S$ .(3)

Combining Equations (2) and (3) yields

$F={\left(\frac{r_0}{r}\right)}^{2}S=\frac{S}{{\rho }^2}$ .(4)

where $\rho =r/r_0$ is the relative heliocentric distance. The quantity ${\left(r_0/r\right)}^2$ may be called the orbital effect. It does not vary more than 3.5 percent (e.g., [9, 25]). Note that $r_0$ is often replaced by the semi-major axis $a$ (see also Equation (51)). We choose, however, $r_0=1\text{AU}$ .

The mean angular velocity of the Earth’s rotation corresponds to ${\Omega }_E=2\pi /d_s\cong 7.2921\times {10}^{-5}{\text{s}}^{-1}$ , where $d_s\cong 86164.09\text{s}$ is the mean sidereal day. Because this angular velocity of rotation is much larger than the angular velocity of the Earth’s revolution around the Sun ranging from about $2.060\times {10}^{-7}{\text{s}}^{-1}$ at perihelion and about $1.926\times {10}^{-7}{\text{s}}^{-1}$ at aphelion, the effect of the Earth’s orbital angular velocity is negligible¹. Therefore, in the case of the Earth, the function $\cos {\Theta }_0$ can be determined using rules of spherical trigonometry leading to (e.g., [1, 3-6, 8-11, 25, 29])

$\cos {\Theta }_0=\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h$ .(5)

Here, ${\delta }_S$ and $h$ are again the geocentric declination and the local hour angle of the Sun, where $h=0$ defines the local solar noon that meets the criterion

$\cos {\Theta }_0=\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S=\cos \left(\varphi -{\delta }_S\right)$ ,(6)

and, hence, ${\Theta }_0=\varphi -{\delta }_S$ . Note that $h$ counted from the meridian at local solar noon, addresses the rotation of the Earth around its polar axis, i.e., it serves to describe the diurnal variation of the insolation at the TOA at a specific location between sunrise, $h=-H$ , and sunset, $h=H$ , where $H$ is the so-called half-day, and $2H$ is the diurnal arc.

The Sun’s apparent geocentric declination is related to the apparent ecliptic longitude of the Sun, ${\lambda }_S$ , counted from the dynamic vernal equinox of the northern hemisphere (NH), hereafter called the March equinox, by (see also Equation (32)) [3, 4, 5, 8, 35]

$\sin {\delta }_S=\sin {\lambda }_S\sin \epsilon$ ,(7)

where $\epsilon \cong 23˚26\text{'}18\text{'}\text{'}.42$ is the true obliquity of the ecliptic in 2025. (This means that the solar irradiance at the TOA is slightly affected by changes in the obliquity.) The apparent ecliptic longitude of the Sun differs from the true longitude of the Earth, $\lambda =\upsilon +\varpi$ , by ${\lambda }_S=\lambda -180˚$ . Here, $\upsilon$ is the true anomaly counted counterclockwise from the perihelion, $\varpi =\omega +{\Omega }_a$ is the longitude of the perihelion relative to the March equinox (see Figure 1), $\omega$ is the argument of the perihelion, and ${\Omega }_a$ is the longitude of the ascending node. Note that ${\delta }_S$ ranges from ${\delta }_S=-\epsilon$ (Tropic of Capricorn) at ${\lambda }_S=270˚$ to ${\delta }_S=\epsilon$ (Tropic of Cancer) at ${\lambda }_S=90˚$ (see Figure 2(b)). Furthermore, for ${\lambda }_S=0˚$ and ${\lambda }_S=180˚$ , the condition of the equinoxes, formula (7) provides ${\delta }_S=0$ [36, 37].

According to Equations (1), (4) and (5), the TSI reaching the TOA at a certain location is given by (e.g., [1, 3-5, 8-10, 25])

$F_{S\downarrow }\left(\varphi ,{\delta }_S,h,\rho \right)=\frac{S}{{\rho }^2}\left(\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h\right)$ .(8)

Because $F_{S\downarrow }\left(\varphi ,{\delta }_S,h,\rho \right)$ is expressed in energy (J) per unit area (m²) and unit time (s), we may write

$F_{S\downarrow }\left(\varphi ,{\delta }_S,h,\rho \right)=\frac{\text{d}{W}_{S\downarrow }\left(\varphi ,{\delta }_S,h,\rho \right)}{\text{d}t}$ ,(9)

where $\text{d}{W}_{S\downarrow }\left(\varphi ,{\delta }_S,h,\rho \right)$ is the amount of solar energy that is flowing through the surface element during the time element dt (e.g., [1-6, 9]). Consequently,

$\frac{\text{d}{W}_{S\downarrow }\left(\varphi ,{\delta }_S,h,\rho \right)}{\text{d}t}=\frac{S}{{\rho }^2}\left(\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h\right)$ .(10)

We discuss the integration of this equation over different time scales in Sections 5 - 8. For readability, we omit the functional dependences expressed by $\left(\varphi ,{\delta }_S,h,\rho \right)$ hereafter.

Note that instead of $W_{S\downarrow }$ , the normalized term $w_{S\downarrow }=W_{S\downarrow }/S$ may be used so that $\text{d}{w}_{S\downarrow }/\text{d}t=\left(\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h\right)/{\rho }^2$ .

Inserting Equation (7) into Equations (5) and (10) yields [8, 38, 39]

$\cos {\Theta }_0=\sin \varphi \sin \epsilon \sin {\lambda }_S+\cos \varphi \sqrt{1-{\sin }^2\epsilon {\sin }^2{\lambda }_S}\cos h$ (11)

and

$\frac{\text{d}{W}_{S\downarrow }}{\text{d}t}=\frac{S}{{\rho }^2}\left(\sin \varphi \sin \epsilon \sin {\lambda }_S+\cos \varphi \sqrt{1-{\sin }^2\epsilon {\sin }^2{\lambda }_S}\cos h\right)$ .(12)

Figure 1. Elements of the Earth’s orbit (with reference to Berger [40]). The orbit of the Earth, E, around the Sun, S, is represented by the ellipse PAE, P being the perihelion and A the aphelion, $a=\overline{OA}$ being the semi-major axis and $b=\overline{OB}$ the semi-minor axis. Furthermore, $\gamma$ is the vernal point, WS and SS are the winter and summer solstices of the NH, respectively. They mirror their present-day locations. The vector $n$ is perpendicular to the ecliptic, and the obliquity, $\epsilon$ , is the inclination of the equator upon the ecliptic; i.e., $\epsilon$ is equal to the angle between the Earth’s axis of rotation and $n$ . The quantity $\upsilon$ is the true anomaly counted counterclockwise from the perihelion. Furthermore, the quantity $\varpi$ is the longitude of the perihelion relative to the dynamic vernal equinox (VE) and is equal to $\xi +\psi$ . The annual general precession in longitude, $\psi$ , describes the absolute motion of $\gamma$ along the Earth’s orbit relative to the fixed stars. The longitude of the perihelion, $\xi$ , is measured from the reference vernal equinox of the standard epoch (e.g., J2000.0) and describing the absolute motion of the perihelion relative to the fixed stars.

Alternatively, we may use $\cos {\delta }_S=\cos \left(\arcsin \left(\sin {\lambda }_S\sin \epsilon \right)\right)$ to obtain [35]

$\frac{\text{d}{W}_{S\downarrow }}{\text{d}t}=\frac{S}{{\rho }^2}\left(\sin \varphi \sin \epsilon \sin {\lambda }_S+\cos \varphi \cos \left(\arcsin \left(\sin \epsilon \sin {\lambda }_S\right)\right)\cos h\right)$ .(13)

Equations (12) and (13) describe the long-term variation of the solar irradiation at a specific latitude as a function of the solar constant, the eccentricity, the obliquity, the Sun’s apparent ecliptic longitude, and the local hour angle. On the scales of multiple thousands of years (kyr), we have to pay attention to Milankovitch’s [4, 5] astronomical theory of climatic variations that ranks as the most important achievement in the theory of climate in the 20^th century [11, 41]. In accord with Berger [40], such long-term changes are denoted as climatic variations. Milankovitch’s astronomical theory regards the change of the eccentricity and the obliquity, and to precession and nutation phenomena due to the perturbations that Sun, Moon, and the principal planets of our solar system exert on the Earth’s orbit (e.g., [9, 35, 38, 40-43]) ideally characterized by Equations (47) to (54). His theory plays a substantial role in the time series analysis of paleoclimate records (see, e.g., [43, 44]). Because of these astronomical phenomena, the insolation at the TOA will vary during such long-term periods. Recently, Smulsky [35] presented a new astronomical theory of climatic variations that addresses three problems: the evolution of orbital motion, the evolution of the Earth’s rotational motion, and the evolution of the insolation controlled by the evolution of these motions.

Figure 2 illustrates the Earth’s heliocentric distance and orbital velocity, the Sun’s apparent declination, and the total solar irradiance versus the Sun’s apparent ecliptic longitude, ${\lambda }_S$ , computed for the tropical year 2009/2010 beginning with the March equinox, TDB = 2454910.9931 (March 20, 2009, 11:50 UT1). The required data were provided by the JPL planetary and lunar ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. In accord with Laue and Drummond [20], Raschke et al. [14], Kopp and Lean [21], and Kopp et al. [22]), we used a solar constant of $S=1361\text{W}\cdot {\text{m}}^{-2}$ .

According to Park et al. [17], the inertial coordinate frame of DE440 is connected to the International Celestial Reference System (ICRS). The current ICRS realization is achieved by Very Long Baseline Interferometry (VLBI) measurements of the positions of extragalactic radio sources (i.e., quasars) defined in the Third Realization of the International Celestial References Frame (ICRF3; Charlot et al. [45]), which is adopted by the International Astronomical Union (IAU). As pointed out by Park et al. [17], the orbits of the inner planets are tied to ICRF3 via VLBI measurements of the Mars-orbiting spacecraft (Konopliv et al. [46]) with respect to quasars with positions known in the ICRF. The ephemeris DE440 is valid for the years 1550 - 2650.

As shown by Figure 2(b), the beginning of the astronomical seasons of the NH refers to the values of the Sun’s apparent declination at the apparent ecliptic longitudes: ${\lambda }_S=0˚$ (spring), ${\lambda }_S=90˚$ (summer), ${\lambda }_S=180˚$ (fall), and ${\lambda }_S=270˚$ (winter) [36, 37]. In 1827, however, Fourier pointed out:

“Dans cette hypothèse du froid absolu de l’espace, s’il est possible de la concevoir, tous les effets de la chaleur, tels que nous les observons àla surface du globe, seraient dus àla présence du soleil. Les moindres variations de la distance de cet astre à la terre occasioneraient des changements très considérables dans les températures, l’excentricité de l’orbite terrestre donnerait naissance à diverses saisons.”

Translated from French into English:

“In this hypothesis of the absolute cold of space, if it is possible to conceive all the effects of heat that we observe on the surface of the globe, would be due to the presence of the Sun. The slightest variations in the distance of this star from the Earth would cause very considerable changes in temperatures, the eccentricity of the Earth’s orbit would give rise to various seasons.”

Fourier’s statements are misleading. The various astronomic seasons are caused by the Sun’s geocentric declination. It is well known, for instance, that the astronomic seasons of the NH differ from those of the southern hemisphere (SH) even for the same Earth’ heliocentric distance. The eccentricity mainly affects the variation of the TSI (see Equations (2) to (4) and Figure 2(d)) from about 1317 W∙m⁻² at aphelion

(${\lambda }_S\cong 104.2˚$ ), the maximum of the Earth’s heliocentric distance $r_a=a\left(1+e\right)\cong 1.521\times {10}^8\text{km}$ , to about 1408 W∙m⁻² at perihelion (${\lambda }_S\cong 282.9˚$ ), the minimum of the Earth’s heliocentric distance $r_p=a\left(1-e\right)\cong 1.471\times {10}^8\text{km}$ , where $a\cong 1.496\times {10}^8\text{km}$ and $e\cong 0.01671$ are the current semi-major axis

and the eccentricity, respectively. In addition, the eccentricity also affects the Earth’s orbit velocity that varies from about 29.29 km∙s⁻¹ at aphelion to about 30.29 km∙s⁻¹ at perihelion (see Figure 2(c)) and, hence, the lengths of the astronomic seasons [15], i.e., March equinox-June solstice, 92.8 d; June solstice-September equinox, 93.6 d; September equinox-December solstice, 89.8 d; December solstice-March equinox, 89.0 d. This means that the summer half-year (spring plus summer) of the NH, is about 7.6 d longer than that of the SH (e.g., [2, 3, 5, 15]).

The Earth’s heliocentric distance varies because of the orbital motion of the Earth around the Sun where, according to Kepler, the Earth’s orbit is nearly elliptical in shape with the Sun at one focus (see Figure 1). The mean plane spanned by the Earth’s orbit is called the ecliptic or the ecliptic plane, strictly formulated as follows (see Glossary of the IAU Division I Working Group\Nomenclature for Fundamental Astronomy (NFA)): The ecliptic is the plane perpendicular to the mean heliocentric orbital angular momentum vector of the Earth-Moon barycenter in the Barycentric Celestial Reference System (BCRS). The ecliptic is also the apparent path of the Sun around the celestial sphere [47].

As mentioned before, the ecliptic is currently inclined by $\epsilon \cong 23˚26\text{'}19\text{'}\text{'}$ to the celestial equator, that is Earth’s equatorial plane projected on the celestial sphere. The ecliptic intersects the celestial equator at two points: The March equinox, $\Upsilon$ , and the September equinox. The March equinox is at the ascending node of the ecliptic on the equator. It is the direction at which the Sun, in its annual apparent path around the Earth, crosses the equator from south to north. It is also referred to as “the First Point of Aries”. The inclination of the ecliptic plane to the plane of the equator, i.e., $\epsilon$ , is the obliquity of the ecliptic [47, 48].

(a) (b)

(c) (d)

Figure 2. (a) Earth’s heliocentric distance and, (b) Sun’s apparent geocentric declination, (c) Earth’s orbital velocity and (d) the total solar irradiance versus the Sun’s apparent ecliptic longitude computed for the tropical year 2009/2010 beginning with the March equinox, TDB = 2454910.9931 (March 20, 2009, 11:50 UT1). The required data were provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. A solar constant of $S=1361\text{W}\cdot {\text{m}}^{-2}$ was used.

As pointed out by Wiener [2, 3], Lambert [49] was the first who dealt with this topic in expanded form, even by considering the average influence of the atmosphere and soil. Wiener also mentioned that Meech developed the formulas and calculated several tables [1]. Meech determined the annual irradiation for 1853 but avoided considering seasonal variations. The latitudinal and seasonal distribution of the normalized daily mean solar irradiance at the TOA for that year is illustrated in Figure 3. Table 1 lists the results of Lambert and Meech used by Wiener for comparison with his own results. To assess these historical results, we added the results of our calculations for 1853 and the tropical years 1874/75 and 2009/10.

Ignoring the effects of the atmosphere, Wiener [2, 3] computed the irradiation at the Earth’s surface for 17 equidistantly distributed days (with respect to $\Delta {\lambda }_S=22.5˚$ ) for the tropical year between March 20, 1874 and March 21, 1875, subdividing the northern and southern hemispheres at intervals of $\Delta \varphi =10˚$ . Because he did not know the solar constant, he considered the normalized solar irradiance. His results are illustrated in Figure 4(a). This latitudinal and seasonal distribution of the daily mean solar irradiance of the Earth without atmosphere, normalized by the solar constant, is nearly identical with that at the TOA because

Figure 3. Latitudinal-annual distribution of the daily mean irradiance at the TOA for 1853 normalized by the solar constant. TDB = 2397854.5 corresponds to January 1, 1853, 00:00 UT1. The Earth’ heliocentric distance and the Sun’s apparent geocentric declination were calculated using the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. The dashed line indicates the apparent geocentric declination of the Sun, and the dotted lines the Equator, the Tropic of Cancer, the Tropic of Capricorn, the Arctic Circle, and the Antarctic Circle. The daily mean values were computed with the Equation (14) using 144 data per day.

the thickness of the intervening atmospheric layer of about ${\vartheta }_a\cong 150\text{km}$ would not notably change the amount of solar irradiance. At the distance at perihelion of $r=1.47098\times {10}^8\text{km}$ (in 2010), for instance, a change by ${\vartheta }_a\cong 150\text{km}$ would lead to ${\left(1-{\delta }_a/r\right)}^{-2}\cong 1.000002$ . This means that the maximum of solar

irradiance at the subsolar point of about 1408 W∙m⁻² would be increased by 0.002 W∙m⁻². This increase is negligible because it is below the accuracy with which the TSI can be measured. Note that the daily average of an arbitrary quantity $\psi \left(\varphi ,\phi ,t\right)$ at a given location is defined by

$\overline{\psi }\left(\varphi ,\phi \right)=\frac{1}{d}{\displaystyle \underset{0}{\overset{d}{\int }}\psi \left(\varphi ,\phi ,t\right)\text{d}t}$.(14)

Milankovitch [5] weighted Wiener’s results with a solar constant of $S=2\text{cal}\cdot {\text{cm}}^{-2}\cdot {\min }^{-1}$ . Later, List [6] re-scaled Milankovitch’s results by using a solar constant of $S=1.94\text{cal}\cdot {\text{cm}}^{-2}\cdot {\min }^{-1}$ . List’s diagram of the latitudinal-annual distribution of the daily solar irradiance at the TOA has been used widely in the literature (e.g., [7, 27, 50, 51]).

For assessing the accuracy of Wiener’s results, we considered the tropical year as well beginning on March 20, 1874, 00:00 UT1 ($TDB=2405602.5$ ), where a subdivision of both hemispheres by $\Delta \varphi =5˚$ and a time step of 600 s were used. The Earth’s heliocentric distance and the Sun’s apparent geocentric declination, both time-dependent, were determined using again the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. Obviously, Wiener’s results agree well with our results, which is remarkable because Wiener did not have a computer. He derived his results using three different methods: 1) graphical, 2) mechanical quadrature, and 3) solution of elliptical integrals of first, second and third kind solved with the aid of Legendre’s tables. Only the results of the latter are listed in Table 1.

Using the zonal average defined by (e.g., [7, 52-55]),

$\widehat{\psi }\left(\theta \right)=\frac{1}{2\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}\psi \left(\theta ,\phi \right)\text{d}\phi }$,(15)

(a) (b)

Figure 4. Latitudinal-annual distribution of the daily mean irradiance of the Earth in the absence of its atmosphere normalized by the solar constant, (a) according to Wiener [2, 3] and (b) this work, where the Earth’s heliocentric distance and the Sun’s apparent geocentric declination were calculated using the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. The dashed lines show the Sun’s apparent geocentric declination, and the dotted lines the Equator, the Tropic of Cancer, the Tropic of Capricorn, the Arctic Circle, and the Antarctic Circle. The daily mean values were computed with Equation (1) using 144 data per day.

where $\psi \left(\theta ,\phi \right)$ is a field quantity like the daily solar irradiance and $\theta =\pi /2-\varphi$ , the global average of $\langle \psi \rangle$ is given by (e.g., [7, 52, 53, 55-58])

$\langle \psi \rangle =\frac{{\displaystyle \underset{\Omega }{\int }\psi \left(\theta ,\phi \right)\text{d}\Omega }}{{\displaystyle \underset{\Omega }{\int }\text{d}\Omega }}=\frac{1}{4\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\displaystyle \underset{0}{\overset{\pi }{\int }}\psi \left(\theta ,\phi \right)\sin \theta \text{d}\theta \text{d}\phi }}=\frac{1}{2}{\displaystyle \underset{0}{\overset{\pi }{\int }}\widehat{\psi }\left(\theta \right)\sin \theta \text{d}\theta }$.(16)

Here, $\Omega =4\pi$ is the solid angle of a sphere, and $\text{d}\Omega =\sin \theta \text{d}\theta \text{d}\phi$ is the differential solid angle, where $\theta$ and $\phi$ are the zenith and azimuthal angles in a right-handed spherical coordinate frame that serve to characterize the location; $\theta$ ranges from 0 (North Pole) to $\pi$ (South Pole), and $\phi$ ranges from 0 to $2\pi$ , where $\phi =0$ may be considered to be equal to the prime meridian, i.e., the zero-longitude, about 100 m east of the transit circle at the Royal Observatory, Greenwich [47].

Table 1. Zonal averages of the normalized daily mean solar irradiance for different parallels of latitude computed by Lambert [49], Meech [1], Wiener [2, 3]. Lambert’s data were adopted from Wiener [2, 3]. Wiener’s and our results are related to the tropical years starting with the March equinoxes of 1874 and 2009.

Latitude in ˚	Lambert	Meech	Wiener	This work
		1853	1874/1875	1853	1874/1875	2009/2010
90	0.12677	0.12674	0.12672	0.12682	0.12681	0.12672
85		0.12775		0.12786	0.12784	0.12775
80		0.13093	0.13096	0.13105	0.13101	0.13092
75		0.13644		0.13655	0.13649	0.13641
70		0.14464	0.14464	0.14477	0.14469	0.14463
66.554	0.15259
65		0.15704		0.15715	0.15706	0.15702
60		0.17368	0.17368	0.17377	0.17368	0.17365
55		0.19128		0.19138	0.19129	0.19127
50		0.20878	0.20876	0.20886	0.20877	0.20876
45	0.22553	0.22552		0.22562	0.22554	0.22553
40		0.24122	0.24122	0.2413	0.24122	0.24122
35		0.25553		0.25563	0.25555	0.25555
30		0.26834	0.26832	0.2684	0.26833	0.26834
25		0.27936		0.27945	0.27939	0.2794
23.446	0.28241
20		0.28857	0.28858	0.28865	0.28861	0.28862
15		0.29584		0.29591	0.29587	0.29589
10		0.30112	0.30112	0.30115	0.30112	0.30113
5		0.30427		0.3043	0.30428	0.3043
0	0.30532	0.30532	0.30532	0.3053	0.30529	0.3053

The zonal averages listed in Table 1 provide the following global averages of the solar irradiance normalized by the solar constant:

Lambert (1779): 0.24805

Meech (1857): 0.24995

Wiener (1879): 0.24972

This work for 1853: 0.24995

This work for 1874/1875: 0.24996

This work for 2009/2010: 0.24996

These results suggest that the global average of the daily mean solar irradiance at the TOA is equal to one quarter of the solar constant [1, 3-5]. Tests for 2010 using an increment of $\Delta \varphi =2.5˚$ provide 0.25001.

2.1. Spitaler’s Criticism

Spitaler [59] disputed the historical findings especially those of Wiener [2, 3] and stated:

“Eine besondere Bedeutung erlangte eine Abhandlung von Chr. Wiener, als das Problem für die Erforschung der klimatischen Verhältnisse der Eiszeit von Bedeutung wurde, indem er fand, daß zur Zeit des Sommersolstitiums die tägliche Bestrahlung des Pols 4/3 derjenigen ist, welche zur selben Zeit am Äquator herrscht. Dieses Ergebnis ist auch in die Hand- und Lehrbücher der Meteorologie übergegangen und hat besonders bei den Nachforschungen über das Polarklima in der Tertiärzeit großes Interesse gefunden.

Ich veranlaßte daher schon vor längerer Zeit meinen ehemaligen Schüler Fr. Hopfner, das Problem der Bestrahlung einer gründlichen mathematischen Revision zu unterziehen, und tatsächlich fand er, daß die Definition der mittleren Bestrahlung eines Breitenkreises mehrdeutig ist und daß man zu gegenseitigen Widersprüchen kommt, sobald man über einen Tag hinausgeht. Damit hatte Hopfner geradezu in ein Wespennest der alten Auffassungen hineingestochen und es wurde auch sofort mit Berichtigungen über ihn hergefallen. Am einfachsten aber hat es M. Milankovitch gemacht, indem er Hopfners Darlegungen kurzer Hand als falsch bezeichnete, ohne es aber zu beweisen.”

Translated from German into English:

“A treatise by Chr. Wiener gained particular importance when the problem became important for the study of the climatic conditions of the Ice Age, when he found that at the time of the summer solstice the daily irradiation of the pole is 4/3 of that at the same time at the equator. This result has also found its way into the handbooks and textbooks of meteorology and has aroused great interest in research into the polar climate in the Tertiary period.

Therefore, some time ago I caused my former student Fr. Hopfner to subject the problem of radiation to a thorough mathematical revision, and indeed he found that the definition of the mean irradiance of a parallel is ambiguous and that contradictions arise once one goes beyond a day. With this, Hopfner had literally stabbed into a wasp’s nest of the old views, and he was immediately attacked with corrections. M. Milankovitch did it the simplest way by calling Hopfner’s explanations wrong without proving it.”

As the diagrams in Figure 5 illustrate, at the time of the respective summer solstice the daily mean irradiance over the poles is about 4/3 of that at the same time at the equator. Furthermore, Milankovitch [5] eventually showed that Hopfner’s explanation is incorrect. In Section 5, we will show that Milankovitch was right.

(a) (b)

Figure 5. The daily mean solar irradiance normalized by the solar constant for different latitudes of (a) the SH and (b) the NH for the tropical year beginning with $TDB=2405602.5$ (March 20, 1874, 00:00 UT1) predicted with Equation (8). The Earth’ heliocentric distance and the Sun’s apparent geocentric declination are based on the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. In accord with Equation (1), daily mean values were computed using 144 data per day.

2.2. Criticism from EIKE Experts

Recently, the findings of Meech [1] and Wiener [2, 3] that the global average of the daily mean solar irradiance at the TOA is equal to one quarter of the solar constant later confirmed by Milankovitch [4, 5], List [6], Fortak [13], Raschke et al. [14], Liou [9], Fu [10], Berger and Yin [15] and many others were questioned by experts of EIKE, a registered association but not an academic institution².

Assuming a long-term bright side of the Earth that would require that the Earth is tidally locked to the Sun (like the Moon to the Earth [60]) and ignoring the obliquity and the eccentricity of its nearly elliptic orbit around the Sun, the EIKE experts argued that the incoming solar radiation at the TOA must only be averaged over the bright side of the Earth, centered on the equator at 25˚E (see Figure 6), to derive a so-called hemispheric Stefan-Boltzmann temperature of the Earth using concentric surface rings of 1˚ width. This consideration would also require a long-term dark side of the Earth centered on the equator at 155˚W. This configuration would mean that the area between the longitudes 115˚E and 65˚W illustrated by Figure 7 would be in darkness. Therefore, the hemispheric average would be equal to half the solar constant. The EIKE experts further argued that only this value of about 684 W∙m⁻² (see Figure 8), reduced by a reflected portion of about 214 W∙m⁻² (that corresponds to an hemispheric albedo in the solar range of ${\alpha }_h\cong 0.31$ ) and a portion of about 80 W∙m⁻² associated with the emission of infrared radiation to space via the atmospheric window must be considered (i.e., 390 W∙m⁻²) to calculate a surface temperature averaged over the bright side of the Earth leading to about 288 K, twice the global average of Gerlich and Tscheuschner [61] for their thought model of the Earth without atmosphere. By comparing this hemispheric average of the near-surface temperature with the global average of the near-surface temperature of about $\langle T_{ns}\rangle \cong 288\text{K}$ (i.e., by comparing apples with oranges), the EIKE experts concluded that their results leave no room for a natural atmospheric greenhouse effect. Any criticism of their physically and astronomically inadequate views was brusquely rejected. In addition, some EIKE authors have deliberately misquoted those scientists’ papers published in peer-reviewed scientific journals, where even diagrams used in these papers were falsified through mutilation.

To stop such awkward attempts to question scientific findings by means of physically and astronomically inadequate arguments, it is imperative to present the facts on which the calculation of incoming solar radiation at the TOA is based.

Figure 6. The bright side of the Earth according to EIKE is divided into concentric rings of one degree each around the base of the Sun at the equator at 25˚E, i.e. 0˚ - 1˚, 1˚ - 2˚, 2˚ - 3˚,..., 89˚ - 90˚ (https://eike-klima-energie.eu/2023/12/31/der-hemisphaerische-stefan-boltzmann-ansatz-ist-kein-reines-strahlungsmodell-teil-1/, retrieved on 05/27/2025). This world view would mean that ${\lambda }_S=0$ or ${\lambda }_S=\pi$ so that Sun’s apparent geocentric declination is ${\delta }_S=0$ .

Figure 7. World Ocean Atlas Climatology (decadal average 1955-2017): Distribution of the annual sea-surface temperature in ˚C at 1˚ × 1˚ [62].

(a) (b)

Figure 8. Comparison of (a) the Earth’s annual global mean energy budget based on the present study of Kiehl and Trenberth [63] in units of W⋅m⁻², and (b) the manipulated diagram presented by EIKE Vice president Limburg on the 15^th International EIKE Conference Climate and Energy, Braunsbedra, Germany, November 25-26, 2022 (https://eike-klima-energie.eu/15-internationale-klima-und-energiekonferenz/, retrieved on 05/27/2025). Translation of the text in purple “Greenhouse effect: A nil sum play … as long as the difference is 86 W⋅m⁻² any value is allowed.” and black: “1/2 albedo reflection.”

3. Further Astronomical aspects

As the Earth is not a sphere but an oblate spheroid and because of the obliquity, i.e., the tilt of the Earth’s rotational axis with respect to the normal vector, $n$ , of the ecliptic plane pointing to the ecliptic pole (see Figure 1), mainly the gravitational forces of the Sun and the Moon cause a torque on the Earth’s equatorial bulge leading to a small temporal change in the angular momentum. This means that the assumption that the angular momentum in a central field is a conservative quantity is not exactly fulfilled. This torque tries to align the Earth’s rotational axis parallel to $n$ [64, 65]. However, like in the case of a spinning toy top on which a torque is acting, the Earth’s rotational axis traces out a cone around the pole of the ecliptic (see Figure 1) in a cycle of about 26,000 years [66]. This behavior is called the luni-solar precession. Because the Sun and the Moon change their positions relative to each other, their gravitational forces also cause a nutation of the Earth’s rotational axis, where the principal period of nutation is 18.6 years [66]. This effect is much smaller in magnitude than the luni-solar precession. Nevertheless, both precession and nutation must be considered to calculate the Terrestrial Intermediate Origin (TIO) and the Celestial Intermediate Origin (CIO), the apparent geocentric position of the Sun, and the sidereal time at Greenwich or another location.

A closed elliptic orbit (ideally characterized by Equations (47) to (54)) requires that the gravitational potential reciprocally depends on $r$ (see the 2^nd term on the right-hand side of Equation (50)). Deviations from that due to the perturbations of the gravity field by other planets lead to an open orbit of a rosette-like shape (see Figure 9). The Earth’s orbit seems to move around the Sun, resulting in a precession of the perihelion. The precession of the orbital plane or the equator of a rotating body is a secular motion. In contrast, nutation is an oscillation of the rotation pole of a freely rotating body that is undergoing torque from external gravitational forces. Nutation of the Earth’s pole is specified in terms of components in obliquity and longitude [47]. According to the recommendation of the International Astronomic Union, Division I Working Group on Precession and the Ecliptic published by Hilton et al. [67], luni-solar precession and planetary precession should be replaced by the precession of the equator and the precession of the ecliptic for general use. Both precession phenomena are still subsumed under the notion of ‘general precession’. The combination of the general precession and the precession of the perihelion is called the climatic precession, and the related parameter $e\sin \varpi$ is called the climatic precession parameter. Figure 10 presents a sketch of the combined effect of these precession phenomena. Today, the North Pole tilts away from the Sun at perihelion during the southern summer. On the contrary, 11,000 years ago, the North Pole tilted towards the Sun at perihelion during the northern summer.

The speed of light is finite. It is about $c\cong 2.9879\times {10}^8\text{m}\cdot {\text{s}}^{-1}$ in vacuum. Thus, sunlight needs, at least, 492 s to travel to the Earth. It is called the light-time. Since the orbital velocity of the Earth is, at least, 29.3 km∙s⁻¹ (see Figure 2(c)), the Earth would move in its orbit, at least, by 14,416 km during the light-time. Thus, it is unavoidable to consider the aberration of light due to the Earth’s motion.

According to Kaplan [68], there are a number of calculations that must be performed to obtain observables. These begin with astrometric reference data: a precomputed solar system ephemeris and, if a star is involved, a star catalog with positions and proper motions listed for a specified epoch. The computations account for the space motion of the object (star or planet), parallax (for a star) or light-time (for a planet), gravitational deflection of light, and the aberration of light due to the Earth’s motions. This means that the computation of the observable “solar irradiance at the TOA” should also account for precision, nutation, Earth’s rotation, polar motion, aberration of light etc. As pointed out by Kaplan [68], there are classical expressions for all these effects (except gravitational deflection), and relativity explicitly enters the procedure in only a few places, usually as added terms to the classical expressions and in the formulas that link the various time scales used. It has become common, then, to view this ensemble of calculations as being carried out entirely in a single reference system; or, two reference systems, barycentric and geocentric, that have parallel axes and differ only in the origin of coordinates (that is, they are connected by a Galilean transformation). For example, the coordinate system defined by the “equator and equinox of J2000.0”, can be thought of as either barycentric or geocentric. The relativistic effects then are interpreted simply as “corrections” to the classical result.

The March equinox plays an important role in astronomy. It is the point of origin for two different commonly used celestial coordinates: equatorial coordinates and ecliptic coordinates (see Figure 11). The line of intersection of the mean plane of the equator and the ecliptic defines the direction of the equinox. Using this direction as the origin, the right ascension, $\alpha$ , is measured in the plane of the equator and the celestial (or ecliptic) longitude, $\lambda$ , is measured in the plane of the ecliptic. Both right ascension and longitude are measured in the positive (or right-handed) sense. Like the hour angle, the right ascension is expressed in time from 0 h to 24 h. The declination, $\delta$ , is measured from the equatorial plane, positive to the north, from 0˚ to 90˚, and the celestial latitude, $\beta$ , is measured from the ecliptic plane, positive to the north, from 0˚ to 90˚ [48].

Figure 9. Open orbit of a rosette-like shape and the precession of the perihelion. Here, $r_P$ is the radius of the circle on which the perihelion is advancing by an angle of $\Delta \vartheta$ , and $r_A$ is the radius on which the aphelion is moving forward by $\Delta \vartheta$ (with reference to [11, 69]).

Figure 10. Combined effect of the precession phenomena (with reference to Crowley and North [44]). The letter P stands for perihelion.

The vector $r$ from the origin, $O$ , to a celestial object may be represented by rectangular coordinates ($x,y,z$ ); for instance, in a Cartesian vector space by $r=xi+yj+zk$ , where $i$ , $j$ , and $k$ are the unit vectors in the direction of $x$ , $y$ , and $z$ . These unit vectors form a right-handed rectangular coordinate frame (aka right-handed trihedron). Coordinate frames with covariant and contravariant bases may also be used. The covariant and the contravariant basis vectors are defined (using Einstein’s summation convention)

by $q_i=\partial r/\partial q^i$ and $q^i=\partial r/\partial q_i$ , $i=1,2,3$ , respectively. Here, $q^i$ and $q_i$ , $i=1,2,3$ , are the contravariant and the covariant coordinates, respectively.

The position vector of the celestial object may also be represented by spherical coordinates, for instance, by $r=r\cos {\rm B}\cos \Gamma i+r\cos {\rm B}\sin \Gamma j+r\sin {\rm B}k$ in which the direction is specified by the longitudinal angle, $\Gamma$ , in the $xy$ -reference plane and the latitudinal angle, ${\rm B}$ , from the reference plane to the position vector. Thus, the Cartesian coordinates are given by $x=r\cos {\rm B}\cos \Gamma$ , $y=r\cos {\rm B}\sin \Gamma$ , and $z=r\sin {\rm B}$ . Therefore, using right ascension, $\alpha$ , and declination, $\delta$ , the equatorial coordinates can be expressed by [70-72] (see Figure 11).

Figure 11. Equatorial and ecliptic reference planes (adopted from [66, 73]). The star denotes the celestial object in question.

$\begin{array}{lll}{x}_{eq}\hfill & =\hfill & r\cos \alpha \cos \delta \hfill \\ y_{eq}\hfill & =\hfill & r\sin \alpha \cos \delta \hfill \\ z_{eq}\hfill & =\hfill & r\sin \delta \hfill \end{array}$ (17)

where the actual distance, $r$ , is given by

$r=+\sqrt{x_{eq}^2+y_{eq}^2+z_{eq}^2}$ .(18)

Since

$\frac{y_{eq}}{x_{eq}}=\tan \alpha$ ,(19)

the right ascension is given by

$\alpha =\arctan \left(\frac{y_{eq}}{x_{eq}}\right)=\text{atan}2\left(y_{eq},x_{eq}\right)$ ,(20)

where atan2 is the four-quadrant inverse tangent. Because of

$\frac{z_{eq}}{r}=\sin \delta \text{and}\frac{z_{eq}}{+\sqrt{x_{eq}^2+y_{eq}^2}}=\tan \delta$ ,(21)

the declination is given by

$\delta =\arcsin \left(\frac{z_{eq}}{r}\right)=\arctan \left(\frac{z_{eq}}{+\sqrt{x_{eq}^2+y_{eq}^2}}\right)$ .(22)

As mentioned before, the equator and the ecliptic are moving because of the effects of perturbing forces on the rotation and motion of the Earth. Hence, the March equinox and the obliquity change with time. Therefore, these celestial coordinate frames must be carefully defined to relate to a standard frame that may be regarded as being fixed in space. Based on the geocentric rectangular coordinates provided by the JPL ephemeris DE440 [17], the right ascension, ${\alpha }_S$ , and the declination, ${\delta }_S$ , of the Sun can, in principle, be computed using Equations (20) and (22), respectively. However, first, an adjustment of the position vector $r$ for the aberration of light due to the Earth’s motion must be considered. Then, the reductions of the geocentric position $r_{GCRS}$ with respect to the Geocentric Celestial Reference System (GCRS) to a position $r_t$ with respect to the equator and the equinox of date for frame bias (GCRS to J2000.0), performed by matrix $B$ , precession performed by matrix $P$ , and nutation performed by matrix $N$ are required [68, 72], i.e.,

$r_t=NPB{r}_{GCRS}=M{r}_{GCRS}$ (23)

With $M=NPB$ . These matrices are discussed, for instance, in the Astronomical Almanac for the Year 2023. As mentioned before, we used the NOVAS F3.1 subroutines for aberration, frame bias, precession and nutation (Kaplan et al. [18, 19]) to calculate $r_t$ .

The ecliptic coordinates $x_{ecl}$ , $y_{ecl}$ , and $z_{ecl}$ can be calculated similarly, but under consideration of the celestial longitude, $\lambda$ , and the celestial latitude, $\beta$ (see Figure 11). The ecliptic coordinates and the equatorial coordinates are related to each other by

$\begin{array}{lll}{x}_{ecl}\hfill & =\hfill & x_{eq}\hfill \\ y_{ecl}\hfill & =\hfill & y_{eq}\cos \epsilon +z_{eq}\sin \epsilon \hfill \\ z_{ecl}\hfill & =\hfill & -y_{eq}\sin \epsilon +z_{eq}\cos \epsilon \hfill \end{array}$ (24)

Replacing $x_{eq}$ , $y_{eq}$ , and $z_{eq}$ according to Equation (17) gives for the ecliptic coordinates with respect to the equinox (either the mean or the true equinox) as the origin for the ecliptic longitude and latitude [70, 72]

$\cos \beta \cos \lambda =\cos {\alpha }_e\cos \delta$ ,(25)

$\frac{y_{ecl}}{x_{ecl}}=\tan \lambda =\frac{\sin \delta \sin \epsilon +\sin {\alpha }_e\cos \delta \cos \epsilon }{\cos {\alpha }_e\cos \delta }$ ,(26)

leading to

$\cos \beta \sin \lambda =\sin \delta \sin \epsilon +\sin {\alpha }_e\cos \delta \cos \epsilon$ ,(27)

and

$\frac{z_{ecl}}{r}=\sin \beta =-\sin {\alpha }_e\cos \delta \sin \epsilon +\sin \delta \cos \epsilon$ .(28)

We similarly obtain

$\frac{y_{eq}}{x_{eq}}=\tan {\alpha }_e=\frac{-\sin \beta \sin \epsilon +\cos \beta \cos \epsilon \sin \lambda }{\cos \beta \cos \lambda }$ ,(29)

leading to

$\sin {\alpha }_e\cos \delta =-\sin \beta \sin \epsilon +\cos \beta \cos \epsilon \sin \lambda$ ,(30)

and

$\frac{z_{eq}}{r}=\sin \delta =\sin \beta \cos \epsilon +\cos \beta \sin \epsilon \sin \lambda$ .(31)

If we assume $\beta \approx 0$ , Equation (31) provides [4, 5]

$\sin \delta =\sin \lambda \sin \epsilon$ (32)

on which Equation (7) is based. Note that $\lambda =\overline{\lambda }+\Delta \psi$ , where $\overline{\lambda }$ is the mean longitude and $\Delta \psi$ is the nutation in the longitude. Furthermore, $\epsilon$ is either the mean or the true obliquity of the ecliptic, where $\overline{\epsilon }$ denotes the mean obliquity of date expressed by [17, 74]

$\overline{\epsilon }=84381".448-46".815T-0".00059T^2+0".001813T^3$.(33)

The true obliquity of date is given by $\epsilon =\overline{\epsilon }+\Delta \epsilon$ , where $\Delta \epsilon$ is the nutation in the obliquity. Both, $\Delta \psi$ and $\Delta \epsilon$ may be parameterized by $\Delta \psi =-19".1996\sin {\Omega }_{a,M}$ and $\Delta \epsilon =9".2025\cos {\Omega }_{a,M}$ [17], where ${\Omega }_{a,M}$ is the ascending node of the Moon’s orbit on the ecliptic given by [17, 74]

${\Omega }_{a,M}=125˚02\text{'}40".280-1934˚08\text{'}10".549T+7".455T^2+0".008T^3$ ,(34)

where $T$ is the Barycentric Dynamical Time (TDB) in centuries with respect to the epoch J2000.0, i.e., $JD=2451545.0$ (January 1, 2000, 12:00 UT1), where the Julian century corresponds to 36525 d. However, we used $\Delta \psi$ and $\Delta \epsilon$ directly provided by the JPL ephemeris DE440. The differences between the parameterized and exactly calculated values are shown for 2010 in Figure 12.

(a) (b)

Figure 12. Nutation (a) in longitude and (b) in obliquity provided by the JPL ephemeris DE440 compared with the approximations $\Delta \psi =-19".1996\sin {\Omega }_{a,M}$ and $\Delta \epsilon =9".2025\cos {\Omega }_{a,M}$ mentioned by Park et al. [17].

Figure 13. Relationships of origins (adopted from Hohenkerk, [75]). See text for discussion.

The equinox right ascension is related to the intermediate right ascension, ${\alpha }_i$ , by [72]

${\alpha }_e={\alpha }_i-E_o\left(T\right)$ ,(35)

where $E_o\left(T\right)=E_{prec}\left(T\right)-E_e\left(T\right)$ is the equation of the origins. Here, $E_{prec}\left(T\right)$ is the accumulated precession and $E_e\left(T\right)$ is the equation of the equinoxes representing the accumulated nutation. The precession portion is given in units of arcseconds by [76]

$E_{prec}\left(T\right)=-0".014506-4612".156534T-1".3915817T^2+0".00000044T^3+0".000029956T^4$ .(36)

The equation of the equinoxes (expressed in arcseconds) is given by [68, 76]

$\begin{array}{l}{E}_e\left(T\right)=\Delta \psi \cos \overline{\epsilon }+0.00264096\sin {\Omega }_{a,M}+0.00006352\sin 2{\Omega }_{a,M}\\ +0.00001175\sin \left(2F-2D+3{\Omega }_{a,M}\right)+0.00001121\sin \left(2F-2D+{\Omega }_{a,M}\right)\\ -0.00000455\sin \left(2F-2D+2{\Omega }_{a,M}\right)+0.00000202\sin \left(2F+3{\Omega }_{a,M}\right)\\ +0.00000198\sin \left(2F+{\Omega }_{a,M}\right)-0.00000172\sin 3{\Omega }_{a,M}-0.00000087T\sin {\Omega }_{a,M}.\end{array}$(37)

The fundamental luni-solar $F$ and $D$ can be found, for instance, in the Astronomical Almanac for the Year 2023, page B47. However, for our purposes, the equation of the equinoxes (expressed in seconds) can be approximated with sufficient accuracy by

$E_e\left(T\right)=\left(\Delta \psi \cos \overline{\epsilon }+0.002641\sin {\Omega }_{a,M}+0.000064\sin 2{\Omega }_{a,M}\right)/15$ ,(38)

The time can be obtained by use of the Earth Rotation Angle (ERA) $\theta$ (see Figure 13), i.e., the angle between the non-rotation origins defining the Celestial Intermediate Origin (CIO) and the Terrestrial Intermediate Origin (TIO) to realize the intermediate reference frame of epoch $t$ , expressed as a function of UT1 by [68, 76]

$\theta \left(D_U\right)=2\pi \left(0.7790572732640+1.00273781191135448D_U\right)$ ,(39)

where $D_U$ is the Julian UT1 date minus 2451545.0. We used the subroutine EROT of NOVAS F3.1 (Kaplan et al. [18, 19]) to calculate ERA $\theta$ . Using ERA, the Greenwich Apparent Sidereal Time (GAST) can be expressed by [76]

$GAST\left(D_U,T\right)=\theta \left(D_U\right)-E_o\left(T\right)$ (40)

and the Greenwich Mean Sidereal Time (GMST) is given by

$GMST\left(D_U,T\right)=\theta \left(D_U\right)-E_{prec}\left(T\right)$ (41)

so that

$GAST\left(D_U,T\right)=GMST\left(D_U,T\right)+E_e\left(T\right)$ .(42)

To compute the CIO locator, $s$ , the vector components of the Celestial Intermediate Pole (CIP), $x_{CIP}$ and $y_{CIP}$ , are required. They can be determined from the matrix elements $x_{CIP}=M_{1,3}$ and $y_{CIP}=M_{2,3}$ (in our notation) so that

$s=-x_{CIP}\cdot y_{CIP}/2+\cdots$ ,(43)

where the dots represent lengthy series (for more details, see Capitaine et al. [77]). Tables for the series for $x_{CIP}$ , $y_{CIP}$ and $s+x_{CIP}\cdot y_{CIP}/2$ are available only in electronic form, at the Strasbourg Astronomical Data Center (CDS) via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/400/1145.

Unless explicitly stated otherwise, we will refer to the Sun’s apparent geocentric right ascension and declination and the apparent ecliptic longitude and latitude.

Perihelion and aphelion can simply be determined using the first and second derivative tests [36, 37]. In accord with Equation (18), the first derivative of the actual geocentric distance with respect to time is given by

$\frac{\text{d}r}{\text{d}t}=\frac{1}{r}\left(x_{eq}\frac{\text{d}{x}_{eq}}{\text{d}t}+y_{eq}\frac{\text{d}{y}_{eq}}{\text{d}t}+z_{eq}\frac{\text{d}{z}_{eq}}{\text{d}t}\right)$ .(44)

As mentioned before, the JPL ephemeris DE440 provides the geocentric coordinates and their rates needed for computing this derivative. The optimum is given if the expression in the parentheses is equal to

zero (i.e., $\text{d}r/\text{d}t\cong 0$ ) [36, 37]. The sign of the second derivative for the time of the optimum marks either the perihelion (${\text{d}}^{2}r/\text{d}{t}^2>0$ ) or the aphelion (${\text{d}}^{2}r/\text{d}{t}^2<0$ ). Because the ephemeris does not deliver ${\text{d}}^2{x}_{eq}/\text{d}{t}^2$ , ${\text{d}}^2{y}_{eq}/\text{d}{t}^2$ , and ${\text{d}}^2{z}_{eq}/\text{d}{t}^2$ , we determined the second derivative numerically from the $\text{d}r/\text{d}t$

curve. Thus, Equation (44) can be used to determine the anomalistic year, i.e., the interval between successive passages of the Earth through the perihelion.

In accord with Equation (21), the Sun’s geocentric declination can be determined using

$\tan {\delta }_S=\frac{z_{eq}}{+\sqrt{x_{eq}^2+y_{eq}^2}}$ .(45)

From Equation (45), we can infer the summer and winter solstice using the first and second derivative tests. The first derivative reads

$\frac{\text{d}{\delta }_S}{\text{d}t}={\cos }^2{\delta }_S\frac{\left(x_{eq}^2+y_{eq}^2\right)\frac{\text{d}{z}_{eq}}{\text{d}t}-z_{eq}\left(x_{eq}\frac{\text{d}{x}_{eq}}{\text{d}t}+y_{eq}\frac{\text{d}{y}_{eq}}{\text{d}t}\right)}{{\left(x_{eq}^2+y_{eq}^2\right)}^{3/2}}$ .(46)

Optimum, minimum (winter solstice), and maximum (summer solstice) can be deduced using $\text{d}{\delta }_S/\text{d}t\cong 0$ , ${\text{d}}^2{\delta }_S/\text{d}{t}^2>0$ , and ${\text{d}}^2{\delta }_S/\text{d}{t}^2<0$ , respectively. The $\text{d}{\delta }_S/\text{d}t$ curve serves to compute the second derivative numerically. If ${\text{d}}^2{\delta }_S/\text{d}{t}^2\cong 0$ , a point of inflection may occur. This condition, however, is

only a necessary condition; the sufficient condition requires that the 2^nd derivatives in the neighborhood about this point have opposite signs. Figure 14 illustrates the Sun’s geocentric declination as well as the first and second derivatives, where the latter is multiplied by a factor of 20.

For assessing the results of our astronomical calculations, the results obtained for perihelion and aphelion for or around the years 1610, 1650, 1710, 1750, 1810, 1850, 1853, 1874, 1910, 1950, 2009, and 2010 are listed in Table 2 and Table 3. These tables include the Earth’s heliocentric distance and the distance between the barycenter of the solar system and the Earth. In addition, the results obtained for the March and September equinoxes and the June and December solstices are listed in Table 4 and Table 5, respectively. The perihelion moves from December 28, 1609, to January 3, 2010, leading to a shift of the beginning of the anomalistic year by about six days. In 1874, the Earth reached the perihelion twice. The aphelion moves from June 27, 1610, to July 6, 2010. Our results also confirm that the period from the March equinox to the subsequent September equinox is approximately 7.6 days longer than the period from this September equinox to the subsequent March equinox.

Figure 14. Sun’s geocentric declination, ${\delta }_S$ , the first derivative, $\text{d}{\delta }_S/\text{d}t$ , and the second derivative ${\text{d}}^2{\delta }_S/\text{d}{t}^2$ multiplied by a factor of 20 for the year 1610 beginning with $TDB=2309100.5$ that corresponds to January 1, 1610, 00:00 UT1. The required data were provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19].

Table 2. Perihelion in the reference frame ICRF of the JPL ephemeris DE440 [17] for or around the years 1610, 1650, 1710, 1750, 1810, 1850, 1853, 1874, 1910, 1950, 2009, and 2010.

Julian day	Date, time in UT1	GMST in hh:mm:ss.f	GAST in hh:mm:ss.f	Heliocentric distance in AU	Solar barycenter distance in AU
2309462.0069	Dec 28, 1610, 12:10	18:36:57.328	18:36:56.289	0.98310912	0.97887891
2324073.0	Dec 29, 1650, 12:00	18:32:05.885	18:32:05.323	0.98314219	0.98648820
2345988.375	Dec 30, 1710, 21:00	03:35:25.115	03:35:25.858	0.98315473	0.98425182
2360598.4514	Dec 30, 1750, 22:50	05:26:56.904	05:26:57.980	0.98324605	0.98260243
2382146.8264	Dec 30, 1809, 07:50	14:23:20.429	14:23:20.807	0.98323720	0.98593708
2382514.7014	Jan 2, 1811, 04:50	11:33:43.231	11:33:43.270	0.98323253	0.98402576
2397123.3889	Dec 31, 1850, 21:20	03:59:46.544	03:59:45.718	0.98328535	0.98343918
2397853.2083	Dec 30, 1852, 17:00	23:37:09.245	23:37:08.233	0.98318820	0.98623264
2398220.2292	Jan 1, 1854, 17:30	00:14:09.991	00:14:09.072	0.98323129	0.98653161
2405525.0903	Jan 1, 1874, 14:10	20:54:14.045	20:54:13.386	0.98322835	0.98292504
2405888.8403	Dec 31, 1874, 08:10	14:48:21.057	14:48:20.682	0.98325690	0.98516225
2418672.9514	Jan 1, 1910, 10:50	17:30:51.088	17:30:50.179	0.98326210	0.98134988
2433284.7639	Jan 3, 1950, 06:20	13:09:13.709	13:09:13.529	0.98322745	0.98543732
2454836.1458	Jan 4, 2009, 15:30	22:27:28.762	22:27:29.569	0.98327303	0.98816647
2455199.5069	Jan 3, 2010, 00:10	07:00:03.783	07:00:04.805	0.98328967	0.98694532

Table 3. Aphelion in the reference frame ICRF of the JPL ephemeris DE440 [17] for the years 1610, 1650, 1710, 1750, 1810, 1850, 1853, 1874, 1910, 1950, 2009, and 2010.

Julian day	Date, time in UT1	GMST in hh:mm:ss.f	GAST in hh:mm:ss.f	Heliocentric distance in AU	Solar barycenter distance in AU
2309277.5486	Jun 27, 1610, 01:10	19:29:42.723	19:29:41.710	1.01686460	1.02096020
2323889.2431	Jun 28, 1650, 17:50	12:17:37.198	12:17:36.488	1.01690760	1.01446020
2345803.9931	Jun 29,1710, 11:50	06:18:28.580	06:18:29.210	1.01684910	1.01636060
2360416.0069	Jul 1, 1750, 12:10	06:47:38.693	06:47:39.762	1.01683060	1.01632690
2382330.7361	Jul 2, 1810, 05:40	00:18:25.258	00:18:25.481	1.01681200	1.01514460
2396942.3889	Jul 3, 1850, 21:20	16:06:10.022	16:06:09.306	1.01676020	1.01749750
2398037.7014	Jul 3, 1853, 04:50	23:34:32.067	23:34:31.063	1.01681430	1.01347020
2405708.4722	Jul 3, 1874, 23:20	18:07:14.026	18:07:13.520	1.01673470	1.01599770
2418857.4375	Jul 4, 1910, 22:30	17:18:12.267	17:18:11.432	1.01675660	1.01777940
2433468.4375	Jul 5, 1950, 22:30	17:23:22.690	17:23:22.690	1.01671540	1.01508630
2455016.5694	Jul 4, 2009, 01:40	20:28:48.935	20:28:49.849	1.01666640	1.01227170
2455383.9792	Jul 6, 2010, 11:30	06:27:21.678	06:27:22.720	1.01670200	1.01392020

One may also consider Keplerian elements to compute the heliocentric distance

$r=\frac{p}{1+e\cos \upsilon }=\frac{a\left(1-e^2\right)}{1+e\cos \upsilon }$ .(47)

Here, $p=L^2/\left(m\alpha \right)$ , $e={\left\{1+2E{L}^2/\left(m{\alpha }^2\right)\right\}}^{1/2}$ is the eccentricity, $\upsilon$ is the true anomaly (see Figure 1), $\alpha ={\gamma }_{G}Mm$ , where ${\gamma }_G$ is the gravitational constant, $M$ is the mass of the Sun, $m$ is the mass of the Earth, and $L=m{r}^2\text{d}\upsilon /\text{d}t=const.$ is the angular momentum considered as invariant with time, i.e., the angular momentum in a central field like Newton’s gravity field is a conservative quantity. The quantity $2p$ is called the latus rectum.

Table 4. March and September equinoxes for the years 1610, 1650, 1710, 1710, 1750, 1810, 1850, 1853, 1874, 1910, 1950, 2009, and 2010. Also listed are the Sun’s apparent geocentric position (true equator and equinox of date) expressed by the right ascension, declination and heliocentric distance.

Julian Day	Date, Time in UT1	Sun’s apparent geocentric position			GMST in hh:mm:ss.f	GAST in hh:mm:ss.f
		Right ascension in hh:mm:ss.f	Declination in ˚: ‘: “.f	Heliocentric distance in AU
2309179.2986	Mar 20, 1610, 19:10	00:00:00.090	00:00:00.489	0.99791725	07:02:21.160	07:02:20.155
2309365.8125	Sep 23, 1610, 07:30	11:59:59.302	00:00:04.694	1.00156840	07:37:42.018	07:37:40.952
2323788.9861	Mar 20, 1650, 11:40	23:59:59.523	−00:00:02.924	0.99769393	23:32:20.879	23:32:20.080
2323975.5069	Sep 23, 1650, 00:10	12:00:00.559	−00:00:03.131	1.00174920	00:17:43.381	00:17:42.716
2345703.5278	Mar 21, 1710, 00:40	00:00:00.821	00:00:04.525	0.99741925	12:32:22.981	12:32:23.507
2345890.0208	Sep 23, 1710, 12:30	12:00:00.509	−00:00:03.822	1.00206120	12:37:38.910	12:37:39.558
2360313.2153	Mar 20, 1750, 17:10	23:59:59.815	−00:00:02.244	0.99715147	05:02:22.774	05:02:23.836
2360499.7014	Sep 23,1750, 04:50	12:00:00.062	−00:00:00.239	1.00225470	04:57:37.062	04:57:38.112
2382227.7639	Mar 21, 1810, 06:20	23:59:59.904	−00:00:00.614	0.99688612	18:12:26.629	18:12:26.924
2382414.2361	Sep 23, 1810, 17:40	12:00:00.451	−00:00:03.768	1.00255920	17:47:37.630	17:47:37.741
2396837.4583	Mar 20, 1850, 23:00	23:59:59.731	−00:00:02.283	0.99673925	10:52:28.139	10:52:27.484
2397023.9167	Sep 23, 1850, 10:00	11:59:59.969	−00:00:00.330	1.00267980	10:07:35.857	10:07:35.047
2397933.1806	Mar 20, 1853, 16:20	23:59:59.412	−00:00:03.434	0.99670902	04:12:27.103	04:12:26.077
2398119.6528	Sep 23, 1853, 03:40	12:00:00.576	−00:00:03.936	1.00273860	03:47:38.107	03:47:37.082
2405603.2778	Mar 20, 1874, 18:40	00:00:00.328	00:00:02.151	0.99664796	06:32:29.717	06:32:29.081
2405789.7222	Sep 23, 1874, 05:20	11:59:59.570	00:00:02.242	1.00280420	05:27:34.150	05:27:33.669
2418752.0	Mar 21, 1910, 12:00	23:59:59.537	−00:00:02.907	0.99640641	23:52:30.461	23:52:29.575
2418938.4375	Sep 23, 1910, 22:30	11:59:59.878	00:00:01.302	1.00301590	22:37:33.251	22:37:32.416
2433361.6944	Mar 21, 1950, 04:40	00:00:00.669	00:00:04.362	0.99624619	16:32:32.046	16:32:31.902
2433548.1111	Sep 23, 1950, 14:40	11:59:59.400	00:00:03.876	1.00322720	14:47:29.908	14:47:29.951
2454910.9931	Mar 20, 2009, 11:50	00:00:00.812	00:00:04.791	0.99596760	23:42:34.276	23:42:35.125
2455097.3889	Sep 22, 2009, 21:20	12:00:00.041	−00:00:00.573	1.00354280	21:27:27.209	21:27:28.120
2455276.2292	Mar 20, 2010, 17:30	23:59:59.494	−00:00:03.113	0.99594656	05:22:32.837	05:22:33.826
2455462.6319	Sep 23, 2010, 03:10	12:00:00.005	00:00:00.987	1.00349940	03:17:27.415	03:17:28.444

The Earth’s elliptic orbit around the Sun, characterized by Kepler’s first law that “the orbit of each planet is an ellipse and the Sun is at one of the two foci”, is a consequence of the state of energy in this central field expressed by [11, 65, 69]

$\frac{L^2}{m\alpha r}=1+{\left\{1+\frac{2E{L}^2}{m{\alpha }^2}\right\}}^{1/2}\cos \upsilon$ ,(48)

where

$E=T_{radial}+U_{eff}\left(r\right)$ (49)

is the total energy, and

$U_{eff}\left(r\right)=L^2/\left(2m{r}^2\right)-\alpha /r$ (50)

is the effective potential comprising the centrifugal potential and the gravitational potential [65, 69], and $T_{radial}=m{\left(\text{d}r/\text{d}t\right)}^2/2$ is the radial kinetic energy (equal to zero for a circle). Equation (48) leads to Equation (47) if $p$ and $e$ are inserted. Choosing $\upsilon =0˚$ leads to $r_p=p/\left(1+e\right)$ , and, hence, $r_p=a\left(1-e\right)$ the minimum of the Earth’ heliocentric distance at the perihelion. Choosing $\upsilon =180˚$ leads to $r_a=p/\left(1-e\right)$ and, hence, $r_a=a\left(1+e\right)$ the maximum of the Earth’ heliocentric distance at the aphelion. Thus, $p$ and $e$ are given by $p=a\left(1-e^2\right)$ and $e=\left(a^2-b^2\right)/a$ as already used in Equation (47). Therefore, the eccen-tricity may be expressed by $e=\left(r_a-r_p\right)/\left(r_a+r_p\right)=\left(r_a-r_p\right)/\left(2a\right)$ [35]. Furthermore, in accord with Equation (47), we define the relative heliocentric distance [4, 5, 35]

$\rho =\frac{r}{a}=\frac{1-e^2}{1+e\cos \upsilon }$ .(51)

Table 5. As in Table 4, but for the June and December solstices.

Julian Day	Date, Time in UT1	Sun’s apparent geocentric position			GMST in hh:mm:ss.f	GAST in hh:mm:ss.f
		Right ascension in hh:mm:ss.f	Declination in ˚:’:”.f	Heliocentric distance in AU
2309272.3472	Jun 21, 1610, 20:20	06:00:02.022	23:29:21.937	1.01678880	14:19:12.306	14:19:11.276
2309455.3889	Dec 21, 1610, 21:20	18:00:00.893	−23:29:24.197	0.98324126	03:20:51.792	03:20:50.750
2323882.0139	Jun 21, 1650, 12:20	06:00:01.612	23:29:10.784	1.01675420	06:19:07.099	06:19:06.372
2324065.0972	Dec 21,1650, 14:20	18:00:01.001	−23:29:12.133	0.98327686	20:20:56.440	20:20:55.828
2345796.5	Jun 22, 1710, 00:00	06:00:00.806	23:28:44.072	1.01667670	17:58:56.057	17:58:56.633
2345979.6528	Dec 22, 1710, 03:40	18:00:00.247	−23:28:42.581	0.98335070	09:41:01.827	09:41:02.569
2360406.1597	Jun 21, 1750, 15:50	06:00:00.924	23:28:17.756	1.01657710	09:48:49.280	09:48:50.338
2360589.3681	Dec 21, 1750, 20:50	17:59:59.921	−23:28:16.441	0.98340820	02:51:08.193	02:51:09.226
2382320.6597	Jun 22, 1810, 03:50	05:59:59.010	23:27:40.658	1.01651970	21:48:41.634	21:48:41.850
2382503.9375	Dec 22, 1810, 10:30	17:59:59.232	−23:27:39.624	0.98347769	16:31:16.975	16:31:16.979
2396930.3333	Jun 21, 1850, 20:00	06:00:00.223	23:27:24.664	1.01645340	13:58:38.219	13:58:37.484
2397113.6458	Dec 22, 1850, 03:30	17:59:58.627	−23:27:24.983	0.98349528	09:31:21.772	09:31:20.938
2398026.0556	Jun 21, 1853, 13:20	05:59:59.637	23:27:31.205	1.01647600	07:18:37.182	07:18:36.160
2398209.3819	Dec 21, 1853, 21:10	17:59:59.868	−23:27:32.370	0.98345850	03:11:24.022	03:11:23.052
2405696.125	Jun 21, 1874, 15:00	05:59:58.777	23:27:27.738	1.01645800	08:58:33.225	08:58:32.672
2405879.4722	Dec 21, 1874, 23:20	17:59:59.630	−23:27:27.226	0.98350526	05:21:24.993	05:21:24.579
2418844.8194	Jun 22, 1910, 07:40	05:59:58.465	23:27:07.380	1.01635730	01:38:27.398	01:38:26.528
2419028.2222	Dec 22, 1910, 17:20	18:00:01.498	−23:27:08.362	0.98358801	23:21:32.310	23:21:31.565
2433454.4792	Jun 21, 1950, 23:30	05:59:58.877	23:26:53.686	1.01635120	17:28:20.770	17:28:20.723
2433637.9236	Dec 22, 1950, 10:10	17:59:59.299	−23:26:52.435	0.98362067	16:11:35.538	16:11:35.646
2455003.7361	Jun 21, 2009, 05:40	05:59:58.852	23:26:21.595	1.01627750	23:38:13.142	23:38:14.010
2455187.2431	Dec 21, 2009, 17:50	18:00:00.389	−23:26:19.384	0.98375783	23:51:42.696	23:51:43.668
2455368.9792	Jun 21, 2010, 11:30	06:00:00.084	23:26:17.578	1.01622740	05:28:13.347	05:28:14.359
2455552.4931	Dec 21, 2010, 23:50	18:00:01.930	−23:26:15.798	0.98370865	05:51:44.543	05:51:45.594

The mean heliocentric distance $r_0$ may be determined in accord with Kepler’s second law that “the radius vector drawn from the Sun’s center to the center of the planet sweeps out equal areas in equal times”.

The period $T$ of one revolution of the Earth around the Sun is $T=2mA/L$ , where $A=\pi ab=\pi a^2{\left(1-e^2\right)}^{1/2}$ is the area of the elliptic orbit, and $b=a{\left(1-e^2\right)}^{1/2}$ is the semi-minor axis.

Thus, we may write

$r^2\frac{\text{d}\upsilon }{\text{d}t}=\frac{2\pi }{T}{a}^2{\left(1-e^2\right)}^{1/2}$ . (52)

Considering $a$ and $e$ as constant with time during this revolution, the integration of Equation (52) yields [9, 11]

${\displaystyle \underset{0}{\overset{2\pi }{\int }}{r}^2\text{d}\upsilon }=\frac{2\pi }{T}{a}^2{\left(1-e^2\right)}^{1/2}{\displaystyle \underset{0}{\overset{T}{\int }}\text{d}t}=2\pi a^2{\left(1-e^2\right)}^{1/2}$ (53)

or

$r_0{}^2=\frac{1}{2\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}{r}^2\text{d}\upsilon }=a^2{\left(1-e^2\right)}^{1/2}$ . (54)

As the eccentricity is about $e\cong 0.01671$ (see Table 6), one may use $r_0=a{\left(1-e^2\right)}^{1/4}\cong a$ . The deviation of $a$ from $r_0$ is about 10445 km.

For determining the true anomaly and, hence, the actual distance of the Earth, the transcendental equation for the eccentric anomaly called the Kepler equation [48, 78, 79],

$f\left(E\right)=E-e\sin E-M=0$ ,(55)

must be solved iteratively, for instance, by a Newton-Raphson method

$E^{\left(n+1\right)}=E^{\left(n\right)}-\frac{f\left(E^{\left(n\right)}\right)}{f^{\prime }\left(E^{\left(n\right)}\right)},n=0,1,2,\cdots$ (56)

where $f^{\prime }\left(E^{\left(n\right)}\right)=1-e\cos E^{\left(n\right)}$ , $M=2\left(t-t_P\right)/U=L-\varpi$ is the mean anomaly, $U=2.978\times {10}^4\text{m}\cdot {\text{s}}^{-1}$ is the mean orbital velocity of the Earth, $t-t_P$ is the time since perihelion, $L$ is the mean longitude, and $\varpi$ is the longitude of the perihelion counted counterclockwise from the moving March equinox (see Figure 1). The procedure may be stopped if after k iteration steps the condition $\left|E^{\left(k\right)}-E^{\left(k-1\right)}\right|<{10}^{-6}\mathrm{deg}$ is fulfilled.

Even though the accuracy of the Kepler equation may be sufficient in several cases, we generally use the JPL ephemeris DE440 [17]. Note that the values of the Keplerian elements and their rates, with respect to the mean ecliptic and equinox of J2000.0, valid for the time interval 1800 AD - 2050 AD are listed in Table 8.10.2 by Standish and Williams [79] (see Table 6).

Table 6. Keplerian elements and their rates for the Earth-Moon Barycenter, with respect to the mean ecliptic and equinox of J2000.0, valid for the time-interval 1800 AD - 2050 AD [79].

$a$ (AU, AU/cty)	$e$ rad, rad/cty	$I$ (deg, deg/cty)	$L$ (deg, deg/cty)	$\varpi$ (deg, deg/cty)	$\Omega$ (deg, deg/cty)
1.00000261	0.01671123	−0.00001531	100.46457166	102.93768193	0.0
0.00000562	−0.00004392	−0.01294668	35999.37244981	0.32327364	0.0

The heliocentric coordinates of the EMB in its orbital plane, $r_{hc}$ , with the $x_{hc}$ -axis aligned from the focus to the perihelion are then given by [79]:

$x_{hc}=r\cos \upsilon =a\left(\cos E-e\right)$ (57)

and

$y_{hc}=r\sin \upsilon =a\sqrt{1-e^2}\sin E$ (58)

with $z_{hc}=0$ . The relation between the true anomaly and the eccentric anomaly is given by

$\tan \frac{\upsilon }{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}$ .(59)

4. The total solar irradiance

Figure 15 shows the variation of the TSI at 1 AU, often called the solar constant, and the related

uncertainty during the period from November 1978 to April 2023. The results mainly follow the sunspot number, but solar activity also plays an important role. As illustrated by this figure, the solar activity was low in 2010, and the TSI was close to $S=1361\text{W}\cdot {\text{m}}^{-2}$ . This fact is one of the reasons why Kramm et al. [53] considered 2010. The reconstruction of the TSI at 1 AU for the past 400 years is illustrated in Figure 16. This reconstruction reflects both the Maunder minimum and the Dalton minimum of the sunspot numbers.

Figure 15. “Community-Consensus TSI Composite” and the uncertainty deduced from four-decade-long space-borne TSI-measurements, created by G. Kopp on May 21, 2025 using the methodology of Dudok de Wit et al. [23] (see https://spot.colorado.edu/~koppg/TSI/TSI_Composite-SIST.txt).

Figure 16. Reconstruction of the total solar irradiance (TSI) for the past 400 years based on “Community-Consensus TSI Composite” [23] and SATIRE-T model [24] (with modifications to fix spurious values prior to 1650). The model is scaled by 0.999996 (−0.0061 W∙m⁻²) to match TSI data record of Kopp and Lean [21]. Extended using Community-Consensus TSI Composite annual averages from 1978 onward; computed by Greg Kopp using “Community-Consensus TSI Composite”.

As mentioned before, the solar constant is the TSI at the Earth’s mean heliocentric distance $r_0$ . However, the heliocentric distance of $r_0=1\text{AU}$ is considered due to the relatively small difference. At $TDB=\text{2455291}\text{.0764}$ (April 4, 2010, 13:50 UT1), $r_0=\text{1}\text{.00014990}$ corresponds to the Sun’s apparent geocentric position ${\alpha }_S={\text{00}}^{\text{h}}{\text{54}}^{\text{m}}{\text{04}}^{\text{s}}\text{.561}$ (right ascension), ${\delta }_S=5˚47\text{'}12".953$ (geocentric declination) and ${\lambda }_S\cong \text{14}˚\text{41}\text{'}\text{3}".\text{386}$ (ecliptic longitude), respectively (see Figure 2). The Earth’s heliocentric distance of $r_0=1\text{AU}$ is reached at $TDB=\text{2455290}\text{.5625}$ (April 4, 2010, 01:30 UT1) that corresponds to the Sun’s apparent geocentric position ${\alpha }_S={\text{00}}^{\text{h}}{\text{52}}^{\text{m}}{\text{11}}^{\text{s}}\text{.976}$ , ${\delta }_S=\text{5}˚\text{35}\text{'}\text{27}".\text{332}$ and ${\lambda }_S=\text{14}˚\text{10}\text{'}\text{40}".\text{793}$ , respectively. This means that the distance of 1 AU is reached near the March equinox. Because the March equinox plays a dominant role in astronomy, one may consider it as an additional reference point for the TSI. At the March equinox of 2010, the Earth’s heliocentric distance was about $r=\text{0}\text{.99594656}\text{AU}$ and, hence, $TSI\cong 1372\text{W}\cdot {\text{m}}^{-2}$ .

5. The Daily course of Insolation

In accord with Equation (10), the solar irradiation reaching the TOA at a given location during a certain period of a day, for instance, between sunrise and sunset is given by (e.g., [1-5, 9, 25, 27, 52])

$W_{S\downarrow }={\displaystyle \underset{t_r}{\overset{t_s}{\int }}\frac{S}{{\rho }^2}\left(\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h\right)\text{d}t}$ ,(60)

where $t_r$ and $t_s$ are the times of the location-dependent sunrise and sunset, respectively. As pointed out by Milankovitch [5], one must be aware of the discontinuity of radiation between sunset and subsequent sunrise, during which the considered element at the TOA is not irradiated. Thus, the integration between two arbitrary points of time, $t_1$ and $t_2$ , of a day must fulfil the condition $t_r\le t_1<t_2\le t_s$ . As the variation of the Earth’ heliocentric distance, $r$ , and the Sun’s geocentric declination, ${\delta }_S$ , during a day is negligible and may be replaced by their values at local solar noon, Equation (60) leads to

$W_{S\downarrow }=\frac{S}{{\rho }^2}\left(\sin \varphi \sin {\delta }_S{\displaystyle \underset{t_r}{\overset{t_s}{\int }}\text{d}t}+\cos \varphi \cos {\delta }_S{\displaystyle \underset{t_r}{\overset{t_s}{\int }}\cos h\text{d}t}\right)$ .(61)

Because

$\frac{\text{d}W}{\text{d}t}=\frac{\text{d}W}{\text{d}h}\frac{\text{d}h}{\text{d}t}$ ,(62)

we may insert the mean angular velocity of the Earth’s rotation, $\Omega =\text{d}h/\text{d}t=2\pi /d_a$ (e.g., [1, 3, 7-9]),

into Equation (61), where $d_a$ is the local apparent solar day, counted from the local solar noon to the following local solar noon. Thus, Equation (61) may be written as

$W_{S\downarrow }=\frac{S}{{\rho }^2}\left(\sin \varphi \sin {\delta }_S{\displaystyle \underset{-H}{\overset{H}{\int }}\frac{\text{d}h}{\Omega }}+\cos \varphi \cos {\delta }_S{\displaystyle \underset{-H}{\overset{H}{\int }}\cos h\frac{\text{d}h}{\Omega }}\right)$ ,(63)

i.e., the integration is only performed from sunrise ($-H$ ) to sunset ($H$ ), where $H$ is the half-day.

The apparent solar time, $t_a$ , varies during the Earth’s annual revolution around the Sun due to the eccentricity of the Earth’s orbit (see, e.g., Kepler’s second law), the obliquity of the ecliptic and small variations resulting from irregularities in the Earth’s rotation on its axis. Its difference from the mean solar time, $t$ , can be determined by the so-called equation of time, $E_t$ , by (e.g., [5, 8, 48])

$t_a=t+E_t$ .(64)

Thus, we have

$\frac{\text{d}{t}_a}{\text{d}t}=1+\frac{\text{d}{E}_t}{\text{d}t}$ .(65)

As $\left|\text{d}{E}_t/\text{d}t\right|\ll 1$ , this term can be neglected so that $\text{d}{t}_a\cong \text{d}t$ . Furthermore, $d_a$ varies with time during the Earth’s annual revolution around the Sun for the same reasons. The mean solar day, $d$ , differs from $d_a$ by

$d_a=d+E_t=d\left(1+\frac{E_t}{d}\right)$ .(66)

(a) (b)

Figure 17. Variation in the equation of time (a) through the year and (b) in dependence of the Sun’s geocentric declination, the so-called analemmic curve (in accord with [36, 37, 66, 73]).

As $\left|E_t/d\right|\ll 1$ (see Figure 17), we may use $d_a\cong d$ . Note that the mean solar day is about 3^m56^s longer than the sidereal day. Furthermore, $E_t$ is independent of latitude but varies over the year and depends on the Sun’s geocentric declination (see Figure 17).

The half-day can be deduced using Equation (5) by setting $\cos {\Theta }_0=0$ (e.g., [1, 3, 5, 9, 11, 15])

$H=\arccos \left(-\frac{\sin \varphi \sin {\delta }_S}{\cos \varphi \cos {\delta }_S}\right)=\arccos \left(-\tan \varphi \tan {\delta }_S\right)$ ,(67)

i.e., local sunrise and sunset mainly depend on $\varphi$ and ${\delta }_S$ . This formula, however, is only valid if the half-day fulfills the condition $H<\pi$ [3] ($\pi$ corresponds to 12 hours, i.e., 43200 seconds). The formula is invalid for those locations within the polar domes bounded by the polar circles. At these locations, the Sun does not set within 24 hours for which $\cos {\Theta }_0$ does not become zero [3]. Thus, the condition $H=\pi$ is that of the polar circle, i.e., $\varphi \cong -66˚33\text{'}42\text{'}\text{'}$ (Antarctic circle) and $\varphi \cong 66˚33\text{'}42\text{'}\text{'}$ (Arctic circle). For $\varphi =0˚$ and for the condition of the equinoxes, ${\delta }_S=0˚$ , one obtains $H=\pi /2$ , i.e., the diurnal arc corresponds to 12 hours. Figure 18 illustrates the latitudinal-annual distribution of daylight for 2010.

Because the Earth’s atmosphere causes a refraction of light, the following equation must be used to more accurately calculate sunrise and sunset for a given location at the Earth’s surface (e.g. [47, 48, 80]):

$H=\arccos \left(\frac{\sin \gamma -\sin \varphi \sin {\delta }_S}{\cos \varphi \cos {\delta }_S}\right)=\arccos \left(\frac{\sin \gamma }{\cos \varphi \cos {\delta }_S}-\tan \varphi \tan {\delta }_S\right)$ ,(68)

where ${\gamma }_S$ is the position of the center of the solar disk with reference to the horizon. Customarily, it is assumed that ${\gamma }_S=-34\text{'}-16\text{'}=-50\text{'}=-0.833˚$ , i.e., the center of the solar disk is located below the horizon. The value of −34’ applies to the atmospheric refraction of light with reference to the horizon. In the case of the Sun, the calculated times generally refer to the rise or set of the upper edge of the solar disk, so that −16’ has been included to consider the radius of the solar disk.

Figure 18. Latitudinal-annual distribution of the daylight hours for 2010. The dashed line is the Sun’s apparent geocentric declination.

Because the angular velocity $\Omega$ is a daily mean value, the solution of Equation (63) is given by [2, 3, 9, 10, 16]

$W_{S\downarrow }=\frac{Sd}{{\rho }^2\pi }\left(H\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\sin H\right)$ (69)

and with $\varphi =\pi /2-\theta$

$W_{S\downarrow }=\frac{Sd}{{\rho }^2\pi }\left(H\cos \theta \sin {\delta }_S+\sin \theta \cos {\delta }_S\sin H\right)$ (70)

Using Equation (7), the Sun’s geocentric declination may be replaced in these equations by the apparent ecliptic longitude of the Sun leading to either [8]

$W_{S\downarrow }=\frac{Sd}{{\rho }^2\pi }\left(H\sin \varphi \sin {\lambda }_S\sin \epsilon +\cos \varphi \sqrt{1-{\sin }^2\epsilon {\sin }^2{\lambda }_S}\sin H\right)$ (71)

or [35]

$W_{S\downarrow }=\frac{Sd}{{\rho }^2\pi }\left(H\sin \varphi \sin {\lambda }_S\sin \epsilon +\cos \varphi \cos \left(\arcsin \left(\sin {\lambda }_S\sin \epsilon \right)\right)\sin H\right)$ ,(72)

where the variation of ${\lambda }_S$ and $\epsilon$ during the day is negligible. Obviously, the condition of the equinoxes, ${\lambda }_S=0˚$ and ${\lambda }_S=180˚$ , leads to ${\delta }_S=0˚$ , $\sin H=1$ and, hence,

$W_{S\downarrow }=\frac{Sd}{{\rho }^2\pi }\cos \varphi$ .(73)

Thus, the amount of solar energy is proportional to the cosine of the latitude, where its maximum is related to the equator and its minimum related to both poles is zero [5].

If the Sun does not rise over the polar dome of the winter hemisphere, $W_{S\downarrow }=0$ . If the Sun does not

set over the polar dome of the summer hemisphere, one obtains under consideration of the condition of the polar circle, $H=\pi$ , [2-5]

$W_{S\downarrow }=\frac{Sd}{{\rho }^2}\sin \varphi \sin {\delta }_S=\frac{Sd}{{\rho }^2}\sin \varphi \sin \epsilon \sin {\lambda }_S$ .(74)

For the poles, $\sin \varphi =\pm 1$ and, hence, $\sin \varphi \sin {\delta }_S=\left|\sin {\delta }_S\right|$ . Since $\cos \varphi =0$ , Equation (5) reduces to

$\cos {\Theta }_0=\pm \sin {\delta }_S$ .(75)

Consequently, the annual variation of the insolation at the TOA over the poles only depends on the Sun’s geocentric declination. Sunrise and sunset at the poles can simply be deduced using the condition of the equinoxes ${\delta }_S=0˚$ or ${\lambda }_S=0˚$ and ${\lambda }_S=180˚$ , respectively. As illustrated in Figure 19, at the North Pole, sunrise takes place at the time of the March equinox and sunset at the time of the subsequent September equinox. At the South Pole, sunrise takes place at the time of this September equinox and the sunset at time of the subsequent March equinox. In the case of the poles, Equation (6) leads to Equation (75). Consequently, the maximum value, $F_{S\downarrow ,\max }$ , and the instantaneous value, $F_{S\downarrow }$ , of the solar irradiance are identical (see Figure 19). In 1906, Wilhelm von Bezold [81] already pointed out that the maximum of daily solar irradiation, as has long been known, falls on the pole of the summer hemisphere and the absolute maximum on the South Pole.

(a) (b)

Figure 19. Comparison of the maximum value, $F_{S\downarrow ,\max }$ , and the instantaneous value, $F_{S\downarrow }$ , of the solar irradiance for (a) the South Pole and (b) the North Pole.

Replacing $\rho$ in Equation (69) by $\rho =\left(1-e^2\right)/\left(1+e\cos \upsilon \right)$ (see Equation (51)) yields

$W_{S\downarrow }=\frac{{\left(1+e\cos \upsilon \right)}^2}{{\left(1-e^2\right)}^2}\frac{Sd}{\pi }\left(H\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\sin H\right)$ ,(76)

where the term

$\frac{{\left(1+e\cos \upsilon \right)}^2}{{\left(1-e^2\right)}^2}=\frac{1+2e\cos \upsilon +e^2{\cos }^2\upsilon }{1-2e^2+e^4}\cong 1+2e\cos \upsilon +e^2{\cos }^2\upsilon$

may be approximated by a series expansion [5, 8]. Note that $e^4\ll \left|2e^2\right|\ll 1$ . Equations (70) - (74) and the following Equations (77) - (79) may be treated analogously.

In accord with Equation (14), we obtain for the location-dependent daily mean solar irradiance (expressed either in J∙m⁻²∙s⁻¹ or in W∙m⁻² [8, 15]),

$\overline{W_{S\downarrow }}=\frac{W_{S\downarrow }}{d}=\frac{S}{{\rho }^2\pi }\left(H\cos \theta \sin {\delta }_S+\sin \theta \cos {\delta }_S\sin H\right)$.(77)

Equations (69) and (77) document that only the insolation during the diurnal arc is to be considered. These equations are only valid if the half-day fulfills the condition $H<\pi$ . Obviously, Spitaler’s argument that “the definition of the mean irradiance of a parallel is ambiguous and that contradictions arise once one goes beyond a day” related to Hopfner’s findings [59] is invalid.

As in the case of Equation (73), the condition of the equinoxes leads to

$\overline{W_{S\downarrow }}=\frac{S}{{\rho }^2\pi }\cos \varphi$.(78)

If the Sun does not rise over the polar dome of the winter hemisphere, $\overline{W_{S\downarrow }}=0$ . On the contrary, if the Sun does not set over the polar dome of the summer hemisphere, one obtains [15]

$\overline{W_{S\downarrow }}=\frac{S}{{\rho }^2}\sin \varphi \sin {\delta }_S=\frac{S}{{\rho }^2}\sin \varphi \sin \epsilon \sin {\lambda }_S$.(79)

If the condition $H<\pi$ is fulfilled, the zonal average of the location-dependent daily mean solar irradiance $\overline{W_{S\downarrow }}$ can be calculated using Equation (15). One obtains

$\widehat{\overline{W_{S\downarrow }}}=\frac{1}{2\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}\overline{W_{S\downarrow }}\text{d}\phi }=\overline{W_{S\downarrow }}$.(80)

The global average of the solar irradiance at the TOA for one day is therefore

$\langle \overline{W_{S\downarrow }}\rangle =\frac{1}{2}{\displaystyle \underset{0}{\overset{\pi }{\int }}\widehat{\overline{W_{S\downarrow }}}\left(\theta \right)\sin \theta \text{d}\theta }=\frac{1}{2}{\displaystyle \underset{0}{\overset{\pi }{\int }}\overline{W_{S\downarrow }}\left(\theta \right)\sin \theta \text{d}\theta }$,(81)

finally expressed by

$\langle \overline{W_{S\downarrow }}\rangle =\frac{1}{2}\frac{S}{{\rho }^2\pi }{\displaystyle \underset{0}{\overset{\pi }{\int }}\left(H\cos \theta \sin {\delta }_S+\sin \theta \cos {\delta }_S\sin H\right)\sin \theta \text{d}\theta }$.(82)

This formula can only be solved in the case of the equinoxes, i.e., for ${\delta }_S=0$ . As mentioned before, the half-day is then $H=\pi /2$ and, hence, $\sin H=1$ . Thus, for the equinoxes, one finally obtains

$\langle \overline{W_{S\downarrow }}\rangle =\frac{1}{2}\frac{S}{{\rho }_e^2\pi }{\displaystyle \underset{0}{\overset{\pi }{\int }}{\sin }^2\theta \text{d}\theta }=\frac{S}{4{\rho }_e^2}$.(83)

Here, ${\rho }_e=r_e/r_0$ , where $r_e$ is the Earth’ heliocentric distance at the time of the respective equinox. For the equinoxes, this result confirms that the global average of the daily mean solar irradiance corresponds to a quarter of the solar constant. This result known since the second half of the 19^th century (Meech [1] and Wiener [2, 3]), later confirmed by Milankovitch [4, 5], List [6], Fortak [13], Raschke et al. [14], Liou [9], Fu [10], Berger and Yin [15] and many others has been disputed by the EIKE experts. For the March

and September equinoxes of 2009, we obtained $r_e=0.99596760\text{AU}$ and $\langle \overline{W_{S\downarrow }}\rangle =343.0\text{W}\cdot {\text{m}}^{-2}$ and $r_e=1.00354280\text{AU}$ and $\langle \overline{W_{S\downarrow }}\rangle =337.9\text{W}\cdot {\text{m}}^{-2}$, respectively. For the March equinox 2010, we obtained $r_e=0.99594656\text{AU}$ and $\langle \overline{W_{S\downarrow }}\rangle =343.0\text{W}\cdot {\text{m}}^{-2}$. For the entire tropical year 2009/2010 and the year 2010, we obtained $\langle \overline{W_{S\downarrow }}\rangle =340.2\text{W}\cdot {\text{m}}^{-2}$. Using the data of the daily averaged solar irradiation at the TOA listed by Kopp in his Table 1 [16] yields for 2023 $\langle \overline{W_{S\downarrow }}\rangle =340.1\text{W}\cdot {\text{m}}^{-2}$. This means that Meech’s and Wiener’s findings are also confirmed by Kopp’s and our results, i.e., by results that were derived with other methods than used by Meech [1] and Wiener [2, 3].

6. The Insolation for any number of days in the year

Results from our numerical predictions of the solar irradiance at the TOA performed for 2010 are shown for various parallels of the NH in Figure 20 and those of the SH in Figure 21 with the equator arbitrarily assigned to the NH. The quantities $r$ and ${\delta }_S$ were again computed using the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. These figures illustrate the diurnal and seasonal variations of the solar irradiance due to the Earth’s rotation and the variation of the Earth’s heliocentric distance and the Sun’s geocentric declination. The height and width of the curves of insolation depend on the parallel of latitude, but these curves vary during the annual course. At a given parallel of latitude, the maximum amount of the solar irradiance depends on both $r$ and ${\delta }_S$ (see Equations (1), (4) and (6)), and the duration of isolation during the diurnal arc, $2H$ , and, hence, on $\varphi$ and ${\delta }_S$ (see Equation (67)). Furthermore, at 75˚N and 85˚N the solar irradiance exhibits a diurnal variation during several days before and after the June solstice, despite the Sun does not set. At the North Pole, the Sun does not set during this period, but the solar irradiance varies weakly. The same is true at 75˚S and 85˚S and at the South Pole during several days before and after the December solstice.

Daily cycles of the location-dependent solar irradiance at the TOA, as illustrated by Figure 20 and Figure21 for several parallels of latitude of both hemispheres, must be calculated using Equation (8), where the discontinuity of insolation between sunset and subsequent sunrise must be considered. Therefore, it makes sense to use this equation in numerical simulations generally. Because the cosine of the hour angle is negative for $\pi /2<h<3\pi /2$ , the insolation must be set to zero if

$\cos {\Theta }_0=\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h<0$ (84)

or (see Equation (11)) [8]

$\cos {\Theta }_0=\sin \varphi \sin \epsilon \sin {\lambda }_S+\cos \varphi \sqrt{1-{\sin }^2\epsilon {\sin }^2{\lambda }_S}\cos h<0$ (85)

or (see Equation (13)) [35]

$\cos {\Theta }_0=\sin \varphi \sin \epsilon \sin {\lambda }_S+\cos \varphi \cos \left(\arcsin \left(\sin \epsilon \sin {\lambda }_S\right)\right)\cos h<0$ .(86)

These are the conditions of the local night-time.

As illustrated by Figure 20 and Figure 21, the derivative $\text{d}{W}_{S\downarrow }/\text{d}t$ given by Equation (10) is also affected by the discontinuity caused by the local night-time. It is therefore indispensable to take this discontinuity into account if an arbitrary period of days during the year is considered as illustrated in Figure 22. This figure shows the diurnal variation of the solar irradiance and the daily mean solar irradiance at the TOA for the latitude of 45˚N during the period of six days around the 2010 March equinox.

As suggested by Milankovitch [5], the consideration of the insolation discontinuity can be done as follows: The integration between the beginning of the first day, $t_b=0^{\text{h}}\text{UT}1$ , and the end of the $i^{th}$ day ($i=1,2,\cdots ,n$ ), $t_e={24}^{\text{h}}\text{UT}1$ , can be described by

$W_{S\downarrow }=S{\displaystyle \underset{t_b}{\overset{t_e}{\int }}\frac{1}{{\rho }^2}\left(\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h\right)\text{d}t}$ .(87)

Because of the discontinuity of insolation we must consider the time interval of each day between local sunrise, $t_{r,i}$ , andlocal sunset, $t_{s,i}$ , $i=1,2,\cdots ,n$ , where $n\le 365$ ($n\le 366$ in leap years). The solar irradiance at the TOA for a certain parallel of latitude is therefore given by

(a) (b)

(c) (d)

(e) (f)

Figure 20. Diurnal variation of the solar irradiance, $F_{S\downarrow }$ , at the TOA in 2010 for different parallels of latitude of the NH around (a) the perihelion, (b) the March equinox, (c) the June solstice, (d) the aphelion, (e) the September equinox, and (f) the December solstice. The numerical predictions are based on the Equation (8), where a solar constant of $S=1361\text{W}\cdot {\text{m}}^{-2}$ was used. The Earth’ heliocentric distance and the Sun’s apparent geocentric declination were predicted using the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19].

$\begin{array}{l}{W}_{S\downarrow }=S({\displaystyle \underset{t_{r,1}}{\overset{t_{s,1}}{\int }}\frac{1}{{\rho }_1^2}\left(\sin \varphi \sin {\delta }_{S,1}+\cos \varphi \cos {\delta }_{S,1}\cos h\right)\text{d}t}\\ +{\displaystyle \underset{t_{r,2}}{\overset{t_{s,2}}{\int }}\frac{1}{{\rho }_2^2}\left(\sin \varphi \sin {\delta }_{S,2}+\cos \varphi \cos {\delta }_{S,2}\cos h\right)\text{d}t}\\ +\cdots \\ +{\displaystyle \underset{t_{r,n}}{\overset{t_{s,n}}{\int }}\frac{1}{{\rho }_n^2}\left(\sin \varphi \sin {\delta }_{S,n}+\cos \varphi \cos {\delta }_{S,n}\cos h\right)\text{d}t}).\end{array}$ (88)

As mentioned before, the variation of the Earth’s heliocentric distance, $r$ , and the Sun’s geocentric declination, ${\delta }_S$ , during a day is negligible. Thus, the values of ${\rho }_i$ , ${\delta }_{S,i}$ and $H_i$ may be replaced by their values at local solar noon of the $i^{th}$ day, where ${\rho }_i$ and ${\delta }_{S,i}$ (or ${\lambda }_{S.i}$ ) can be taken, for instance, from the Astronomical Almanac (Section C Sun).

(a) (b)

(c) (d)

(e) (f)

Figure 21. As in Figure 20, but for the SH.

Figure 22. Diurnal variation of the solar irradiance, $F_{S\downarrow }$ , at the TOA and the daily mean solar irradiance for the latitude of 45˚N during the period of six days around the 2010 March equinox.

In accord with Equation (61), we may write

$\begin{array}{l}{W}_{S\downarrow }=S\left(\frac{1}{{\rho }_1^2}(\sin \varphi \sin {\delta }_{S,1}{\displaystyle \underset{-H_1}{\overset{H_1}{\int }}\frac{\text{d}h}{\Omega }}+\cos \varphi \cos {\delta }_{S,1}{\displaystyle \underset{-H_1}{\overset{H_1}{\int }}\cos h}\frac{\text{d}h}{\Omega }\right)\\ +\frac{1}{{\rho }_2^2}\left(\sin \varphi \sin {\delta }_{S,2}{\displaystyle \underset{-H_2}{\overset{H_2}{\int }}\frac{\text{d}h}{\Omega }}+\cos \varphi \cos {\delta }_{S,2}{\displaystyle \underset{-H_2}{\overset{H_2}{\int }}\cos h\frac{\text{d}h}{\Omega }}\right)\\ +\cdots \\ +\frac{1}{{\rho }_n^2}\left(\sin \varphi \sin {\delta }_{S,n}{\displaystyle \underset{-H_n}{\overset{H_n}{\int }}\frac{\text{d}h}{\Omega }}+\cos \varphi \cos {\delta }_{S,n}{\displaystyle \underset{-H_n}{\overset{H_n}{\int }}\cos h\frac{\text{d}h}{\Omega }}\right)),\end{array}$ (89)

where $\Omega ={\Omega }_1={\Omega }_2=\cdots ={\Omega }_n=2\pi /d_a$ is assumed. Following Equation (69), the solution is given by

$W_{S\downarrow }=\frac{Sd}{\pi }{\displaystyle \sum _{i=1}^n\frac{1}{{\rho }_i^2}\left(H_i\sin \varphi \sin {\delta }_{S,i}+\cos \varphi \cos {\delta }_{S,i}\sin H_i\right)}$ .(90)

Figure 23. The daily mean irradiance at the TOA for different parallels of latitude determined by diurnal averaging (see Equation (14)) using 144 values per day, and Equation (91) for the 365 days of 2010, respectively.

Figure 24. Latitudinal-annual distribution of the daily mean irradiance at the TOA for the tropical year 2009/2010 starting with $TDB=2454910.5$ (March 20, 2009, 00:00 UT1), where a solar constant of $S=1361\text{W}\cdot {\text{m}}^{-2}$ was used. The Earth’s heliocentric distance and the Sun’s apparent geocentric declination are based on the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. The dashed line indicates the Sun’s apparent geocentric declination and the dotted lines the Equator, the Tropic of Cancer and the Tropic of Capricorn and the Arctic Circle and the Antarctic Circle.

In accord with the definition of the daily average (Equation (14)), we finally obtain

$\overline{W_{S\downarrow }}=\frac{S}{\pi }{\displaystyle \sum _{i=1}^n\frac{1}{{\rho }_i^2}\left(H_i\sin \varphi \sin {\delta }_{S,i}+\cos \varphi \cos {\delta }_{S,i}\sin H_i\right)}$.(91)

As illustrated in Figure 23, the 2010 daily mean solar irradiance at the TOA for different parallels of latitude obtained by diurnal averaging based on Equation (14)) using 144 data per day, and the Equation

(91) for the 365 days of 2010 are nearly identical, where $\overline{W_{S\downarrow }}\left(\varphi =5˚\text{N}\right)=414.2\text{W}\cdot {\text{m}}^{-2}$, $\overline{W_{S\downarrow }}\left(\varphi =25˚\text{N}\right)=380.4\text{W}\cdot {\text{m}}^{-2}$, $\overline{W_{S\downarrow }}\left(\varphi =45˚\text{N}\right)=307.1\text{W}\cdot {\text{m}}^{-2}$ and $\overline{W_{S\downarrow }}\left(\varphi =65˚\text{N}\right)=213.8\text{W}\cdot {\text{m}}^{-2}$. If the

diurnal variation of the solar irradiance is not to be accounted for explicitly, Equation (78) can be used with a sufficient degree of accuracy to determine the daily mean values.

As already done before (see Figure 3), the variation of the irradiation at the TOA may be illustrated by the latitudinal-time-dependent distribution of the tropical year 2009/2010 (Figure 24). In accord with Smulsky [35], we also derived the insolation for the summer half-year (astronomic seasons of Spring and Summer) 2009 and the winter half-year (astronomic seasons of Fall and Winter) 2009/2010, and the halved insolation for the tropical year 2009/2010 (Figure 25).

7. The annual Course of Insolation

Figure 26 illustrates the daily mean solar irradiance at the TOA for different latitudes of the SH and NH, where the equator is arbitrarily assigned to the NH. Based on the Equation (8), the calculations were performed for 2010 as well. Figure 27 that is based on the data listed by Kopp in his Table 1 [16] shows the results for 2023. Obviously, the differences between both figures are negligible. This means that the interannual variation of annual insolation is marginal on the decadal time-scale.

As Figure 26 and Figure 27 illustrate, the maximum of the daily mean solar irradiance during the summer of the SH is up to 6.7 % greater than that during the summer of the NH, due to the variation in the

Figure 25. Smulsky-diagram [35] of the insolation for the 2009 summer half-year, the 2009/2010 winter half-year and the halved insolation for the tropical year 2009/2010. The open circles illustrate Milankovitch’s values for the several latitudes ranging from $\varphi =0˚$ to $\varphi =90˚$ in steps of $\Delta \varphi =5˚$ listed in his Table VII. Milankovich’s results were re-scaled to the current solar constant of $S=1361\text{W}\cdot {\text{m}}^{-2}$ .

Earth’s heliocentric distance. However, as mentioned before, the time span between the March equinox and the subsequent September equinox is about 7.6 days longer than that between this September equinox and the subsequent March equinox. The solar irradiation due to this temporal deficit of the SH is largely compensated by the “excess” solar radiation, so that for the same latitudes both hemispheres receive almost the same annual irradiation despite the differences in irradiation in the same seasons [82] (see also Table 2).

(a) (b)

Figure 26. The daily mean irradiance, $\overline{W_{S\downarrow }}$ , at the TOA for different parallels of latitude of (a) the SH and (b) the NH for 2010 predicted with the Equation (8) using a solar constant of $S=1361\text{W}\cdot {\text{m}}^{-2}$ . The Earth’s heliocentric distance and the Sun’s apparent geocentric declination are based on the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. The daily mean values were computed with Equation (1) using 144 data per day.

(a) (b)

Figure 27. As in Figure 26, but using Kopp’s data for 2023 listed in his Table 1 [16].

Because of the annual course of the Sun’s geocentric declination, the daily mean solar irradiance has a pronounced, annual course, too [82]. Due to the Sun’s two zenith positions over the course of the year, the amount of daily mean solar irradiance at the equator reaches two local maxima and two local minima (see also Figure 28(b)): The former occurs about ten days before the March equinox and about thirteen days after the September equinox; the latter is near of the solstices. The absolute minimum of 384.8 W∙m⁻² occurs three days after the June solstice and the absolute maximum of 438.1 W∙m⁻² occurs ten days before the March equinox. Because the average is about 415.5 W∙m⁻², the variation with respect to this average is relatively small. Also the area around the equator, up to about 12˚N and 12˚S latitude, shows a double period of daily mean solar radiation, but the two maxima move closer together [82]. Beyond these latitudes, poleward, only one maximum and one minimum occur, i.e., there is only one period. At 20˚N and 20˚S latitude, the difference between the summer maximum and the winter minimum is already three to four times larger than the difference at the equator. The difference grows rapidly with increasing latitudes [82]. At the south and north pole, the respective summer maximum of daily solar irradiance exceeds that of the respective winter minimum (0 W/m²) by 559.4 W∙m⁻² and 524.2 W∙m⁻², respectively. This means that around the time of a hemisphere’s summer solstice, the daily mean solar radiation at the pole exceeds the daily mean solar radiation occurring at the Equator. Furthermore, at December solstice, the daily mean solar radiation at the South Pole is almost 7% higher than that at the North Pole around June solstice.

As already mentioned, the annual variation of the irradiation at the TOA may be illustrated by the latitudinal-time-dependent distribution of either (a) the daily insolation (expressed in MJ∙m⁻²) or (b) the daily mean irradiance (expressed in W∙m⁻²). Results for 2010 are illustrated in Figure 29.

8. seasonal and interdecadal Variations of Insolation during the past 400 years

As mentioned before, the JPL planetary und lunar ephemeris DE440 covers the years 1550 - 2650 [17]. However, the reconstruction of the TSI at 1 AU (called the solar constant) illustrated in Figure 16 only goes back to 1610. Thus, for our study, we calculated the seasonal variations in the daily mean solar irradiance at the TOA for 1610, 1650, 1710, 1750, 1810, 1850, 1910 and 1950 with respect to those of 2010 shown in Figure 26. The respective values of the solar constant were taken from the data on which Figure 16 is based. The differences for several parallels of latitude of the SH and NH are illustrated in Figure 30. Compared to the daily mean solar irradiance of 2010 in which the solar activity was low, there are differences in the seasonal variations of the daily mean solar irradiance in the range of a few Watts per square meters for all years considered here.

(a) (b)

Figure 28. As in Figure 24, but only for latitudes from $\varphi =-15˚$ to $\varphi =15˚$ .

(a) (b)

Figure 29. Latitudinal-annual distribution of (a) the daily insolation (expressed in MJ∙m⁻²) and (b) the daily mean irradiance (expressed in W∙m⁻²) at the TOA for 2010 starting with TDB = 2455197.5 (January 1, 2010, 00:00 UT1). The Earth’s heliocentric distance and the Sun’s apparent geocentric declination were predicted using a solar constant $S=1361\text{W}\cdot {\text{m}}^{-2}$ and the data provided by the JPL ephemeris DE440 [17] and the NOVAS F3.1 subroutines [18, 19]. The dashed lines indicate the Sun’s apparent geocentric declination, and the dotted lines the Equator, the Tropic of Cancer and the Tropic of Capricorn and the Arctic Circle and the Antarctic Circle.

Even though the reconstruction of the TSI does not vary more than 2.2 W∙m⁻², we found, particularly in the polar regions of both hemispheres, that there are larger seasonal differences in daily mean solar radiation, mainly with different signs, resulting in a butterfly-like distribution of these differences across latitudes and seasons. The largest deviations in the daily mean solar irradiance with respect to 2010 occur in 1910 (SH: from −7.3 W∙m⁻² to 7.4 W∙m⁻²; NH: from −7.4 W∙m⁻² to 7.3 W∙m⁻²) followed by 1810 (SH: from −5.0 W∙m⁻² to 5.5 W∙m⁻²; NH: from −5.8 W∙m⁻² to 5.2 W∙m⁻²), 1950 (SH: from −4.4 W∙m⁻² to 4.4 W∙m⁻²; NH: from −4.5 W∙m⁻² to 4.3 W∙m⁻²) and 1710 (SH: from −2.7 W∙m⁻² to 3.5 W∙m⁻²; NH: from −4.6 W∙m⁻² to 3.4 W∙m⁻²). In the case of these years, we may conclude that the daily mean solar irradiance is

(a) (b)

(c) (d)

(e) (f)

(g) (h)

(i) (j)

(k) (l)

(m) (n)

(o) (p)

Figure 30. The differences of the daily mean solar irradiance for the years 1610 (a, b), 1650 (c, d), 1710 (e, f), 1750 (g, h), 1810 (i, j), 1850 (k, l), 1910 (m, n) and 1950 (o, p) with respect to the year 2010 for several parallels of latitude of the SH (a, c, e, g, i, k, m, o) and NH (b, d, f, h, j, l, n, p). The equator is arbitrarily assigned to the NH.

most sensitive to secular changes in the astronomic elements from the polar circles polewards and least sensitive around the equator. The smallest differences in the daily mean solar irradiance with respect to 2010 occur in 1750 (SH: from −1.9 W∙m⁻² to 1.1 W∙m⁻²; NH: from −2.0 W∙m⁻² to 1.2 W∙m⁻²). In this case, the largest positive differences occur around the equator.

To interpret our results, the difference $\Delta r=r_{2010}-r_{\chi }$ in the Earth’s heliocentric distance and the difference $\Delta {\delta }_S={\delta }_{S,2010}-{\delta }_{S,\chi }$ in the Sun’s apparent geocentric declination are shown in Figure 31. Here, $\chi$ marks the year considered. Because the $\Delta {\delta }_S$ -curves for some years look very different for some years than for others, we compared our results with those provided by the JPL Horizons on-line solar system (https://ssd.jpl.nasa.gov/horizons/) that are based on the JPL planetary and lunar ephemeris DE441. We found that the results are identical. The results of our comparison are illustrated in Figure A1 (see Appendix). In addition, we compared the $\Delta {\delta }_S$ -curves with those obtained from the results provided by the older JPL planetary and lunar ephemeris E430. Also in this case, the results are identical.

(a) (b)

Figure 31. Differences with respect to 2010: (a) $\Delta r=r_{2010}-r_{\chi }$ in the Earth’s heliocentric distance, and (b) $\Delta {\delta }_S={\delta }_{S,2010}-{\delta }_{S,\chi }$ in the Sun’s apparent geocentric declination. Here, $\chi$ stands for 1610, 1650, 1710, 1750, 1810, 1850, 1910 and 1950.

The differences in the Sun’s apparent geocentric declination with respect to 2010 were the largest for 1910, followed by 1810 and 1950 (Figure 31(b)). At four times within the annual course, the Sun’s apparent geocentric declination is nearly the same for 1750 and 2010. On the contrary, 1950,1910, 1850, 1810, 1710, 1650, and 1610 had zero differences in the Sun’s geocentric declination with respect to 2010, only twice in the annual course. However, the timing of zero differences differed among these years. It seems that the multiple occurrences of the same Sun’s geocentric declination for 2010 and 1750 explain the quite different behavior of the differences in daily mean solar irradiance obtained for 1750 compared to the other cases shown in Figure 30. Like 1750, the 1610 and 1650 differences in Sun’s apparent geocentric declination with respect to 2010 showed a secondary maximum and a secondary minimum (Figure 31(b)). These years’ differences in daily mean solar irradiance with respect to 2010 display a quite different distribution within the annual course compared to years with just one minimum and one maximum (Figure 30 and Figure 31(b)). Differences in mean daily solar irradiance with respect to 2010 were the largest in years, in which the difference in the Sun’s apparent geocentric declination with respect to 2010 had only one minimum and one maximum.

Differences in the Earth’s heliocentric distance with respect to 2010 are negative between perihelion and aphelion and positive from the aphelion to the perihelion (Figure 31(a)). Except for the solstices, the Earth’s heliocentric distance of 1610 differed the most from that of 2010 and was the largest around the equinoxes. Differences with respect to 2010 increased gradually backward in time.

Furthermore, on the multi-decadal and century-scales, the smooth changes of astronomical elements yield the largest differences in daily mean solar irradiance from the polar circles polewards. In both hemispheres, these spatial differences start increasing around the time of the perihelion at 70˚ followed by the start of the increase about a month later at 80˚. The differences in daily mean solar irradiation are opposite in sign on the NH and SH (Figure 30). The multi-decadal and inter-centurial variation is smallest for the Equator. These findings mean that the polar regions are the most sensitive to changes in the Earth’s heliocentric distance, the Sun’s geocentric declination, and the solar constant. Furthermore, the time of zero differences in daily mean solar irradiance along all parallels of latitudes can differ between the NH and SH. In some years, such a time of zero difference doesn’t occur (Figure 30).

9. Summary and CONCLUSIONS

In our paper, we assessed the prediction of insolation at the TOA for different time scales ranging from the diurnal course to the annual course and arbitrary periods of days in between. We also reviewed the historical work on solar irradiance performed by Meech [1] and Wiener [2, 3] and assessed the accuracy of their findings using our results for the year 1853 and the tropical years 1874/1875 and 2009/2010. In addition, we investigated the variation of the daily mean irradiance at the TOA under smooth changes of astronomical elements like the Earth’s heliocentric distance and the Sun’s apparent geocentric declination during the past four centuries by predicting the daily mean solar irradiance over the annual course of the years 1610, 1650, 1710, 1750, 1810, 1850, 1910 and 1950 and comparing the results with those obtained for 2010.

Because the heliocentric distance of the Earth, r, and the apparent geocentric declination of the Sun, ${\delta }_S$ , are affected by the variation of the Earth’s orbit around the Sun due to smooth changes in the astronomical conditions caused by precession of the equator and the ecliptic, nutation, changes in the obliquity (and long-term changes in eccentricity), trustworthy astronomical calculation methods are indispensable for long-term studies. Therefore, our predictions for all years considered in our study are based on the ICRS geocentric rectangular coordinates and their rates provided by the JPL ephemeris DE440 [17] and the adjustment of these data to the equator and the equinox of date using the subroutines for aberration, frame bias, precession and nutation of the Naval Observatory Vector Astrometry Software (NOVAS), Version F3.1 [18, 19]. Our predictions were performed for 90˚S to 90˚N latitudes using an increment of 5˚, where a 10-minute time step was used. For 2009 and 2010, the value of the TSI at 1 AU of $S=1361\text{W}\cdot {\text{m}}^{-2}$ was chosen. For all other years considered in our study, we took the data from the reconstruction of the TSI for the past four centuries that are based on the “Community-Consensus TSI Composite” [23] and SATIRE-T model [24] (with modifications to fix spurious values prior to 1650).

We evaluated the historical findings on the solar irradiance by Meech [1] and Wiener [2, 3] using our predictions for the year 1853 and the tropical years 1874/1875 and 2009/2010. Our results mainly confirm the results of Meech and Wiener. Spitaler’s criticism “the definition of the mean irradiance of a parallel is ambiguous and that contradictions arise once one goes beyond a day” related to Hopfner’s findings is, therefore, invalid. This means that Milankovitch’s comment to Hopfner’s findings is, indeed, correct.

Based on our predictions for the entire tropical year 2009/2010 and the year 2010, we obtained for the

global average of the solar irradiance at the TOA $\langle \overline{W_{S\downarrow }}\rangle =340.2\text{W}\cdot {\text{m}}^{-2}$. Using the data of the daily averaged solar irradiation at the TOA listed by Kopp in his Table 1 [16] yields for 2023 $\langle \overline{W_{S\downarrow }}\rangle =340.1\text{W}\cdot {\text{m}}^{-2}$. Thus,

Meech’s and Wiener’s findings that the global average of the daily mean solar irradiance is one quarter of the solar constant are also confirmed by Kopp’s and our results. In addition, we also confirmed the results of Kopp [16] for 2023.

Our results for all years considered in our study also demonstrate that the method for calculating solar irradiation propagated by the experts of the German association EIKE is based on pure ignorance. Apparently, the scientific literature on solar irradiation at the TOA is completely rejected by these EIKE experts.

Even though the reconstruction of the TSI does not vary more than 2.2 W∙m⁻², we found, particularly in the polar regions of both hemispheres, that there are notable seasonal differences in daily mean solar radiation, mainly with different signs, resulting in a butterfly-like distribution of these differences across latitudes and seasons. The largest deviations in the daily mean solar irradiance with respect to 2010 occur in 1910 (SH: from −7.3 W∙m⁻² to 7.4 W∙m⁻²; NH: from −7.4 W∙m⁻² to 7.3 W∙m⁻²) followed by 1810 (SH: from −5.0 W∙m⁻² to 5.5 W∙m⁻²; NH: from −5.8 W∙m⁻² to 5.2 W∙m⁻²), 1950 (SH: from −4.4 W∙m⁻² to 4.4 W∙m⁻²; NH: from −4.5 W∙m⁻² to 4.3 W∙m⁻²) and 1710 (SH: from −2.7 W∙m⁻² to 3.5 W∙m⁻²; NH: from −4.6 W∙m⁻² to 3.4 W∙m⁻²). In the case of these years, we may conclude that the daily mean solar irradiance is most sensitive to smooth changes in the astronomic elements from the polar circles polewards and least sensitive around the equator.

The smallest differences in the daily mean solar irradiance with respect to 2010 occur in 1750 (SH: from −1.9 W∙m⁻² to 1.1 W∙m⁻²; NH: from −2.0 W∙m⁻² to 1.2 W∙m⁻²). In this case, the largest positive differences occur around the equator.

Our results for the years considered in our study reveal a remarkable nonlinear relationship between the annual course of daily mean solar irradiance and the combined effects of Earth’s heliocentric distance and the apparent geocentric declination of the Sun. Our results for 1610 exhibit the largest positive and negative deviations in the Earth’s heliocentric distance with respect to 2010. These deviations would cause differences in the daily mean solar irradiance with respect to 2010 that range from −4.8 W∙m⁻² to 6.1 W∙m⁻². However, these differences only vary from −2.3 W∙m⁻² to 2.0 W∙m⁻² for the SH and from −3.5 W∙m⁻² to 2.3 W∙m⁻² for the NH.

As mentioned before, the largest deviations in the daily mean solar irradiance with respect to 2010 occur in 1910. In this case, the Sun’s apparent geocentric declination shows the largest positive and negative deviations and the second smallest difference in the Earth’s heliocentric distance, both with respect to 2010. Thus, we may conclude that over the annual course, (1) the Sun’s apparent geocentric declination plays a more important role for deviations in daily mean solar irradiation than the Earth’s heliocentric distance; and (2) the impacts of the Earth’s heliocentric distance and the Sun’s geocentric declination are non-linear. Given the magnitude of differences in solar irradiance caused by secular changes of these astronomical elements, they must be considered on the multi-decadal to multi-centennial scales.

Our results demonstrate that the polar regions could be affected by differences in the daily mean solar irradiance due to smooth changes in the astronomic elements. However, since the zonal averages must be weighted either by $\sin \theta$ or $\cos \varphi$ in global averaging (see Equation (16)), differences in the daily mean solar irradiance due to smooth changes in the astronomic elements play an insignificant role on the global scale.

Acknowledgements

We thank the JPL team around Drs. James G. Williams and William M. Folkner for making the planetary and lunar ephemeris DE 440 available and the NOVAS team of the US Naval Observatory around Dr. George Kaplan for making the NOVAS F3.1 software available. Furthermore, we thank Drs. George Kaplan, US Naval Observatory, and Jon D. Giorgini, JPL Horizons, for helpful comments and suggestions.

Appendix

A1. A comparison between DE440 and DE441 results

(a) (b)

(c) (d)

(e) (f)

(g) (h)

(i) (j)

(k) (l)

(m) (n)

(o) (p)

Figure A1. Comparison of our results for the differences for $\Delta r=r_{2010}-r_{\chi }$ and $\Delta {\delta }_S={\delta }_{S,2010}-{\delta }_{S,\chi }$ with those provided by the JPL Horizons on-line solar system (https://ssd.jpl.nasa.gov/horizons/) that are based on the JPL ephemeris DE441. Here, $\chi$ stands for 1610 (a, b), 1650 (c, d), 1710 (e, f), 1750 (g, h), 1810 (i, j), 1850 (k, l), 1910 (m, n) and 1950 (o, p).

A2. Uncertainty Analysis

As requested by one of the reviewers, we inserted an uncertainty analysis. In accord with Mölders et al. [83] and Iorio [84], we consider principles of Gaussian error propagation (GEP) for an uncertainty analysis. In doing so, the equation to predict a quantity $\psi$ , here the solar irradiance at the TOA given by Equation (8), is derived for the astronomical quantities $r$ , $S$ and ${\delta }_S$ . The standard deviation of the predicted quantity can be calculated from the individual derivations ($\partial \psi /\partial {\chi }_i$ ) and the chosen uncertainty ${\sigma }_{{\chi }_i}$ of the astronomical quantities ${\chi }_i$ ($r$ , $S$ and ${\delta }_S$ ) by

${\sigma }_{\psi }=\sqrt{{\displaystyle \sum _{i=1}^n{\left(\frac{\partial \psi }{\partial {\chi }_i}\right)}^2{\sigma }_{{\chi }_i}^2}}$ ,(92)

where $n$ is the number of astronomical quantities considered. Thus, in accord with Equation (8), we may write

${\sigma }_{F_S}={\left({\left(\frac{\partial F_{S\downarrow }}{\partial r}\right)}^2{\sigma }_r^2+{\left(\frac{\partial F_{S\downarrow }}{\partial S}\right)}^2{\sigma }_S^2+{\left(\frac{\partial F_{S\downarrow }}{\partial {\delta }_S}\right)}^2{\sigma }_{{\delta }_S}^2\right)}^{1/2}$ ,(93)

The three derivations read

$\frac{\partial F_{S\downarrow }}{\partial r}=-\frac{2}{r}{F}_{S\downarrow }$ ,(94)

$\frac{\partial F_{S\downarrow }}{\partial S}=\frac{F_{S\downarrow }}{S}$ ,(95)

and

$\frac{\partial F_{S\downarrow }}{\partial {\delta }_S}=F_{S\downarrow }\frac{\sin \varphi \cos {\delta }_S-\cos \varphi \sin {\delta }_S\cos h}{\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S\cos h}$ .(96)

Note that the quantity $F_{S\downarrow }/S$ in Equation (95) is the normalized solar irradiance. Its latitudinal-annual distribution is illustrated in Figure 4.

If we consider the local solar noon at which the solar irradiance reaches its local maximum, we have $h=0$ and, hence,

$\frac{\partial F_{S\downarrow }}{\partial {\delta }_S}=F_{S\downarrow }\frac{\sin \varphi \cos {\delta }_S-\cos \varphi \sin {\delta }_S}{\sin \varphi \sin {\delta }_S+\cos \varphi \cos {\delta }_S}=F_{S\downarrow }\frac{\sin \left(\varphi -{\delta }_S\right)}{\cos \left(\varphi -{\delta }_S\right)}=F_{S\downarrow }\tan \left(\varphi -{\delta }_S\right)$ .(97)

Thus, Equation may be written as

${\zeta }_{F_S}=\frac{{\sigma }_{F_S}}{F_{S\downarrow }}={\left(4\frac{{\sigma }_r^2}{r^2}+\frac{{\sigma }_S^2}{S^2}+{\tan }^2\left(\varphi -{\delta }_S\right){\sigma }_{{\delta }_S}^2\right)}^{1/2}$ ,(98)

where ${\zeta }_{F_S}$ is the relative uncertainty. We choose the following uncertainties: ${\sigma }_r=\pm 10000\text{km}\cong \pm 6.6846\times {10}^{-5}\text{AU}$ , ${\sigma }_S=\pm 1\text{W}\cdot {\text{m}}^{-2}$ and ${\sigma }_{{\delta }_S}=\pm 30"=\pm 0.0083˚$ . The terms $4{\sigma }_r^2/r^2$ and ${\sigma }_S^2/S^2$ are negligible. Shortly after the March equinox the difference $\varphi -{\delta }_S$ amounts to 89˚. Thus, ${\tan }^2\left(\varphi -{\delta }_S\right){\sigma }_{{\delta }_S}^2$ amounts to $6.89\times {10}^{-5}$ . This value is much larger than the values of the other terms. At the June-solstice the difference $\varphi -{\delta }_S$ amounts to 66.56˚. Thus, the value of ${\tan }^2\left(\varphi -{\delta }_S\right){\sigma }_{{\delta }_S}^2\cong 1.12\times {10}^{-7}$ is comparable with the two other terms.

NOTES

¹This is not true, for instance, at the perihelion transit of Mercury. Near perihelion, Mercury’s orbital angular velocity slightly exceeds that of its rotation, and the apparent planetocentric motion of the Sun is retrograde resulting into brief secondary sunrise and sunset at the longitude of 90˚W [31-34].

²President of EIKE is Dr. Holger Thuss, a historian, listed as a policy expert by the Heartland Institute. The Heartland Institute also co-sponsored some of EIKE’s international conferences on climate and energy.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Meech, L.W. (1856) On the Relative Intensity of the Heat and Light of the Sun upon Different Latitudes of the Earth. Smithsonian Contributions to Knowledge, 9, 64 p.
[2]	Wiener, C. (1877) Ueber die Stärke der Bestrahlung der Erde durch die Sonne in den verschiedenen Breiten und Jahreszeiten. Zeitschrift fuer Mathematik und Physik, 22, 341-368.
[3]	Wiener, C. (1879) Ueber die Stärke der Bestrahlung der Erde durch die Sonne in den verschiedenen Breiten und Jahreszeiten. Meteorologische Zeitschrift, XIV, 113-130.
[4]	Milankovitch, M. (1920) Théorie mathématique des phénomènes thermiques produits par la radiation solaire. Gauthier-Villars.
[5]	Milankovitch, M. (1941) Kanon der Erdbestrahlungen und seine Anwendung auf das Eiszeitenproblem. Royal Serbian Academy, Section of Mathematical and Natural Sciences, 33, 633 p.
[6]	List, R.J. (1951) Smithsonian Meteorological Tables. 6th Edition, Robert J. List. Smithsonian Institution.
[7]	Peixoto, J.P. and Oort, A.H. (1992) Physics of Climate. American Institute of Physics.
[8]	Berger, A., Loutre, M. and Tricot, C. (1993) Insolation and Earth’s Orbital Periods. Journal of Geophysical Research: Atmospheres, 98, 10341-10362. https://doi.org/10.1029/93jd00222
[9]	Liou, K.N. (2002) An Introduction to Atmospheric Radiation. 2nd Edition, Academic Press.
[10]	Fu, Q. (2003) Radiation (SOLAR). In: Holton, J.R., Ed., Encyclopedia of Atmospheric Sciences, Elsevier, 1859-1863. https://doi.org/10.1016/b0-12-227090-8/00334-1
[11]	Kramm, G. and Dlugi, R. (2011) Scrutinizing the Atmospheric Greenhouse Effect and Its Climatic Impact. Natural Science, 3, 971-998. https://doi.org/10.4236/ns.2011.312124
[12]	Smulsky, J.J. and Krotov, O.I. (2014) New Computing Algorithm of the Earth’s Insolation. Applied Physics Research, 6, 56-82. https://doi.org/10.5539/apr.v6n4p56
[13]	Fortak, H. (1971) Meteorologie. Deutsche Buch-Gemeinschaft.
[14]	Raschke, E., Vonder Haar, T.H., Bandeen, W.R. and Pasternak, M. (1973) The Annual Radiation Balance of the Earth-Atmosphere System during 1969-70 from Nimbus 3 Measurements. Journal of the Atmospheric Sciences, 30, 341-364. https://doi.org/10.1175/1520-0469(1973)030<0341:tarbot>2.0.co;2
[15]	Berger, A. and Yin, Q. (2012) Astronomical Theory and Orbital Forcing. In: Matthews, J.A., Ed., The SAGE Handbook of Environmental Change: Volume 1, SAGE Publications Ltd, 405-425. https://doi.org/10.4135/9781446253045.n19
[16]	Kopp, G. (2023) Daily Solar Flux as a Function of Latitude and Time. Solar Energy, 249, 250-254. https://doi.org/10.1016/j.solener.2022.11.022
[17]	Park, R.S., Folkner, W.M., Williams, J.G. and Boggs, D.H. (2021) The JPL Planetary and Lunar Ephemerides DE440 and De441. The Astronomical Journal, 161, 105. https://doi.org/10.3847/1538-3881/abd414
[18]	Kaplan, G.H., Bartlett, J., Monet, A., Bangert, J. and Puatua, W. (2009) User’s Guide to NOVAS Version F3.0—Naval Observatory Vector Astrometry Software. United States Naval Observatory.
[19]	Kaplan, G.H., Bartlett, J., Monet, A., Bangert, J. and Puatua, W. (2011) User’s Guide to NOVAS Version F3.1—Naval Observatory Vector Astrometry Software. United States Naval Observatory.
[20]	Laue, E.G. and Drummond, A.J. (1968) Solar Constant: First Direct Measurements. Science, 161, 888-891. https://doi.org/10.1126/science.161.3844.888
[21]	Kopp, G. and Lean, J.L. (2011) A New, Lower Value of Total Solar Irradiance: Evidence and Climate Significance. Geophysical Research Letters, 38. https://doi.org/10.1029/2010gl045777
[22]	Kopp, G., Fehlmann, A., Finsterle, W., Harber, D., Heuerman, K. and Willson, R. (2012) Total Solar Irradiance Data Record Accuracy and Consistency Improvements. Metrologia, 49, S29-S33. https://doi.org/10.1088/0026-1394/49/2/s29
[23]	Dudok de Wit, T., Kopp, G., Fröhlich, C. and Schöll, M. (2017) Methodology to Create a New Total Solar Irradiance Record: Making a Composite Out of Multiple Data Records. Geophysical Research Letters, 44, 1196-1203. https://doi.org/10.1002/2016gl071866
[24]	Wu, C., Krivova, N.A., Solanki, S.K. and Usoskin, I.G. (2018) Solar Total and Spectral Irradiance Reconstruction over the Last 9000 Years. Astronomy & Astrophysics, 620, A120. https://doi.org/10.1051/0004-6361/201832956
[25]	Kondratyev, K.Y. (1969) Radiation in the Atmosphere. Academic Press.
[26]	Emilio, M., Kuhn, J.R., Bush, R.I. and Scholl, I.F. (2012) Measuring the Solar Radius from Space during the 2003 and 2006 Mercury Transits. The Astrophysical Journal, 750, Article 135. https://doi.org/10.1088/0004-637x/750/2/135
[27]	Haltiner, G.J. and Martin, F.L. (1957) Dynamical and Physical Meteorology. McGraw-Hill Book Company.
[28]	Möller, F. (1973) Einführung in die Meteorologie. Bibliographisches Institut.
[29]	Iqbal, M. (1983) An Introduction to Solar Radiation. Academic Press.
[30]	Vardavas, I. and Taylor, F. (2007) Radiation and Climate. Oxford University Press. https://doi.org/10.1093/acprof:oso/9780199227471.001.0001
[31]	Soter, S. and Ulrichs, J. (1967) Rotation and Heating of the Planet Mercury. Nature, 214, 1315-1316. https://doi.org/10.1038/2141315a0
[32]	Morrison, D. (1970) Thermophysics of the Planet Mercury. Space Science Reviews, 11, 271-307. https://doi.org/10.1007/bf00241524
[33]	Vasavada, A. (1999) Near-surface Temperatures on Mercury and the Moon and the Stability of Polar Ice Deposits. Icarus, 141, 179-193. https://doi.org/10.1006/icar.1999.6175
[34]	Bauch, K.E., Hiesinger, H., Greenhagen, B.T. and Helbert, J. (2021) Estimation of Surface Temperatures on Mercury in Preparation of the MERTIS Experiment Onboard Bepicolombo. Icarus, 354, Article ID: 114083. https://doi.org/10.1016/j.icarus.2020.114083
[35]	Smulsky, J.J. (2016) Fundamental Principles and Results of a New Astronomic Theory of Climate Change. Advances in Astrophysics, 1, 1-21. https://doi.org/10.22606/adap.2016.11001
[36]	Yallop, B.D. and Hohenkerk, C.Y. (1992) Astronomical Phenomena. In: Seidelmann, P.K., Ed, Explanatory Supplement to the Astronomical Almanac, University Science Books, 475-503.
[37]	Yallop, B.D., Hohenkerk, C.Y., and Bell, S.A. (2013) Astronomical Phenomena. In: Urban, S.E. and Seidelmann, P.K., Eds., Explanatory Supplement to the Astronomical Almanac, University Science Books, 505-528.
[38]	Berger, A. (1978) Long-Term Variations of Daily Insolation and Quaternary Climatic Changes. Journal of the Atmospheric Sciences, 35, 2362-2367. https://doi.org/10.1175/1520-0469(1978)035<2362:ltvodi>2.0.co;2
[39]	Lindzen, R.S. and Pan, W. (1994) A Note on Orbital Control of Equator-Pole Heat Fluxes. Climate Dynamics, 10, 49-57. https://doi.org/10.1007/bf00210336
[40]	Berger, A. (1988) Milankovitch Theory and Climate. Reviews of Geophysics, 26, 624-657. https://doi.org/10.1029/rg026i004p00624
[41]	Monin, A.S. and Shishkov, Y.A. (2000) Climate as a Problem of Physics. Uspekhi Fizicheskih Nauk, 170, 419-445. https://doi.org/10.3367/ufnr.0170.200004d.0419
[42]	Bretagnon, P. (1974) Termes a longues periodes dans le systeme solaire. Astronomy & Astrophysics, 30, 141-154.
[43]	Lindzen, R.S. (1994) Climate Dynamics and Global Change. Annual Review of Fluid Mechanics, 26, 353-378. https://doi.org/10.1146/annurev.fl.26.010194.002033
[44]	Crowley, T.J. and North, G.R. (1991) Paleoclimatology. Oxford University Press.
[45]	Charlot, P., Jacobs, C.S., Gordon, D., Lambert, S., de Witt, A., Böhm, J., et al. (2020) The Third Realization of the International Celestial Reference Frame by Very Long Baseline Interferometry. Astronomy & Astrophysics, 644, A159. https://doi.org/10.1051/0004-6361/202038368
[46]	Konopliv, A.S., Park, R.S. and Folkner, W.M. (2016) An Improved JPL Mars Gravity Field and Orientation from Mars Orbiter and Lander Tracking Data. Icarus, 274, 253-260. https://doi.org/10.1016/j.icarus.2016.02.052
[47]	Urban, S.E. and Seidelmann, P.K. (2013) Explanatory Supplement to the Astronomical Almanac. University Science Books.
[48]	Seidelmann, P.K. (1992) Explanatory Supplement to the Astronomical Almanac. University Science Books.
[49]	Lambert, J.H. (1779) Pyrometrie, oder vom Maaße des Feuers und der Wärme. Haude und Spener.
[50]	Fritz, S. (1951) Solar Radiant Energy and Its Modification by the Earth and Its Atmosphere. In: Malone, T.F., Ed., Compendium of Meteorology, American Meteorological Society Boston, 13-33.
[51]	Wallace, J.M. and Hobbs, P.V. (2006) Atmospheric Science: An Introductory Survey. Academic Press.
[52]	Peixóto, J.P. and Oort, A.H. (1984) Physics of Climate. Reviews of Modern Physics, 56, 365-429. https://doi.org/10.1103/revmodphys.56.365
[53]	Kramm, G., Dlugi, R. and Mölders, N. (2017) Using Earth’s Moon as a Testbed for Quantifying the Effect of the Terrestrial Atmosphere. Natural Science, 9, 251-288. https://doi.org/10.4236/ns.2017.98026
[54]	Wiin-Nielsen, A. and Chen, T.C. (1993) Fundamentals of Atmospheric Energetics. Oxford University Press.
[55]	Hantel, M. and Haimberger, L. (2016) Grundkurs Klima. Springer-Verlag.
[56]	Spitaler, R. (1886) Die Wärmeverteilung auf der Erdoberfläche, in Denkschriften der mathem.-naturwiss. Klasse der Kaiserl. Akademie der Wissenschaften in Wien, 1-20.
[57]	Kramm, G., Dlugi, R., and Zelger, M. (2009) Comments on the “Proof of the Atmospheric Greenhouse Effect” by Arthur P. Smith. arXiv: 0904.2767v3. http://arxiv.org/abs/0904.2767v3
[58]	Riley, K.F., Hobson, M.P., and Bence, S.J. (1998) Mathematical Methods for Physics and Engineering. Cambridge University Press.
[59]	Spitaler, R. (1940) Die Bestrahlung der Erde durch die Sonne und die Temperaturverhältnisse in der quatären Eiszeit. Abhandlungen der deutschen Gesellschaft der Wissenschaften und Künste in Prag, Mathematisch-Naturwissenschaftliche Abteilung, III, 78.
[60]	Pumfrey, S. (2024) The Bright Side of the Moon: The Historical Journey of Our Tidally Locked Satellite. Journal of Physics: Conference Series, 2877, Article ID: 012034. https://doi.org/10.1088/1742-6596/2877/1/012034
[61]	Gerlich, G. and Tscheuschner, R.D. (2009) Falsification of the Atmospheric CO2 Greenhouse Effects within the Frame of Physics. International Journal of Modern Physics B, 23, 275-364. https://doi.org/10.1142/s021797920904984x
[62]	Locarnini, R.A., Mishonov, A.V., Baranova, O.K., Boyer, T.P., Zweng, M.M., Garcia, H.E., Reagan, J.R., Seidov, D., Weathers, K.W., Paver, C.R. and Smolyar, I.V. (2019) World Ocean Atlas 2018, Volume 1: Temperature. US Department of Commerce, National Oceanic and Atmospheric Administration.
[63]	Kiehl, J.T. and Trenberth, K.E. (1997) Earth’s Annual Global Mean Energy Budget. Bulletin of the American Meteorological Society, 78, 197-208. https://doi.org/10.1175/1520-0477(1997)078<0197:eagmeb>2.0.co;2
[64]	Budó, A. (1990) Theoretische Mechanik. VEB Deutscher Verlag der Wissenschaften.
[65]	Greiner, W. (1977) Theoretische Physik, Band 1, Mechanik I. Verlag Harry Deutsch.
[66]	Seidelmann, P.K. (2013) Introduction to Positional Astronomy. In: Urban, S.E. and Seidelmann, P.K., Eds., Explanatory Supplement to the Astronomical Almanac, University Science Books, 1-44.
[67]	Hilton, J.L., Capitaine, N., Chapront, J., Ferrandiz, J.M., Fienga, A., Fukushima, T., et al. (2006) Report of the International Astronomical Union Division I Working Group on Precession and the Ecliptic. Celestial Mechanics and Dynamical Astronomy, 94, 351-367. https://doi.org/10.1007/s10569-006-0001-2
[68]	Kaplan, G.H. (2005) The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models. United States Naval Observatory.
[69]	Mittelstaedt, P. (1970) Klassische Mechanik. Bibliographisches Institut.
[70]	Hohenkerk, C.Y., Yallop, B.D., Smith, C.A. and Sinclair, A.T. (1992) Celestial Reference Systems. In: Seidelmann, P.K., Ed, Explanatory Supplement to the Astronomical Almanac, University Science Books, 95-198.
[71]	Yoder, C.F. (2013) Astrometric and Geodetic Properties of Earth and the Solar System. In: Ahrens, T.J., Ed., Global Earth Physics: A Handbook of Physical Constants, American Geophysical Union, 1-31. https://doi.org/10.1029/rf001p0001
[72]	Hohenkerk, C.Y. (2013) Positions. In: Urban, S.E. and Seidelmann, P.K., Eds., Explanatory Supplement to the Astronomical Almanac, University Science Books.
[73]	Seidelmann, P.K. and Wilkins, G.A. (1992) Introduction to Positional Astronomy In: Seidelmann, P.K., Ed., Explanatory Supplement to the Astronomical Almanac, University Science Books, 1-38.
[74]	Folkner, W.M., Williams, J.G., Boggs, D.H., Park, R.S. and Kuchynka, P. (2014) The Planetary and Lunar Ephemerides DE430 and DE431, in IPN Progress Report. Jet Propulsion Laboratory, California Institute of Technology, 81.
[75]	Hohenkerk, C.Y. (2006) Implementation of the New Nomenclature in the Astronomical Almanac. Proceedings of the International Astronomical Union, 2, 476-476. https://doi.org/10.1017/s1743921307011507
[76]	McCarthy, D. and Guinot, B. (2013) Time. In: Urban, S.E. and Seidelmann, P.K., Eds., Explanatory Supplement to the Astronomical Almanac, University Science Books.
[77]	Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.T. (2003) Expressions for the Celestial Intermediate Pole and Celestial Ephemeris Origin Consistent with the IAU 2000A Precession-Nutation Model. Astronomy & Astrophysics, 400, 1145-1154. https://doi.org/10.1051/0004-6361:20030077
[78]	Mölders, N. and Kramm, G. (2014) Lectures in Meteorology. Springer International Publishing.
[79]	Standish, E.M. and Williams, E.M. (2013) Orbital Ephemerides of the Sun, Moon, and planets. In: Urban, S.E. and Seidelmann, P.K., Eds., Explanatory Supplement to the Astronomical Almanac, University Books, 305-345.
[80]	Unsöld, A. (1967) Der neue Kosmos. Springer-Verlag.
[81]	von Bezold, W. (1906) Über Strahlungsnormalen und Mittellinien der Temperatur. In: Pernter, J.M., Von Hann, J., and Hellmann, G., Eds., Meteorologische Zeitschrift: Hann-Band zum vierzigjährigen Redaktionsjubiläum J. Hann’s von Freunden und Kollegen Gewidmet, F. Vieweg, 279-287.
[82]	Defant, A. and Obst, E. (1923) Lufthülle und Klima. F. Deuticke, Leipzig.
[83]	Mölders, N., Jankov, M. and Kramm, G. (2005) Application of Gaussian Error Propagation Principles for Theoretical Assessment of Model Uncertainty in Simulated Soil Processes Caused by Thermal and Hydraulic Parameters. Journal of Hydrometeorology, 6, 1045-1062. https://doi.org/10.1175/jhm455.1
[84]	Iorio, L. (2019) Calculation of the Uncertainties in the Planetary Precessions with the Recent EPM2017 Ephemerides and Their Use in Fundamental Physics and Beyond. The Astronomical Journal, 157, Article 220. https://doi.org/10.3847/1538-3881/ab19bf