Earth Orbit v 2 . 1 : a 3-D visualization and analysis model of Earth ’ s orbit , Milankovitch cycles and insolation

Introduction Conclusions References Tables Figures


Introduction
The astrophysical characteristics of our star, the Sun, determine to first order the continuously habitable zone around it (Kasting et al., 1993;Kasting, 2010), in which rocky planets are able to maintain liquid water on their surface and sustain life.The surface temperature of a planet depends to first order upon the incoming flux of solar radiation (insolation) to its surface.Additionally, energy for our metabolism (and most modern economic activities) is obtained exclusively from the Sun via the process of oxygenic photosynthesis performed by green terrestrial plants and marine phytoplankton.The high oxygen content of Earth's atmosphere, necessary for the evolution of placental mammals (Falkowski et al., 2005), is due to billions of years of photosynthesis and the geological burial of reduced carbon equivalents (Falkowski and Godfrey, 2008;Falkowski and Isozaki, 2008;Kump et al., 2010).Thus, the Sun is central to climate formation and stability and to our evolution and continued existence as a species.
T. S. Kostadinov and R. Gilb: Earth Orbit v2.1 The temporal and spatial patterns of insolation and their variability on various scales determine climatic stability over geologic time, as well as climate characteristics such as diurnal, seasonal and pole to equator temperature contrasts, all of which influence planetary habitability.Insolation can change due to changes in the luminosity of the Sun itself.This can happen due to the slow increase of solar luminosity that gives rise to the faint young Sun paradox (Kasting, 2010;Kump et al., 2010), or it can happen on much shorter timescales such as the 11-year sunspot cycle (Fröhlich, 2013;Hansen et al., 2013).
Importantly, insolation is also affected by the orbital elements of the planet.According to the astronomical theory of climate, quasi-periodic variations in Earth's orbital elements cause multi-millennial variability in the spatio-temporal distributions of insolation, and thus provide an external forcing and pacing to Earth's climate (Milankovitch, 1941;Berger 1988;Berger and Loutre, 1994;Berger et al., 2005).These periodic orbital fluctuations are called Milankovitch cycles, after the Serbian mathematician Milutin Milanković who was instrumental in developing the theory (Milankovitch, 1941).Laskar et al. (2004) provide a brief historical overview of the main contributions leading to the pioneering work of Milanković.There are three Milankovitch orbital parameters: orbital eccentricity (main periodicities of ∼ 100 and 400 kyr (1 kyr = one thousand years)), precession (quantified as the longitude of perihelion relative to the moving vernal equinox, main periodicities ∼ 19 and 23 kyr) and obliquity of the ecliptic (main periodicity 41 kyr) (Berger, 1978a).Obliquity is strictly speaking a rotational, rather than an orbital parameter; however, we refer to it here either as an orbital or Milankovitch parameter, for brevity.
The pioneering work by Hays et al. (1976) demonstrated a strong correlation between these cycles and paleoclimatological records.Since then, multiple analyses of paleoclimate records have been found to be consistent with Milankovitch forcing (e.g., Imbrie et al., 1992;Rial, 1999;Lisiecki and Raymo, 2005).Notably, the glacial-interglacial cycles of the Quaternary have been strongly linked to orbital forcing, particularly summertime insolation at high northern latitudes (Milankovitch, 1941;Berger 1988;Berger and Loutre, 1994;Bradley, 2014, and references therein).Predicting the Earth system response to orbital forcing (including glacier growth and melting) is not trivial, and there are challenges in determining which insolation quantity (i.e., integrated over what time and space scales) is responsible for paleoclimate change, e.g., peak summer insolation intensity, or overall summertime-integrated insolation at northern latitudes (Imbrie et al., 1993;Lisiecki et al., 2008;Huybers, 2006;Huybers and Denton, 2008;Bradley, 2014).Moreover, some controversies related to the astronomical theory remain, notably the 100 kyr problem, or the so-called mid-Pleistocene transition.This refers the fact that the geological record indicates that the last ca.one million years have been dominated by 100 kyr glacial-interglacial cycles, a gradual switch from the previously dominant 41 kyr periodicity.This transition cannot be explained by orbital forcing alone, as there was actually a decrease in the 100 kyr eccentricity band variance during this period (e.g., Imbrie et al., 1993;Loutre et al., 2004;Berger et al., 2005;Bradley, 2014, and references therein).Current consensus focuses on the explanation that the mid-Pleistocene transition is due to factors within the Earth system itself, rather than astronomical factors -e.g., internal climate system oscillations, nonlinear responses due to the continental ice sheet size, or CO 2 degassing from the Southern Ocean (Bradley, 2014, Sects. 6.3.3 and 6.3.4, and references therein).Finally, alternative astronomical influences on climate have also been proposed, such as the influence of the orbital inclination cycle (Muller and McDonald, 1997).
The Milankovitch cycles are due to complex gravitational interactions between the bodies of the solar system.Astronomical solutions for the values of the Milankovitch orbital parameters have been derived by Berger (1978a) and Berger (1978b), referred to henceforth as Be78 (valid for 1000 kyr before and after present), and Laskar et al. (2004), referred to henceforth as La2004 (valid for 101 000 kyr before present to 21 000 kyr after present).Here the present is defined as the start of Julian epoch 2000 (J2000), i.e., the Gregorian calendar date of 1 January, 2000 at 12:00 UT (Universal Time, that is, mean solar time at the Prime Meridian + 12 h) (Meeus, 1998).There are several other solutions as well, for example, Berger and Loutre (1992) and Laskar et al. (2011).These astronomical solutions are crucial for paleoclimate and climate science, as they enable the computation of insolation at any latitude and time period in the past or future within the years spanned by the solutions (Berger and Loutre, 1994;Berger et al., 2010;Laskar et al., 2004), and subsequently the use of this insolation in climate models as forcing (e.g., Berger et al., 1998).Climate models are an important method for testing the response of the Earth system to Milankovitch forcing.
While most Earth science students and professionals are well aware of Earth's orbital configuration and the basics of the Milankovitch cycles, the details of both and the way the Milankovitch orbital elements influence spatio-temporal patterns of insolation on various time and space scales remain elusive.It is difficult to appreciate the pivotal importance of Kepler's laws of planetary motion in controlling the effects of Milankovitch cycles on insolation patterns.The three-dimensional nature of Earth's orbit, the vast range of space and timescales involved, and the geometric details are complex, and yet those same factors present themselves to computer modeling and 3-D visualization.Here, we present "Earth Orbit v2.1": an astronomically precise and accurate 3-D visualization and analysis model of Earth's orbit, Milankovitch cycles, and insolation.The model is envisioned for both research and pedagogical applications and offers 3-D visualizations of Earth's orbital geometry, Milankovitch parameters and the ensuing insolation forcing.It is developed in MATLAB ® and has an intuitive, user-friendly graphical user interface (GUI) (Fig. 1).Users are presented with a choice between the Be78 and La2004 astronomical solutions for eccentricity, obliquity and precession.A "demo" mode is also available, which allows the three Milankovitch parameters to be varied independently of each other (and exaggerated over much larger ranges than the naturally occurring ones), so users can isolate the effects of each parameter on orbital geometry, the seasons, and insolation.Users select a calendar date and the Earth is placed in its orbit using Kepler's laws; the calendar can be started on either vernal equinox (20 March) or perihelion (3 January).A 3-D orbital configuration visualization, as well as spatio-temporal surface and line plots of insolation and insolation anomalies (with respect to J2000) on various scales are then produced.Below, we first describe the model parameters and implementation.We then detail the model user interface, provide instructions on its capabilities and use, and describe the output.We then present successful model validation results, which are comparisons to existing independently derived insolation, solar declination and season duration values.Finally, we conclude with brief analysis of sources of uncertainty.Throughout, we provide examples of the pedagogical value of the model.
Various insolation solutions and visualizations exist (Berger, 1978a;Rubincam, 1994 (however, see response of Berger, 1996); Laskar et al., 2004;Archer, 2013;Huybers, 2006).Notably, the AnalySeries software (Paillard et al., 1996;Paillard, 2014) shares many of the functionalities presented here and offers many additional ones, such as paleoclimatic time-series analysis and many more choices for insolation computation.Importantly, the model presented here was developed independently from AnalySeries (or other similar efforts) and computes insolation from first principles of orbital mechanics (Kepler's laws) and irradiance propagation, using exclusively its own internal geometry.The only model inputs are the three Milankovitch orbital parameters, either real astronomical solutions (Be78 or La2004) or userentered demo values.No insolation computation code from the above-cited existing solutions has been used, so comparison with these solutions constitutes independent model verification, referred to here as validation, because we consider the La2004 and Meeus (1998) solutions the geophysical truth (Sect.5).
The unique contribution of our model consists of the combination of the following features: a. central to the whole model is a user-controllable, 3-D pan-tilt-zoom plot of the actual Earth orbit; b. an interactive user-friendly GUI that serves as a singleentry control panel for the entire model and makes it suitable for use by non-programmers and friendly to didactic applications; c. the Milankovitch cycles are incorporated explicitly and insolation is output according to real or user-selected demo orbital elements, which d. allows users to enter exaggerated orbital parameters independently of each other and isolate their effects on insolation, as well as view the orbit with exaggerated eccentricity; e. the source code is published and advanced users can check its logic, as well as modify it and adapt it; and f. the software is platform-independent.
The issue of climate change has come to the forefront of Earth science and policy and it is arguably the most important global issue of immediate and long-term consequences (e.g., IPCC, 2013).Earth's climate varies naturally over multiple timescales, from decadal to hundreds of millions of years (e.g., Kump et al., 2010).It is thus crucial to understand natural climate forcings, their timescales, and the ensuing response of the Earth system.In addition, detailed understanding of the Sun's daily path in the sky and the patterns of insolation have become important to increasing numbers of students and professionals because of the rise in usage of solar power (thermal and photovoltaic).We submit that the model presented here can enhance understanding of all of these important subject areas.

Key definitions, model parameters and implementation
The model input parameters, and their values and units, are summarized in Table 1.The following definitions, discussion and symbols are consistent with those of Berger et al. (2010).The reader is referred to their Fig. 1.According to Kepler's first law of planetary motion, Earth's orbit is an ellipse, and the Sun is in one of its foci (e.g., Meeus, 1998).Orbital eccentricity, e (Table 1), is a measure of the deviation of Earth's orbital ellipse from a circle and is defined as e = 1 − b 2 /a 2 , where a is the semi-major axis (Table 1) and b is the semi-minor axis of the orbital ellipse (e.g., Berger and Loutre, 1994).The semi-major axis is equal to about 1 AU (Meeus, 1998;Standish et al., 1992) and determines the size of the orbital ellipse and thus the orbital period of Earth; it is considered a fixed constant in the model, as its variations are extremely small (Berger et al., 2010;Laskar et al., 2004, their Fig. 11).Various orbital period definitions are possible; here, the sidereal period is used as a model constant (Meeus, 1998).Thus, Kepler's third law of planetary motion is implicit in these two constant definitions and is not included explicitly elsewhere in model logic.The obliquity of the ecliptic, ε, is the angle between the direction of Earth's axis of rotation and the normal to the orbital plane, or the ecliptic (Table 1).Eccentricity and obliquity are two of the three Milankovitch orbital parameters.The third Milankovitch orbital parameter, precession, is the most challenging for instruction and visualization.There are two separate kinds of precession that combine to create a climatic effect -precession of the equinoxes (also termed axial precession), and apsidal precession, that is, precession of the perihelion in the case of Earth's orbit.Axial precession refers to the wobbling of Earth's axis of rotation that slowly changes its absolute orientation in space with respect to the distant stars.The axis or rotation describes a cone (one in each hemisphere) in space with a periodicity of about 26 000 years (Berger and Loutre, 1994).This is the reason why the star α UMi (present-day Polaris, or the North Star), has not and will not always be aligned with the direction of the North Pole.Also, due to axial precession, the point of vernal equinox in the sky moves with respect to the distant stars and occurs in successively earlier zodiacal constellations.Axial precession is clockwise as viewed from above the North Pole, hence the north celestial pole describes a counterclockwise motion as viewed by an observer looking in the direction of the north ecliptic pole.Precession of the perihelion refers to the gradual rotation of the line joining aphelion and perihelion, with respect to the distant stars (or the reference equinox of a given epoch) (Berger, 1978a;Berger and Loutre, 1994).
Axial precession and precession of the perihelion combine to modulate the relative position of the equinoxes and solstices (i.e., the seasons) with respect to perihelion, which is what is relevant for insolation and climate.This climatically relevant precession is implemented in the model and is quantified via the longitude of perihelion, ω, which is the angle between the directions of the moving fall equinox and perihelion at a given time, measured counterclockwise  Laskar et al. (2004) in the plane of the ecliptic (Berger et al., 2010).Because both perihelion and equinox move, the longitude of perihelion will have a different (shorter) periodicity than one full cycle of axial wobbling alone (Berger and Loutre, 1994).
The direction of Earth's radius vector when Earth is at fall equinox (∼ 22 September) is referred to as the direction of fall equinox above.This is the direction with respect to the distant stars where the Sun would be found on its annual motion on the ecliptic on 20 March -that is, at vernal equinox.
In other words, that is the direction of the vernal point in the sky (Berger et al., 2010; their Fig. 1 and Appendix B), the origin of the right ascension coordinate.This distinction between vernal equinox and the direction of the vernal point can cause confusion, especially since the exact definition of longitude of perihelion can vary (e.g., c.f. Berger, 1978a;Berger et al., 1993Berger et al., , 2010;;Berger and Loutre, 1994;Joussaume and Braconnot, 1997) and the longitude of perihelion can also be confused with the longitude of perigee, ω = ω + 180 • , which is the angle between the directions of vernal equinox and perihelion, measured counterclockwise as viewed from the direction of the North Pole, in the plane of the orbit (Berger et al., 2010).Here, we use the terminology and definitions of Berger et al. (2010).
The magnitude of the climatic effect of precession is modulated by eccentricity.In the extreme example, if eccentricity were exactly zero, the effects of precession would be null.Climatic precession, e sin ω, is the parameter that quantifies precession and determines season lengths, the Earth-Sun distance at summer solstice (Berger and Loutre, 1994) and various important insolation quantities (Berger et al., 1993, their Table 1).This interplay between eccentricity and precession presents an important way to introduce both concepts pedagogically and to test student comprehension.
The solar "constant", S o , is defined here as the total solar irradiance (TSI) on a flat surface perpendicular to the solar rays at a reference distance of exactly 1 AU (Table 1).As Berger et al. (2010) note, due to eccentricity changes, the mean distance from the Earth to the Sun over a year is not constant on geologic timescales.It also matters how this mean distance is defined -for example, over time (mean anomaly) vs. over angle (true anomaly).True and mean anomaly are defined below in Sects.2.1 and 2.2, respectively.If S o is defined to be the irradiance from the Sun at the mean Earth-Sun distance, then it is indeed not a true constant.As used here, S o is a true model constant as long as the luminosity of the Sun itself is assumed constant.The default value is chosen to be 1366 W m −2 (Fröhlich, 2013).Recent evidence suggests that the appropriate value may actually be about 1361 W m −2 (Kopp and Lean, 2011).Users can change the value of S o independently of other model inputs in order to study the effects of changes in absolute solar luminosity -for example, in order to simulate the faint young Sun (e.g., Kasting, 2010) or the sunspot cycle (e.g., Hansen et al., 2013).

Model coordinate system; Sun-Earth geometry parameterization; solar declination
According to Kepler's first law of planetary motion, Earth orbits the Sun in an ellipse, and the Sun is in one of the ellipse's foci.The heliocentric equation of the orbital ellipse in polar form is given by (Meeus, 1998;his Eq. 30.3) In the above, the Sun is at the origin of the coordinate system; a is the semi-major axis of the orbital ellipse; e is eccentricity; ν is true anomaly; and r is Earth's instantaneous radius vector, that is, the vector originating at the Sun and ending at the instantaneous planetary position.Let r designate the length of Earth's radius vector henceforth.True anomaly, ν, is the angle between the directions of perihelion and the radius vector, subtended at the Sun and measured counterclockwise in the plane of the orbit (e.g., Meeus, 1998, his Ch. 30;Berger et al., 2010, their Fig. 1).The true longitude of the Sun (or simply true longitude) is equal to Earth's true anomaly plus the longitude of perigee (Berger et al., 2010, their Eq. 6).True longitude is the angle Earth has swept from its orbit, subtended at the Sun, since it was last at vernal equinox, and it is equivalent to the angle the Sun has travelled along the ecliptic in the same time.Mean longitude is the longitude of the mean Sun, in an imaginary perfectly circular orbit of the same period, that is, mean longitude is proportional to the passage of time, much like mean anomaly (see Sect. 2.2 below).
In the Earth Orbit v2.1 model, given a user-selected calendar date, true anomaly, ν, is determined by solving the inverse Kepler equation (see Sect. 2.2 below).The Earth's radius vector is then solved for using Eq. ( 1) above.Because the main model coordinate system is heliocentric Cartesian, the (r, ν) pair of polar coordinates is then transformed to Cartesian (x, y) for plotting.The model's main coordinate system has its origin at the center of the Sun, its positive x axis pointing in the direction of perihelion, its positive y axis pointing 90 degrees counterclockwise in the plane of the ecliptic, and its positive z axis perpendicular to it in the direction of the north ecliptic pole.The Earth is initially parameterized as a sphere in its own geocentric Cartesian coordinate system in terms of its radius and geographic latitude and longitude (corresponding to the two angles of a spherical coordinate system).The Earth's coordinate system's x and y axes are in the plane of the Equator (shown as a black dotted line, Fig. 2), and its z axis is pointing towards the true North Pole and is coinciding with Earth's axis of rotation; these axes are also plotted in black dotted lines; the z axis is lengthened so that it pierces Earth's surface at the Poles, and the North Pole is labeled.Earth is plotted as a transparent mesh so that important orbital elements can be seen through it at various zoom levels (Fig. 2).The color scale of Earth's mesh is just a function of latitude and no day and night sides are explicitly shown.Earth's radius is not to scale with the orbit itself or with the Sun's radius.Thus, the center of Earth has its true geometric orbital position (and is the tip of its instantaneous radius vector); however, the surface of the sphere in the model is arbitrary and must not be interpreted as the true surface onto which insolation is computed, for example.The insolation computations (Sect.2.3) are geocentric.The Sun is also plotted (not to scale) as a sphere centered at the origin of the main model coordinate system.
The Earth is oriented properly in 3-D with respect to the orbital ellipse by using a rotation matrix to rotate its coordinate system.The 3-D rotation matrix is computed using Rodrigues' formula (Belongie, 2013) for 3-D rotation about a given direction by a given angle.The direction about which Earth is rotated is determined by a vector which is always in the orbital plane (k component is zero), and the i and j components are determined by the longitude of perihelion.The angle by which Earth is rotated is determined by obliquity.Thus, the rotation matrix is a function of two of the three Milankovitch parameters and is a valuable and useful instructional tool/concept for lessons in geometry, mathematics, astronomy, physical geography, and climatology.At this point the Earth is correctly oriented in 3-D space with respect to its orbit and the distant stars.Earth is then translated to its proper instantaneous position on its orbit by addition of its radius vector to all relevant Earth-bound model elements (which are then plotted in the main heliocentric coordinate system).
Declination is one of the two spherical coordinates of the equatorial astronomical coordinate system.It is measured along a celestial meridian (hour circle) and is defined as the angle between the celestial Equator and the direction toward the celestial object (Meeus, 1998).Solar declination varies with the seasons, due to obliquity.It is zero at the equinoxes, reaches a maximum of +ε at summer solstice and a minimum of −ε at winter solstice.Solar declination determines the length of day and the daily path of the Sun in the sky at a given latitude (i.e., its altitude and azimuth above the horizon as a function of time).Thus, solar declination determines instantaneous and time-integrated insolation.In turn, solar declination and its evolution over the course of a year are a function of the orbital elements; thus it provides the mathematical and conceptual link between the Milankovitch orbital elements and insolation and climate.Here, we compute instantaneous solar declination using the angle between the direction of the North Pole and Earth's radius vector, calculated using their dot product.Thus, we explicitly compute solar declination from the geometry of the model and it is a model emergent property rather than prescribed a priori; therefore, this also applies to insolation computations (Sects.2.3 and 5).

Implementation of Kepler's second law of planetary motion
The heliocentric position of a planet in an elliptical orbit at a given instant of time is given in terms of its true anomaly, ν -see Eq. (1) and Sect.2.1 above.True anomaly can also be thought of as the angle (subtended at the Sun) which the planet has "swept" from its orbit since last perihelion passage.Kepler's second law of planetary motion states that the planet will "sweep" equal areas of its orbit in equal intervals of time and governs the value of true anomaly as a function of time (e.g., Meeus, 1998;Joussaume and Braconnot, 1997).
At non-zero eccentricity, ν is not simply proportional to time since last perihelion passage (time of flight) expressed as a fraction of the orbital period in angular units.The latter quantity is called mean anomaly, M. Kepler's second law is used to relate M and ν, using an auxiliary quantity called eccentric anomaly, E. E and M are related by Kepler's equation (Meeus, 1998;Chapter 30): where e is orbital eccentricity.When E is known, ν can be solved for using (Meeus, 1998; Chapter 30): The forward Kepler problem consists of solving for time of flight, M, given the planetary position, ν.This is straightforward by first solving for E in Eq. ( 3) and using it to solve for M in Eq. ( 2).However, in the most intuitive case, which is implemented here, the user enters a desired date, and the position of the planet has to be determined from the date (i.e., time of flight/mean anomaly M is given, and true anomaly has to be determined).This is referred to as the inverse Kepler problem and amounts to solving for E in Eq. ( 2) and then for ν in Eq. (3).Solving for E is not straightforward, as no analytical solution exists.Numerous numerical methods exist for the solution of the inverse Kepler problem.Here, the binary search algorithm of Sinnott (1985) is used, as given in Meeus (1998).It has the advantage of being computationally efficient, which becomes important when time series of insolation is the desired model output.It also has the distinct advantages of being valid for any value of eccentricity and converging to the exact solution to within the user machine's precision.

Implementation of Insolation Computation
Instantaneous insolation at the top of the atmosphere (TOA) can be computed as where r is the length of the radius vector of Earth expressed in AU, and h is the altitude of the Sun above the horizon (e.g., Berger et al., 2010).Equation ( 4) is an expression of the inverse square law and the Lambert cosine law of irradiance.
The radius vector length is computed in the model for the chosen date (and not for every instant) using Eq. ( 1). S o is the TSI at r o = 1 AU by definition (Sect.2).In this equation insolation, S, is defined as the total (spectrally integrated) solar radiant energy impinging at the TOA on a unit surface area parallel to the mathematical horizon at a given latitude at a given instant.S carries the units of S o , here W m −2 .S needs to be integrated over time and/or space in order to compute insolation quantities of interest.Here, the main discrete time step over which S is computed and output is one 24 h period (i.e., daily insolation).Daily insolation is a function of latitude, date, and S o .The date is associated with a given true anomaly for a given calendar start date and orbital configuration (Joussaume and Braconnot, 1997;Sect. 2.3.1).This determines the current solar declination and the length of the radius vector of Earth (i.e., the Sun-Earth distance).The user inputs the desired latitude, date and TSI, and the rest of the quantities are computed from the model geometry.Solar declination and the latitude determine the daily evolution of solar altitude, h, as a function of time, as follows (e.g., Meeus, 1998): sin h = sin δ sin ϕ + cos δ cos ϕ cos t . (5) In the above equation δ is solar declination, ϕ is geographic latitude on Earth, and t is the hour angle of the Sun.δ is assumed constant for the day of interest, and t is a measure of the progress of time (e.g.Berger et al., 2010).Note that this assumes the time derivative of the solar hour angle is equal to one, i.e. it ignores the time derivative of the Equation of Time (or equivalently, the annual variability in the time derivative of the right ascension of the Sun is ignored).Half the day length, t s , (i.e., the time between local solar noon and sunset), is determined by setting h = 0 • in Eq. (5): In Eq. ( 6) t s is expressed in terms of hour angle of the Sun in angular units.Equation ( 5) is integrated over time (under the assumptions here, equivalently, over hour angle) from solar noon to sunset in order to compute the time-average of the sine of the solar altitude for the given date and latitude: Equation ( 7) is integrated numerically with a very small time step of about 10 s.Because the altitude of the Sun is symmetric about solar noon, it is sufficient to integrate only from solar noon to sunset time.Daily insolation is then computed by using the time-averaged sin h quantity in Eq. ( 4).The results are scaled by multiplying by the actual day length and dividing by 24 h.The resulting quantity represents the mean daily insolation over a full day, which is the standard value used in climate and paleoclimate science (e.g., Laskar, 2014).If this daily insolation is multiplied by 24 h (in seconds), total energy receipt for that day (in J m −2 ) can be calculated.At high latitudes, there are periods of the year with no sunset or no sunrise.These cases depend on the relationship of latitude and solar declination (e.g., Berger et al., 2010).They are handled separately by either integrating Eq. ( 7) over 24 h, or, in the case of no sunrise (polar night), assigning a value of exactly 0 W m −2 to daily insolation.

Integrating insolation over longer time periodscaveats
Because of the varying eccentricity and longitude of perihelion, there is no fixed correspondence between true anomaly (or true longitude) and any one single calendar date, even if one were to define a fixed calendar start date.et al., 2010).If one wishes to make insolation comparisons between different orbital configurations, one must strictly define a calendar start date, and even then insolation will be in phase for different geological periods only for that date (Joussaume and Braconnot, 1997).Thus, the question of "what is insolation on June 20" is ill posed, unless one defines strictly what is meant by the date of 20 June.The problem persists if one wishes to compare insolation integrated over periods of time longer than a day, because over geologic timescales, absolute values and the interval of true longitudes "swept" between two classical calendar dates are not constant.Thus, there are two ways to define a calendar -the classical or fixed-day calendar, in which month lengths follow the present-day configuration and the date of vernal equinox is fixed, or a fixed-angular calendar, which defines months beginning at certain true longitudes (function of true anomaly and the precession phase, see also Sect.2.1) and they can therefore have a different number of days depending on the orbital configuration (Joussaume and Braconnot, 1997;Chen et al., 2011).The time intervals between solstices and equinoxes also varies, because of varying eccentricity and because these intervals happen in different places in the orbit with respect to perihelion.Thus season lengths vary over geologic time.Earth Orbit v2.1 outputs season length in the main GUI to emphasize this important fact (Fig. 1).Earth Orbit v2.1 uses the classical calendar dates (24 h periods) as the user time input, rather than true anomaly or true longitude.This choice is much more intuitive to non-experts, and best serves the educational purposes of the model.The user has as a choice of calendar start date (Sect.3) and true solar longitude is output (Sect.4.2; Fig. 1) to remind users of the above considerations.The effect of calendar choice on insolation phases and comparisons and on climate models is discussed at length by Joussaume andBraconnot (1997), Timm et al. (2008) and Chen et al. (2011).The time step of integration can also influence the results of insolation computations, for example, if annual insolation is averaged with a 5-day step, results are substantially different from the case when a 1-day step is used (not shown).For this reason, the model computes annually averaged insolation at a given latitude by using 1-day steps of integration.Finally, we note that the daily insolation computations of the model are robust and validated for real values of the orbital parameters (Sect.5.1, but see also Sect.6); however, the model currently has limited functionality for making comparisons of insolation integrated over longer time periods over different geologic scales.In order to make such comparisons, the use of the elliptical integrals method of Berger et al. (2010) is recommended, as well as the Laskar et al. (2004) methods, both of which come with accompanying software (Berger, 2014 andLaskar, 2014).In addition, users are referred to the latest version of the AnalySeries software package (Paillard et al., 1996;Paillard, 2014) for additional insolation and time-series options.All of the above can also be used for verification of the output of the model presented here.

Model user interface
The Earth orbit model is provided as Supplement (see section Code availability & license).The model is developed and runs in MATLAB ® .All model control is realized via a single, user-friendly GUI panel (Fig. 1).Users are presented with a choice between the Be78 and the La2004 astronomical solutions for eccentricity, obliquity and precession.A "demo" mode is also available.If a real astronomical solution is chosen, users are asked to input a year before or after present (defined as J2000, i.e., 1 January, 2000 at 12:00 noon UT, see Introduction) for which they wish to run the model.The GUI only allows users to choose years within the respective solution's validity: the Be78 solution is available for 1000 kyr before and after present (J2000), whereas the La2004 solution is available for 101 000 kyr in the past and 21 000 kyr in the future.The La2004 solutions are provided by Laskar (2014) (specifically at http://www.imcce.fr/Equipes/ASD/insola/earth/La2004/index.html) in tabulated form in 1 kyr intervals.The Be78 solutions are obtained by transcribing code from NASA GISS (see Acknowledgements).The model looks up the values of eccentricity, obliquity and precession for the chosen year and solution (using linear interpolation between tabulated years if necessary), and these values are used in subsequent visualizations and analyses.If the user chooses the "demo" mode, they select, independently of each other, the values of the Milankovitch parameters, which can be greatly exaggerated.In this way users can isolate the effects of each parameter on orbital geometry, the seasons, and insolation.The "demo" mode is central to the pedagogical value and applications of the model because it allows users to build and visualize an imaginary orbit of, for example, very high eccentricity while keeping obliquity fixed.Moreover, the model will output all subsequent parameters, such as solar declination, day length and radius vector length, based on this exaggerated imaginary orbit.
Users input the desired calendar date, geographic latitude on Earth (positive degrees in the Northern Hemisphere and negative degrees in the Southern Hemisphere), and desired value of TSI.The calendar date defaults to the current date, latitude defaults to 43 • N, and TSI defaults to 1366 W m −2 (Sect.2).Two choices of calendar start date are available: either fix vernal equinox to be at the beginning of 20 March (default), or fix perihelion to be at the beginning of 3 January.The availability of this choice complicates interpretation of model output; however it has high instructional value.It illustrates that the choice of calendar start date and a calendar system is a human construct, accepted by convention; it is based on the actual year and day length but is relative.This can also help test knowledge of the concepts explained in

Graphical output
The main output of the model is a 3-D plot of Earth's orbital configuration.Figure 2a illustrates the orbital configuration using the contemporary values of the Milankovitch parameters (the La2004 solution for J2000 is shown), for 16 September.The current phase of the precession cycle is such that Northern Hemisphere (NH) winter solstice occurs shortly before perihelion (longitude of perihelion is ∼ 102.9 • ).This results in Northern Hemisphere spring and summer being longer than the respective fall and winter (as shown in the GUI, Fig. 1). Figure 2b illustrates the orbital configuration also on 16 September, but for 10 kyr in the future (also using the La2004 solution).Since this represents about a half of a precession cycle, the timing of the seasons is approximately 180 • out of phase with respect to the contemporary configuration (the longitude of perihelion is ∼ 279.2 • , and Northern Hemisphere summer occurs near perihelion and is the shortest season).Because we chose to fix the calendar start date such that vernal equinox is always on 20 March, and the eccentricity is fairly low, the date 16 September still occurs near the fall equinox, like in the contemporary example.However, because the length of time passing between vernal equinox and fall equinox is now shorter, 16 September almost coincides with fall equinox, unlike the contemporary case.Of course obliquity and eccentricity have also changed 10 kyr in the future, but unlike the longitude of perihelion, their changes are small in absolute terms, and thus this cannot be readily visualized by comparing Fig. 2a and b.This is one reason why it is very useful to have the ability to choose arbitrary independent values of the Milankovitch parameters in the demo mode, constructing an imaginary orbit.Figure 2c illustrates one example of such an imaginary orbit with greatly exaggerated eccentricity (0.6) and obliquity (45 • ) and longitude of perihelion of 225 • , that is, very different from the J2000 values.This imaginary orbit illustrates that the date 1 July can occur in the fall, due to the large eccentricity and the specific phase of precession chosen.Spring lasts only ∼ 20 days in this configuration because it occurs during perihelion passage, where the planet is much faster according to Kepler's second law, as compared to aphelion passage (fall lasts ∼ 229 days in this configuration).Summer lasts about 58 days.Thus, July 1 occurs during the fall season, counterintuitively.Importantly, such an exaggerated eccentricity means that the planet is very close to the Sun during perihelion, and some really high insolation values can occur even at modest solar declinations (e.g., for 29 March, at 43 • N, solar declination is ∼ 27 • , day length is ∼ 16 h, and daily insolation is 3307 W m −2 , far exceeding any contemporary value anywhere on Earth).The reason is that the Sun-Earth distance then is only 0.4 AU, and the distance factor becomes a first-order effect on insolation, whereas it is a second-order factor in the real Earth orbit configuration (angle being the first-order factor, see Eq. 4).
The plots of Fig. 2 have pan-tilt-zoom capability, so users can view the orbital configuration from many perspectives; this is at the core of the pedagogical value of the model.The plot is updated with the current parameter selections by pressing the "Plot/Update Orbit" button.Finally, note that the apparent eccentricity of the orbits also changes with the view angle and the projection onto a 2-D screen.This should not be confused with the intrinsic orbital eccentricity, which can be also judged by the relative distance of the orbital foci (marked with an "x") from the ellipse's center (the intersection of the semi-major and semi-minor axes, red lines in Fig. 2) Users are presented with several options of plotting insolation as a function of time and latitude.First, insolation can be plotted for a single year (using the currently selected Milankovitch parameters) as a function of day of year and latitude (Fig. 3a, upper panel).Insolation anomalies with respect to the J2000 La2004 orbital configuration are also plotted, using S o = 1366 W m −2 (Fig. 3a, lower panel).Anomalies are especially useful when analyzing the effect of changes in insolation on the glacial-interglacial cycles.For example, the anomalies at 65 • N during summer months 115 kyr before present (Fig. 3a, lower panel) suggest the inception of glaciation (e.g., Joussaume and Braconnot, 1997), as these areas were receiving about 35-40 W m −2 less insolation than they are receiving now.The data in these plots is computed with steps of 5 days and 5 degrees of latitude.Multi-millennial insolation time series can also be plotted in a 3-D surface plot as a function of year since J2000 and day of year, at the selected latitude.Users select the start and end years for the time series.The data for these plots are computed for steps of 1 kyr and one day (for day of year).An example of the output is provided in Fig. 3b.
Several time-series line plots are also produced.Insolation time series are plotted for the currently selected latitude; both the currently selected date and the annual average are shown (Fig. 4a).A multi-panel plot (Fig. 4b) allows comparison of the three Milankovitch parameters.Precession is visualized as the longitude of perihelion, as well as the climatic precession parameter, e sin ω (Berger and Loutre, 1994).A separate GUI button allows users to optionally produce timeseries plots of several paleoclimatic data sets (Fig. 4c).The top panel shows the EPICA CO 2 (Lüthi et al., 2008a, b) and deuterium temperature (Jouzel et al., 2007a, b) time series which go back to ∼ 800 kyr before present.The bottom panel of Fig. 4c shows two benthic oxygen isotope (δ 18 O) data  2001) data (Zachos et al., 2008).These data sets go back to 5320 kyr and 67 000 kyr before present, respectively.To first order, higher δ 18 O values are associated with higher continental ice sheet volumes and lower benthic ocean water temperatures (Zachos et al., 2001).For this reason, the y axis of the lower panel of Fig. 4c is inverted, so that higher values of EPICA CO 2 and temperature (generally warmer climates) from the upper panel of Fig. 4c can be easily associated with lower δ 18 O values (also generally warmer climates).These paleoclimatic data are included for convenience of the user and no further interpretation or analyses are provided.Users are cautioned that the interpretation of these paleoclimatic signals and their uncertainties, time resolution and chronology (age models) is fairly complex (e.g., Bradley, 2014, and data source references) and beyond the scope of this work.They are provided here for illustrative purposes only, e.g., this enables users to easily visualize the last few glacial-interglacial cycles (and the mid-Pleistocene transition to 100 kyr cyclicity, see Introduction), or to visually correlate these paleoclimatic time series with the corresponding Milankovitch parameter and insolation curves.

Numerical/Ancillary output
Ancillary data (and their units) are output in the main GUI window (Fig. 1) and are updated every time the Earth orbit plot (Fig. 2) is re-drawn (Sect.4.1) (i.e., every time the "Plot/Update Orbit" button is pressed).Variables that are output in the main GUI are as follows: solar declination, insolation at the TOA for the chosen date and latitude, day length, Sun-Earth distance, length of seasons (as defined in the North Hemisphere (NH)), the longitude of perigee, and true and mean longitude of the Sun.As a reminder, the longitude of perigee is the angle between the directions of vernal equinox and perihelion and true longitude is the angle Earth has swept from its orbit, subtended at the Sun, since it was last at vernal equinox; mean longitude is proportional to time instead (for detailed definitions, see Sects. 2 and 2.1 above).Users also are given the option of saving the data used to make the insolation plots in Fig. 3  for 21 March and less than 2 W m −2 for 20 June, respectively (Fig. 5a and b, solid lines with dots), which corresponds to less than 0.5 % of the absolute values (Fig. 5c and d, solid lines with dots).Importantly, these differences are generally much smaller or of the same order of magnitude as the corresponding differences between the Be78 and La2004 astronomical solutions as computed by Earth Orbit v2.1 (Fig. 5, dashed lines).Furthermore, these differences are generally smaller than the uncertainty resulting from varying estimates of the TSI (e.g., Fröhlich, 2013, vs. Kopp andLean, 2011, see Sect.2); also, these differences are smaller than the total contemporary anthropogenic radiative forcing on climate due to fossil fuel emissions (IPCC, 2013;their Fig. SPM. 5).
The Earth Orbit v2.1 model uses its own internally constructed orbital geometry and first principles equations to compute insolation.There is no additional a priori prescribed constraint to the model other than the orbital elements astronomical solution and the semi-major axis and orbital period (Sects. 2 and 2.3; Table 1).Therefore the validation presented here is an independent verification of the model's geometry and computations, taking the Laskar (2014) values as truth.Section 6 discusses sources of model uncertainty which can explain some of the small differences observed.

Solar declination validation and season length validation
Solar declination was validated against the algorithms of Meeus (1998).The model year is neither leap, nor common (Table 1) and is thus not equivalent to any single Gregorian calendar year.In order to validate declination at all dates, the Meeus (1998) algorithm was used to compute solar declinations for 12:00 UT on each date of four years (2009)(2010)(2011)(2012)2012 being leap) and average the declinations for each date (not day of year, Fig. 6).These averages were then compared to the solar declination output by the model for that date.Results indicate differences are always less than ∼ 0.2 • (Fig. 6, black line).By construction, model solar declination on 20 March will always be exactly zero degrees.
In reality, the exact instance of vernal equinox varies year to year, so these validation differences are expected.Importantly, the differences between the model and the 4-year averaged Meeus (1998) declinations are consistently smaller than the daily rate of change of declination (Fig. 6, green curve), as computed from the Meeus (1998) data.Additionally, these differences are of a similar magnitude to the standard deviation of declination between these four years for each date (Fig. 6, red curve).

Sources of uncertainties
Assumptions and approximations in the model and the underlying astronomical solutions propagate to uncertainties in the model outputs, such as declination and insolation.Some of these assumptions were already discussed, such as calculating insolation for a given calendar date vs. true longitude interval (fixed-date vs. fixed-angle calendars), and choosing integration steps for insolation time series (Sect.2.3.1).The calendar bias discussed in detail in Sect.2.3.1 means that if one compares insolation over geologic time on a given classical calendar date (e.g., 16 September), which occurs a given number of 24 h periods after the fixed vernal equinox, one is not necessarily comparing insolation at the same true longitude.The same argument is valid for an arbitrary interval of time longer than a day and shorter than a full orbital cycle.This calendar bias creates the artificial north-south tilt observed in insolation anomalies (Chen et al., 2011), which is also exhibited by the Earth Orbit v2.1 model output (Fig. 3a, second panel).This is expected because Earth Orbit v2.1 uses the classical calendar dates, which are more user friendly.
Next, we draw the users' attention to a few additional sources of uncertainty.Determination of some of these uncertainties is outside the scope of this work; however, users can run sensitivity analyses using the model in order to quantify them.Importantly, uncertainties in the astronomical solutions that are used as input to the model will propagate to insolation computations.There are differences between the different astronomical solutions (e.g., Fig. 5).Accuracy is highest near the present time and degrades further into the past or future (Laskar, 1999;Laskar et al., 2004).Chaotic components of planetary orbital motions introduce an uncertainty that increases by an order of magnitude every ten million years, making it impossible to obtain astronomical solutions for the Milankovitch parameters over period longer than a few tens of millions of years (Laskar et al., 2004).As a reminder, the Be78 solution is valid for one million years in the past or future, whereas the La2004 solution is valid from 101 million years before present to 21 million years in the future; however, solutions for times further back in time than 50 million 6. Solar declination validation: difference between solar declination as computed by the internal geometry of the Earth orbit model (using the La2004 orbital parameters for J2000) and mean actual declination from the years 2009, 2010, 2011 and 2012 as computed for 12:00:00 UT every day with the algorithms in Meeus (1998) (black solid line).The rate of change of declination (green solid line) and the standard deviation of declination for each date for the four years (red solid line, N = 4 for each data point) are also shown for reference.The model computations were performed with the calendar start date fixed at vernal equinox of 20 March.The date 29 February 2012 was removed from the analysis, so the abscissa corresponds to a given date, that is, dates, not days of year were averaged for a given mean solar declination across the four years.Abscissa ticks represent the 15th of each month.If 00:00 UT is used for the Meeus computations instead, differences (black curve) have a different pattern and are larger, but never exceed ∼ 0.4 degrees (not shown).See Sect.5.2 for details. years before present should be treated with caution (Laskar et al., 2004).
Due to the gravitational interaction of Earth and other solar system bodies, in particular Jupiter, Venus and the Moon, high-frequency variability (timescales of years to centuries) of the Milankovitch parameters is superimposed on the longterm low-frequency Milankovitch cycles.An example of such variability is the nutation in obliquity with a period of ∼ 18 years.These high-frequency fluctuations also lead to insolation changes.Bertrand et al. (2002) used results from the VSOP82 planetary position solution (Bretagnon, 1982) and a simple climate model to demonstrate that the amplitudes of these high-frequency variations and the effect on insolation and surface temperature is negligible (equivalent to model noise) as compared to the 11-year Sun cycle or the low-frequency trends.
The model prescribes the sidereal year as the orbital period (Table 1), which is slightly longer than the tropical year (Meeus, 1998).The difference is on the order of 0.01 days.The use of these two different period definitions leads to negligible differences in solar declination on a given date (for the J2000 La2004 orbital parameters; not shown), much smaller than the validation differences of Fig. 6.We conclude that the choice of orbital period does not influence the insolation computations significantly.
A single value for solar declination and the radius vector length is used in the computation of daily insolation (Sect.2.3).In reality, these quantities change continuously, instead of having discrete values.This is likely to introduce small errors in insolation that will generally be smaller in magnitude than the difference in daily insolation between successive days.Importantly, day length and daily insolation values near perihelion at very high eccentricities (that can occur only in Demo mode, not in the real Milankovitch cycles) should be treated with caution due to significant violation of the assumptions in applying Eqs. ( 6) and ( 7) (see Sect. 2.3).In such cases the radius vector length and the declination of the Sun may change significantly over the course of one 24 h period, and the hour angle of the Sun changes significantly more slowly than assumed (its time derivative is less than one); this is not dealt with rigorously in this implementation.
Sunrise and sunset times used in the insolation computation are referred to the center of the disk of the Sun and the mathematical horizon at the given latitude.Note also that irradiance is given at the top of the atmosphere (TOA), but all computations are geocentric, rather than topocentric, which should lead to negligible insolation differences.Since the model year is not an integral number of days, if total annual insolation is computed by summing daily insolation values, the 19 March insolation needs to be scaled by 1.256363 to reflect the fact that this day is 24 × 1.256363 h long in the model (Berger et al., 2010).Here, we average daily insolation to output average annual insolation, so this correction is not applied.

Concluding remarks
We presented Earth Orbit v2.1, an interactive 3-D analysis and visualization model of the Earth orbit, Milankovitch cycles, and insolation.The model is written and runs in MATLAB ® and is controlled from a single integrated userfriendly GUI.Users choose a real astronomical solution for the Milankovitch parameters or user-selected demo values.The model outputs a 3-D plot of Earth's orbital configuration (with pan-tilt-zoom capability), selected insolation time series, and numerical ancillary data.The model is intended for both research and educational use.We emphasize the pedagogical value of the model and envision some of its primary uses will be in the classroom.The user-friendly GUI makes the model very accessible to non-programmers.It is also accessible to non-experts and the primary and secondary education classroom, as minimal scientific background is required to use the model in an instructional setting.Disciplines for which the model can be used span mathematics (e.g., spherical geometry, linear algebra, curve and

Figure 2 .
Figure 2. (A) Present (J2000) orbital configuration for 16 September, using the La2004 solution and the calendar start date fixed at vernal equinox on 20 March.The orbital ellipse is shown in blue, the semi-major and semi-minor axes (perpendicular to each other) are in red and the lines connecting the solstices and equinoxes (also perpendicular to each other) are shown in black.The perihelion point, as well as the equinoxes and solstices are labeled (NH = Northern Hemisphere).The Sun is shown as a semi-transparent yellowish sphere centered at one of the orbital ellipse's foci, both of which are marked with an "x" along the semi-major axis.The Earth is plotted with its center on the corresponding place along the orbit, and the angle it has swept since last perihelion passage (the true anomaly angle) is filled in semitransparent light green.Earth's Equator is plotted as a solid black line, and its axis of rotation is plotted as a dotted black line , with the North Pole marked.The spheres of the Earth and the Sun are not to scale, the rest of the figure is geometrically/astronomically accurate and to scale.This plot is in 3-D and has pan-tilt-zoom capability in the Earth Orbit v2.1 model.The corresponding GUI with numerical ancillary output is shown in Fig. 1 (for latitude 43 • N and S o = 1366 W m −2 ).(B) Real orbital configuration for 16 September, 10 kyr in the future, using the La2004 solution and a 20 March equinox as calendar start date.(C) Demo (imaginary) orbital configuration for 1 July (vernal equinox fixed at 20 March), eccentricity = 0.6, obliquity = 45 • , longitude of perihelion = 225 • .The geometry is consistent with Berger et al. (2010), their Fig. 1, although it is being viewed in (A) and (C) from the direction of fall equinox, as opposed to from the direction of vernal equinox in their figure.The apparent eccentricity of the three orbits in Fig.2is also due to the view angle of the 3-D plot and the respective projection onto a 2-D monitor/paper; the intrinsic eccentricity can be judged by tilting the plot or observing the relative distance from the two foci (the Sun being at one of them) to the center of the ellipse, the intersection of the semi-major and semi-minor axes (red lines).

Figure 3 .
Figure 3. (A) A day of year-latitude insolation plot for 115 kyr before present (J2000) (upper panel) and the corresponding anomaly from J2000 (lower panel), using S o = 1366 W m −2 .(B) Insolation time series at 65 • N as a function of day of year, spanning 200 kyr before and after present (J2000).Negative years are in the past.Both (A) and (B) use the La2004 solution.
Figure 4. (A) Insolation time-series plot spanning 200 kyr before and after present (J2000) at 65 • N on 20 June (blue) and annual average (red); (B) time-series plots of Milankovitch orbital parameters spanning 500 kyr before and after present.Panels from top to bottom display eccentricity, obliquity, and longitude of perihelion and climatic precession; both (A) and (B) use the La2004 solution.(C) Time-series plots of paleoclimatic data spanning one million years before present: EPICA ice core CO 2 (blue) and deuterium temperature (green) (upper panel) and the Lisiecki and Raymo (2005) (blue) and Zachos et al. (2001) (red) compilations of benthic oxygen isotope (δ 18 O) data (lower panel).Note the y axis of the δ 18 O plot is inverted.Negative years for all Fig. 4 panels are in the past.

Figure 5 .
Figure 5. Absolute differences (W m −2 , solid lines with dots) between our insolation solution (using the La2004 astronomical parameters) and the Laskar (2014) insolation solution (also using the La2004 astronomical solutions; insolation provided by his Windows pre-compiled package at http://www.imcce.fr/Equipes/ASD/insola/earth/binaries/index.html) for 21 March (A) and 20 June (B).Differences between the Be78 and La2004 astronomical solutions (insolation for both computed by our model) are shown for comparison with dotted lines.Data are shown for three different latitudes -20 • S (red), 45 • N (green), and 65 • N (blue).(C) same as in (A) but displaying percent insolation difference, (D) same as in (B) but displaying percent insolation difference.Earth Orbit v2.1 insolation computations use the model's own orbital geometry with no additional a priori input other than the Milankovitch parameter solutions of La2004.Negative years are in the past.See Sect.5.1 for details.

Table 1 .
Summary of constant and variable model input parameters.
Sect. 2.3.1.The effect of the different choice of calendar start date is most apparent at exaggerated eccentricities and/or at longitudes of perihelion that are very different from the contemporary value.Insolation time-series output (Sect.4) is only computed for the calendar being fixed to vernal equinox on 20 March.
www.geosci-model-dev.net/7/1051/2014/Geosci.Model Dev., 7, 1051-1068, 2014 in ASCII format.The first row and column of these files list the abscissa and ordinate values of the data, respectively.Daily insolation is the most important model output from climate science perspective and is the fundamental discrete time unit at which the model calculates energy receipt at the TOA.Daily insolation was validated against the results of Laskar et al. (2004), as provided in Laskar (2014) (specifically, the precompiled Windows package at http://www.imcce.fr/Equipes/ASD/insola/earth/binaries/index.html).In both the Earth orbit model and the Laskar software, the La2004 solution for the orbital parameters was used, and the default model solar constant (Table1) was used.