New Long-Term Climate Oscillations

The astronomical theory of climate change is based on the solution of differential equations describing Earth’s orbital and rotational motions. The equations are used to calculate the change in insolation over the Earth’s surface. As a result of the author’s solution of the orbital problem, the periods and amplitudes of Earth-orbit variations and their evolution have been refined. Unlike previous studies, the equations of Earth’s rotational motion are solved completely. The Earth’s rotational axis precesses relative to a direction different from the direction of the orbit’s axial precession, and oscillates with periods of half a month, half a year and 18.6 years. Also, its oscillations occur with irregular periods of several tens of thousands of years and more. All these motions lead to oscillations of the obliquity in the range of 14.7˚ to 32.1˚, which prove to be 7 - 8 times larger than obtained by a previous theory. In the same proportion, the Earth’s insolation oscillations increase in amplitude, with insolation extremes occurring in other epochs than those in the previous theory. The amplitudes and the onset times of the extremes correlate with known paleoclimate changes. Thirteen insolation periods of paleoclimate variation over an interval of 200 thousand years of the orbit’s ascending node; i—orbit’s inclination; φ р —perihelion angle; Т—time in million years reckoned from December 30, 1949; Т е1 , Т Ω , and Т i —the least oscillation periods for the eccentricity, the ascending node, and the orbit inclination expressed in a thousand years; Т p —perihelion rotation period averaged over a


INTRODUCTION
Long-term climate oscillations are analyzed in the Astronomical theory of climate change (or alternately named, Astronomical theory of ice ages). The theory is based on the solution of the following three problems: 1) what are the changes in the Earth's orbit ? 2) what are the changes in the Earth's axis of rotation? 3) what are the changes in the amount of solar radiation over the Earth's latitude based on the first two changes? The original version of the theory was developed by M. Milankovitch [1] in the first quarter of the 20 th century. Subsequently, the proposed approach was improved by other researchers [2][3][4][5][6]. At the end of the 20 th century, activities aimed at revisiting the above problems were initiated [7]. As the result of a more precise solution to the problem of Earth's rotational motion, a second version of the climate 3. RESULTS

Earth's Motions and Their Variations
The Earth moves in an elliptical orbit around the Sun, which is located at one focus of an ellipse ( Figure 1). The shortest Earth-Sun distance in perihelion is denoted by R p , and the largest distance in aphelion by R a . The period of Earth's motion with respect to motionless space connected with the Solar system is P sd = 365.25636042 days [16,17]. The quantity P sd is called the sidereal period of the Earth's rotation around the Sun. The Earth's orbital motion proceeds in an anticlockwise direction, based on the orbit being viewed from the Earth's North Pole N. The normal to the orbital plane is denoted as S  , and is called the orbital axis.
With respect to motionless space, the Earth rotates around its axis, N  at an angular velocity of ω E = 7.292115•10 −5 1/sec in an anticlockwise direction coincident with the direction of the Earth's orbital motion. The value of ω E corresponds to a full revolution performed by the Earth in 0.99726968 days. The Earth's rotational axis N  is inclined to the orbital axis S  at an angle equal in the contemporary epoch, to ε = 23.43˚. This inclination is called the obliquity. During the orbital motion of Earth, the orientation of its rotational axis N  remains unchanged in space ( Figure 1). That is why at two points of the orbit, at times March 20, (20.03), and September 22, (22.09), the axis N  is normal to the Earth-Sun direction. With respect to Earth, the Sun is in its equatorial plane. That is why the southern and northern hemispheres receive identical amounts of solar radiation and the day appears to be equal in duration to night. These points are called the day of vernal equinox, (20.03), and the day of autumnal equinox, (22.09). At the time, June 21, (21.06), the axis N  is least inclined to the Earth-Sun line, and the northern hemisphere is therefore illuminated with maximum solar radiation. At the time of December 21,(21.12), the axis, N  Natural Science  [14,15,18,19].
is at maximum inclination to the Earth-Sun line. This is the reason the southern hemisphere is at maximum illumination at that time, and polar night approaches at high latitudes in the northern hemisphere.
Since the time spent on reaching and leaving the extreme angles lasts for several days, these points are called respectively, the summer solstice day (21.06), and the winter solstice day (21.12).
The inclination of Earth's axis, N  to the orbital axis, S  leads to the variation in duration of sunlight, both during the year and on the same day at different latitudes. On summer solstice day, (Figure 1, 21.06), we have a polar day in the whole region between the North Pole and the Arctic Circle. Then, as the latitude decreases, the day gets shorter, reaching a 12-hour duration at the equator, and we have a polar night established below the Antarctic Circle. Contrary to this, on the day of winter solstice, (21.12), in the territory between the North Pole and the Arctic Circle we have a polar night. Then the day starts increasing in duration. At the equator, the day lasts for 12 hours, and a polar day sets in below the Antarctic Circle. As we approach the equinoctial points of 20.03, and 22.09, the difference between the days in latitude decreases in value. The day's duration along all latitudes becomes identical at 12 hours.
As the Earth moves along its orbit, the alteration of seasons occurs. The duration of the seasons is defined by the Earth's motion over certain orbital segments. From the vernal equinox day, 20.03, until the summer solstice day, 21.06, the duration of spring is 92.7 days. Over the summer segment, the duration is 93.7 days. Over the autumnal segment, the duration is 89.9 days. Over the winter segment, the duration is 89.0 days.
The Earth's orbital and rotational motions define the variation in climate in the current epoch. However, the motions vary in time, and the climate undergoes change. The position of Earth's orbit precesses in space. The Earth's orbital axis S  (Figure 1), rotates, or in other words, it precesses about the direction of M  , which is motionless in space. The precession proceeds clockwise with a period of 68.7 thousand years. In a clockwise direction, the Earth's axis N  precesses about the direction of 2 M  , also motionless in space. The precession period is 25.74 thousand years. Besides this, the axes S  and N  execute oscillations, each with respect to its own precessional axis M  and 2 M  , respectively. In addition, the shape of the orbit (its eccentricity, whose value varies from 0 to 0.064; the current value being equal to 0.016) and the perihelion position, both undergo variations. Today, the perihelion is over the winter segment ( Figure   1), when winter sets in the northern hemisphere. Since the Earth's orbital perihelion rotates in anticlockwise direction at a mean period of 147 thousand years, its position in other epochs can be at any point in the Earth's orbit.
The changes of the Earth's orbital and rotational motions lead to changes in its climate. Below, the latter changes are considered in more detail. Natural Science

Evolution of the Earth's Orbital Motion
The evolution of the Earth's orbital motion is considered in a motionless frame, xyz whose origin is located at the center, O of celestial sphere 1 ( Figure 2). Note that, depending on the particular problem of interest, the point O can be located either at the center of mass of the Solar system, at the center of the Sun, or at the center of the Earth. As a result of the interaction of Solar-system bodies, the Earth's equatorial plane 2 and its orbital plane 3 both alter their positions, which are denoted with digits 2 and 3 in the epoch Since the annual motion of the Sun over celestial sphere 1 with respect to the Earth, proceeds along circles 5 or 3, the planes of the circles are additionally called, the planes of moving and motionless ecliptic, respectively. The frame xyz is related to the plane of the motionless Equator 2. The moving plane 5 of the Earth's orbit is defined by the angle φ Ω = γ 0 γ 2 of the ascending-node position γ 2 and by the inclination of the plane i.
The Earth moves around the Sun along an open trajectory which is close to the shape of an ellipse. At one point in the trajectory, the perihelion, the Earth approaches the Sun to the shortest distance R p , and at the opposite point, the aphelion, it moves away from the Sun to the largest distance R a . In Figure 2, the projection of the perihelion onto the celestial sphere 1 is denoted with the character B and its position is coordinated with the angle φ p = γ 2 B. The shape of the orbit is defined by its eccentricity The oscillations of the angles φ Ω and I reflect the rotation process of the orbit's axis S  (Figure 2), that is, the normal to the orbit plane, with a period of 68.7 thousand years about the motionless vector M  [20]. The latter vector is the sum of the angular momentum vectors of all bodies in the Solar system. The are manifested in the behavior of angles φ Ω and i in Figure 3.
Thus, the evolution of the Earth's orbit proceeds as a result of the following four motions: 1) precession of the orbit's axis S  ; 2) oscillations of the orbit's axis S  ; 3) oscillations of the orbit's eccentricity e; and 4) rotation of the orbit in its own plane (perihelion rotation). Figure 3 shows the evolution of the Earth's orbital parameters e, φ Ω , i and φ p over a span of one million years. The shortest eccentricity oscillation period is T e1 = 94.5 thousand years, and the longer ones are the periods T e2 = 413 thousand years and T e3 = 2.31 million years [19,20]. Both the angle, φ Ω of the orbit's ascending node and the orbit's inclination, i oscillate with a period T Ω = T i = 68.7 thousand years. The increasing behavior of the angle, φ p toward future ( Figure 3) reflects the non-uniform rotation of the perihelion in anticlockwise direction with a mean period of T p = 152 thousand years (over a time interval of one million years). As it is seen from the graph, there exists an epoch of a reverse, or clockwise, perihelion's motion.

Evolution of the Earth's Rotational Motion
Evolution of the Earth's rotational motion is treated in the motionless frame x e y e z e ( Figure 2) connected Figure 3. Evolution of Earth's orbital parameters over a time interval of 1 million years: е-eccentricity; φ Ω -angle of the orbit's ascending node; i-orbit's inclination; φ р -perihelion angle; Т-time in million years reckoned from December 30, 1949; Т е1 , Т Ω , and Т i -the least oscillation periods for the eccentricity, the ascending node, and the orbit inclination expressed in a thousand years; Т p -perihelion rotation period averaged over a time interval of 1 million years. Natural Science with the Earth orbit's motionless plane, 3. The inclination angle and the precession angle ψ = γ 0 γ 1 both define the position of equator, 4 moving with respect to the orbit's motionless plane, 3. The precession angle ψ decreases toward future in an oscillatory manner at a mean rate of  Figure 4. Over the latter time interval, the quantity, Δψ varies from −0.184 to 0.233 radian, so that the full oscillation swing amounts to 0.417 radian.
The inclination difference, Figure 4, is given with respect to the initial angle θ 0 = 0.40904645. The quantity Δθ oscillates similar to, Δψ, yet in a narrower range from −0.0845 to 0.0855, so that the oscillation range here amounts to 0.17 rad. Thus, the oscillation amplitude for angle θ is 2.45 times smaller than that for angle ψ. Besides, the oscillations of Δθ do not coincide in phase with the oscillations of Δψ; they are shifted along the time axis T by −7.5 thousand years.

Evolution of the Earth's Orbital Motion Relative to Its Rotational Motion
The orbital-motion parameters i, φ Ω and φ p and the rotational-motion parameters ψ and θ are defined by the obliquity, ε and perihelion angle φ pγ of the orbit's moving plane 5 with respect to the moving equator 4 ( Figure 2). The oscillation spectrum of φ pγ is rather broad since the angles φ p , φ Ω , i, ψ and θ contribute to the oscillations of this angle. The average variation of the angle, φ pγ proceeds according to the law φ pγm = φ p − (2π·T/P pr ). Figure 5 shows the variation of ε over five different time intervals In [13]. Over short time intervals, the oscillations of θ and ε are roughly identical. Indicated in the graphs are the main oscillation periods T ni and amplitudes (θ ai and ε a4 ) of the inclination angle: half-month period T n2 and half-year periods T n3 and T n4 = 18.6 years. Those oscillations are called nutation oscillations. The precession angle ψ exhibits similar oscillation periods, the amplitudes being two or three times greater.
Over the time interval In = 0.1 year, half-monthly oscillations are observed and diurnal oscillations with period T n1 = 0.9973 day can be traced; over the interval In = 1 year, half-year oscillations emerge; over the interval In = 10 years, an oscillatory trend with a T n4 = 18.6-year period is observed, with oscillations at this frequency prevailing over the time interval of In = 100 years.
Over the time interval In = 100 years, it is seen that the calculated obliquity ε 1) oscillates about its mean value, 2) received by S. Newcomb [16] and J. Simon et al. [17]. The oscillation amplitude ε a4 = 9.2" of the period, T n4 = 18.6 years also coincides with the observations. In astronomy, this quantity is called the nutation constant.

Evolution of the Earth's Obliquity and Insolation over a Span of 1 Million Years
As it is evident from Figure 5, over the time interval of In = 10 thousand years, a coincidence of the new obliquity ε 1 with the data 2 and 3 yielded by the first version of the Astronomical theory of climate change [1][2][3][4][5][6] is observed over a span of 2000 years.
It should be noted that these authors solved the problem in question over a large time interval: Sharaf and Budnikova [3] for 30 million years, and Laskar et al. [6] for an even longer period. We showed [18] that their results do not fundamentally differ from the results of Milankovitch [1], Berger and Loutre [4], Edvardsson et al. [5]. Therefore, we compare our results with these authors who are typical representatives of the former Astronomical theory.
As can be seen from Figure 5, after 2000 years, obliquity, ε calculated within the new version 1 of the theory shows clear deviation. As it is seen from the graphs of Figure 6, over the time interval of 1 million years the oscillations of ε as yielded by the second version of the theory proceed in the range of 14.7˚ to 32.1˚, whereas the same range in the previous theory was from 22.08˚ to 24.45˚; in other words, the range of oscillations in the second version of the theory proves to be seven times greater. This difference is due to the fact that in the second version of Astronomical theory, the Earth's rotation problem was treated in full, without simplifications. The solution of this problem and various checks of obtained data were analyzed at length in the publications [7,13,21] and will be covered in the subsection 4.3.
The amount of solar radiation reaching the Earth's surface, also called the Earth's insolation, is defined by the parameters e, ε and φ pγ . Figure 6 gives a comparison of the insolation changes 65N s Q occurring during the summer caloric half-year at the 65-deg northern latitude in the second version of the theory, (curve 1) [21] with the changes as calculated by the first version of the theory [6]. Here, the amplitude of insolation oscillations is also seven times greater than that in the previous theory. Besides, the insolation extremes occur at other times, and the oscillation periods are different. Note that the astronomical Natural Science summer and winter half-years measured from the vernal equinox day to the autumnal equinox day and visa versa differ in duration for different epochs. That is why it is caloric half-year, equal in duration, that are considered here.
In order to compare climates in other epochs with the current climate, we consider the insolations at equivalent latitudes I. For calculating of I, we consider the Earth's latitude φ characterized by receiving same amount of summer solar radiation, Q s as in the current epoch. Figure 6 shows the insolation oscillations at equivalent latitudes I over a time interval of 1 million years. The lowest values I ≈ 90˚ indicate that at latitude 65˚N in summertime, there was less solar radiation on the pole than now. The highest values, such as I ≈ 23˚ at the time −0.031 million years denote epochs in which, in summertime, the amount of solar radiation having reached the Earth at latitude 65˚N exceeds the amount of solar radiation having fallen onto it presently in the tropics, i.e. in the equatorial area. Such profound insolation oscillations lead to substantial climate oscillations. As it is seen from curve 2, the oscillations of I in the previous theory were less significant. Natural Science GJ/m 2 , respectively. In those epochs, the obliquities were ε = 23.44˚, 32.10˚ and 14.8˚, respectively. The summer insolation Q s (Figure 7, dashed lines) in contemporary epoch, 1 exhibits a minimum value on the poles, and it reaches a maximum value at the tropics φ = ε, and in the vicinity of the equator, this insolation exhibits a minimum. As we move from the cold (line 3) to the warm epoch 2, the summer insolation Q s on the poles increases by a factor of 2.07. At latitude 65˚N this insolation increases by a factor of 1.57. Since, on the average, this latitude well represents the variation of insolation at high latitudes, the above insolation was adopted by Milankovitch [1] as reference one in the climate characterization procedure. In warm epoch 2, the summer insolation Q s has an equatorial minimum in the southern hemisphere, and in cold epoch 3, in the northern hemisphere.

Variation of Insolation over the Earth's Latitude
The winter insolation Q w (Figure 7) on the poles is zero, and it monotonically increases in the equatorial region. In the equatorial region, the insolation Q w exhibits a maximum at a latitude, φ at which the summer insolation Q s shows a minimum. Over the period from cold epoch 3 to warm epoch 2, the winter insolation Q w exhibits most pronounced variations at middle latitudes. In the latter situation, for epochs 2 and 3 under consideration, e.g. at the latitude φ = 40˚, the change of the winter insolation is 1.38 times greater in the northern hemisphere than in the southern hemisphere. In cold epoch 3, the winter insolation at all latitudes is greater than that in warm epoch 2. In other words, during the cold epochs the winter seasons are warmer than those during the warm epochs.
The annular insolation Q T (Figure 7) monotonically increases from the poles toward the equator. At the equator, the annular insolation exhibits a maximum, with the annular insolation being symmetrical with respect to the equator. In other words, the amounts of heat per year are identical in both hemispheres. and at the latitude, φ = 45˚ the annular insolation experiences no changes. In the equatorial region, the changes of Q T are reciprocal to its changes at the high latitudes: in cold epoch 3, the amount of heat per year exceed that in the warm epoch. In the latter situation, the change of insolation Q T is four times smaller than that in the high-latitude region. That is why the main changes of the annular insolation occur at high latitudes.

Periods and Gradations of Earth's Climate Changes
Over the previous interval of 200 thousand years (see Figure 8), 13 climatic periods O I , 1 I , 2 I , 12 I were identified [19,22]. As a result of the comparison of these periods with paleoclimate data for Western Siberia over 50 thousand years, it was found that the periods 3 I , 2 I , 1 I , O I refer respectively to the Ermakov ice age, Karginsky warming, Sartan glaciation, and Holocene optimum. Those events also correspond to ice ages and interglacial periods in Europe and North America. Also, the following gradations of the warm and cold climate were introduced ( Figure 8): moderately warm, warm, and extremely warm climate levels, and moderately cold, cold, and extremely cold climate levels. During the past period of 1 million years (see Figure 9), the Earth has experienced six extremely cold (e.c.) periods and four extremely warm (e.w.) periods. The total number of cold (c.) and warm (w.) periods was 16 warm periods and 16 cold periods. Other periods were moderately cold (m.c.) ones and and 2c-the first and second boundaries of cold levels; m.w., w., e.w-moderately warm, warm, and extremely warm levels; m.c., c., e.c. -moderately cold, cold, and extremely cold levels.

Statements of the Problems and the Differences among Them
In the previous Astronomical theory of climate change, variations of Earth's orbital elements were received on the basis of the theory of secular perturbations. This is an approximate analytical method for solving the interaction problem for Solar-system bodies. In that theory, the changes of the equator plane 4 ( Figure 2) were analyzed approximately. With respect to this plane, the angle ε of the inclination of the orbital plane 5 and the perihelion position (being, in our notation, the angle φ pγ in Figure 2) were determined.
In the new Astronomical theory the interaction problem for Solar-system bodies (for simplicity, we will call this the orbital problem) was solved, using no simplifications, with a high-precision numerical method using a specially developed Galactica system [9][10][11]. In the latter case, we consider a change in the orbital plane (item 5 in Figure 2) relative to the fixed space represented by the equatorial plane 2 on a certain epoch. Here, the inclination, i of the orbital plane 5 and the ascending-node angle φ Ω differ from the angles figuring in the theory of secular perturbations. As it was noted above, in this theory the angle ε between the moving equatorial plane 4 and the moving orbital plane 5 was used.
Unlike in the theory of secular perturbations, we also numerically solved a second problem; the one about the Earth's rotation, governed by its own differential equations. As a result, we obtain the laws of variation of the inclination angle θ and the precession angle ψ of the moving equatorial plane 4 with respect to the motionless orbital plane 3 ( Figure 2).
A third complex geometric problem was then solved analytically for determining both the obliquity ε of the moving equatorial plane 4 with respect to the moving orbital plane 5 and the perihelion angle φ pγ .
The solution of the problem over short time intervals of the order of one thousand years is available in celestial mechanics. Over time intervals of several million years, during which those two planes and the orbit perihelion executed many irregular rotations in both directions, no solution was known.
A fourth problem, which yields the variation of the insolation versus the change of Earth's orbital and its axis parameters, thus offering an insolation theory, was presented in its complete form by M. Milankovitch in the first quarter of the twentieth century. We have also solved this problem in a new way [15]. Here, a new mathematical algorithm for elliptic motion, more understandable to non-specialists and useful for computer calculations, was employed.
All equations, including the differential equations for the orbital and rotational motions, were derived using a novel approach. Since the resulting data were other quantities in different representations, new methods for their analysis were developed. All that activity was accompanied with the development of computer codes written in various programming languages.

Adequacy of the Solution of the Orbital Problem
In connection with the new solutions of all four problems, at due stages, all the problems were subjected to checking. For checking the adequacy in the solution of the orbital problem, nine trustworthiness criteria were developed. Some of those criteria were included into the Galactica program; that is why control was exercised right in the course of solving the problem. For all bodies having an observation base (these are the planets from Mercury to Neptune, and also the Moon, the solutions over the time interval of several thousand years were compared with the secular variations of orbital parameters. The coincidence was found to be excellent [20,23].
Over intervals of hundred thousand and million years, the orbital parameters were compared to the results obtained by the previous researchers [3,4,6]. The data were also found to be coincidental. Each of subsequent authors took into account the experience gained by previous researchers and made the theory of secular perturbations more refined. The theory was confirmed by making necessary comparisons. The Natural Science later the works were published, the longer was the time interval over which the corresponding results were coincident with our data [20].
As it was noted previously, the perturbation theory is an approximate method for solving the orbital problem. After 20 million years, the solution obviously started departing from the actual data: the orbits of individual planets started increasing in size and, later, the planets could leave the Solar system [24]. We have solved the orbital problem over a time interval of 100 million years. All the orbital parameters of the planets and Moon executed steady oscillations, and there existed no tendency toward the change of those oscillations [20].

Adequacy of the Solution of the Earth's Rotation Problem
Points concerning the reliability of the solution of the Earth's rotation problem were analyzed in detail in publications [7,13,14]. Within the adopted solution method, all necessary checks were performed. For instance, the problem was solved in succession, implying the action of one of ten bodies (the Sun, the eight planets, and the Moon) [12]. The obtained oscillation periods of the Earth's axis were confirmed with general conclusions made on the basis of the theorem of change in angular momentum and also, with the results of other authors [25]. With the actions due to all bodies, the problem was solved over different time intervals; as it was shown previously in sec. 4 and sec. 5, the obtained data proved to be coincident with observations. Integration of the equations over a time interval of 200 thousand years was performed with different initial conditions and integration steps. The latter resulted in no changes of the oscillation periods, oscillation amplitudes, and the onset times of the extremes.
From the graph with the interval I n = 100 yr in Figure 5, one can see that the middle of the oscillations of the calculated 1 obliquity ε coincides with the observed average angle 2. That average angle 2 complied with all ancient observations made over a period of 2.5 thousand years of their history. The obtained oscillations with 9.2'' amplitude and 18.6-year period were found to be precisely coincident with the observed oscillations. In the graph of Figure 5, with the interval I n = 10,000 yr., it is seen that the deviation of the calculated angle, ε from the linear trend established for observational data, starts manifesting itself after 2.5 thousand years. In Figure 5, solutions for future time intervals are shown. The solutions for previous time intervals have a similar general appearance.
Following 2.5 thousand years, solutions obtained by previous researchers also started departing from the linear trend. From this time on, differences between solutions 2 and 3 by those authors and our solution 1 was observed. Over a longer time interval, namely, 200 thousand years, the Earth's rotation problem was solved initially towards the future [26]. The solution for, ε was found to feature a different oscillation structure and other onset times for extremes, and what was most important, the amplitude of new oscillations exceeding the amplitude of the previous oscillations by a factor of 7 -8 was revealed. From this time on, a solution check for the Earth's rotation problem was initiated. This check lasted for three years.
The Earth's rotation problem is one of the most complex problems in mechanics. Its solution can depend on the fundamental assumptions made while deriving the equations, on the choice of initial conditions, and on the procedure of solution reduction to the Earth's moving orbit. That is why a cardinal check of obtained results would be their obtaining, without solving differential equations for rotational motion.
While studying the orbits of the bodies, we found that the evolution of the Moon's orbital axis was similar to that of Earth's axis of rotation. That result had led us to a compound model for Earth rotation in which, part of the Earth's mass was uniformly distributed among peripheral bodies rotating around a central body, along a circular orbit. Under the action of the Moon, the Sun, and the planets, the orbits of the peripheral bodies started changing. It should be noted that under the axis of the orbit is meant a perpendicular to its plane. The evolution of the orbital axis of one of the bodies, modeled the evolution of the Earth's axis. Such modeling of the Earth's rotational motion included performing several solution stages of the orbital problem being treated with the Galactica program. In the initial series of our studies [20,27], three models were studied, and the possibility of modeling the evolution of the Earth's axis was confirmed. In those models, the precession periods of the orbit axes were 170 and 2604 years, whereas the average pe-Natural Science riod of the Earth's precession axis is known to be 25,740 years. Was subsequently yet developed 11 models, while has not been reached the required period of precession [7,13]. In the 13 th model, the orbital radius of the peripheral bodies was equal to the Earth's radius, the rotation period of the bodies was 0.142 hour, and the interaction between the model's bodies was amplified by a factor of 9.6 in comparison with the gravitational interaction. Thus, the bodies of the 13 th model rotated 170 times faster than the Earth does. For studying the evolution of such models, the Galactica program was developed; in this model, the possibility of changing the interaction between the model's bodies was provided.
The solution of the problem for the 13 th compound model of the Earth, over a time interval of 300 years [7,13], has yielded all characteristics of the dynamics of the Earth's axis, including the half-monthly and half-year oscillations of the angles ε and ψ and the oscillations of those angles with a 18.6-year period.
The amplitudes of those oscillations have also proved to be coincident with the solution results of the Earth's rotation problem. Such a coincidence in the results of the model problem with the results of the direct problem occurs once over a period of 3 thousand years. Further solution of the problem is hampered by the necessity to reduce the integration step to values for which the computing time turns out to be too long. Thus, over a time interval of 3000 years the compound model of the Earth has confirmed the obtained results of integration of the differential equation for Earth's rotation. The latter indicates that the assumptions and simplifications adopted in the derivation of the equations, the derivation of the equations itself, the solution method, and the transforming of the integrated data into the final form have also been confirmed.
A second check consisted in using an alternative integration method. In the program DfEqAl1.for solving the Earth's rotation problem, a fourth-order Runge-Kutta integration method as implemented by Krut'ko et al. [28] was used. Over a time interval of 200 thousand years, a growth of daily oscillations of the derivatives ψ and θ  was identified. Next, a code for solving the DfEqADP8.for program was developed, which uses the Dormand-Prince method, i.e. the eight-order Runge-Kutta integration method [29]. In integrating the equations of rotational motion over a period of 200 thousand years all previously obtained results have found confirmation. In the latter case, the amplitude of the daily oscillations of ψ and θ  showed no increase and remained at one and the same level. Thus, the particular method for integrating the equations does not affect the obtained results, and the application of a more accurate method confirms them.
A third check consisted in employing another method for solving the problem. The differential equations for rotational motion contained the coordinates of ten bodies acting on the Earth (eight planets, the Sun, and the Moon). Those coordinates were determined while solving the orbital problem with the Galactica program. However, the data array over large time intervals obtained with an integration step equal to 10 −5 -10 −4 years took an overly large memory space. That is why we have developed a mathematical model for the Solar system [30] that yielded the coordinates of the bodies at a desired time, calculated by the formulas for elliptic motion in which the ellipse parameters at each time, were determined from the data initially calculated by the Galactica program. In the course of the solution process of the problem, the mathematical model for the Solar system was subjected to a thorough check. Nonetheless, there still existed a probability that, over large time intervals, the insignificant differences between the results of the mathematical model for the Solar system and the coordinate values obtained by the Galactica problem, could affect the evolution of the rotational parameters ε and ψ. We have developed a new program, glc3rte2.for, for joint solution of the orbital problem and the Earth's rotation problem. In the latter program, in one step over time, the Galactica program solved the orbital problem and, then, the Dormand-Prince method was used to solve, in the same step, the Earth's rotation problem. With the help of the new program, we have obtained solutions of those problems over different time intervals, including the interval of 200 thousand years. All previously obtained results have found confirmation. The latter check also confirms the mathematical model for the Solar system over large time intervals. Table 1 gives a quantitative comparison of precession periods P prN , and the minimum and maximum obliquities (ε min and ε max , respectively), accurate to five significant digits. Within this accuracy level, the first two methods have yielded identical results. As it is evident from Table 1, the data calculated by the Table 1. Comparison of results obtained by three methods for integration of rotational-motion equations over a period of 200 thousand years: RK-4-Runge-Kutta method of the fourth order; DP-8-Runge-Kutta method of the eighth order in Dormand-Prince realization; Gal-the bodies' coordinates are determined by the Galactica program, and the rotational-motion equations are solved by the DP-8 method. third method, differed in terms of precession period P prN in the fourth digit, and in terms of ε in the fifth digit. Since the latter method was a most accurate one, the latter values provided refined results in comparison with the data obtained by the first two methods. Thus, the various tests and verification of the initial solution method of the Earth's rotation problem and also, independent solutions of the same problem by three other methods, have confirmed the fact that the Earth's rotational axis executes oscillations at an amplitude of 7 -8 times greater than that obtained in the previous solutions.

Adequacy of the Solution of the Third and Forth Problems
As it was noted previously, the third problem on the determination of the obliquity ε between the moving equatorial plane 4 and the moving orbital plane 5, and the perihelion angle φ pγ (Figure 2), was a geometrically complex problem because a multitude of revolutions on the axes N  and S  , and the perihelion B, were executing about different axes in different directions. For instance, the Earth axis, N  executed 777 revolutions per 20 million years of problem solution. In the problem, it was required to determine, from the obtained solution of the orbital problem for i, φ Ω , and φ р (Figure 2) and the Earth's rotation problem for θ and ψ, the relative-motion characteristics, ε and φ pγ . Those transformations were derived; yet, they involved inverse trigonometric functions, known to be many-valued. The expressions themselves are cumbersome, with some their parts presenting imaginary or infinitely growing expressions. Those singularities needed to be identified for the revealing factors underlying their behavior and for developing algorithms for elimination of these factors. Initially, this problem was solved by means of spherical geometry. However, because of the complexity of involved logical concepts, there was no firm belief in adequacy of the obtained solution. Fortunately, an idea suggesting a second method had emerged. The axial vectors, N  and S  were replaced with their projections onto the axes of a Cartesian coordinate system. Then, trigonometric means were employed to derive the required transformations. The two transformation systems have allowed us to reveal errors in each of the systems and fix them, unless both systems started yielding the same results in the examined 20-million-year time interval.
The new algorithm developed for calculating the insolation in the fourth problem was checked [19] by performing insolation calculations using the orbital data by Milankovitch [1], J. Laskar et al. [6]. The new algorithm has yielded results, the same as those yielded by the Milankovitch algorithm.

Physical Cause of the Difference between the New and Previous Theory
As it was noted previously, according to the new solutions for the angles θ and ψ, the Earth's rotation axis N  rotates around the vector  . After substitution of actual values, we have obtained a value P rl = −41.1 kyr, that is, the very oscillation period of the angles θ M2 and ψ M2 . In the previous Astronomical theory of climate change, the main oscillation period of the obliquity, ε was equal to 41 thousand years. As was noted previously, that theory was based on a simplified solution of the Earth's rotation problem. The obtained simplified solution has resulted in identical positions of the vectors 2 M  and M  ( Figure 2) and finally, in the obtained period of 41 thousand years.
In the new theory (Figure 2), the vectors 2 M  and M  of the precessional axes have different orientations. That is why the moments of forces by which all bodies act on the Earth exhibit a wide range of oscillations. As a result, the oscillation range of θ with respect to the motionless orbit, (3) also increases. In addition, the oscillation range of the angle, ε between the moving orbital axis S  and the moving axis N  increases as well. All in all, in the new theory, the oscillation amplitude of angle, ε turns out to be exceeding the same amplitude in the previous theory by a factor of 7 -8.
Thus, the oscillation period of obliquity ε in the previous theory of Earth-axis dynamics, was the period P rl of the relative precession of Earth's rotational axis N  and its orbit S  . That is why it is a fact that the precessional axes M  and 2 M  of the orbital axis S  and the Earth axis N  were assumed coincident, was the physical cause behind obtaining erroneous results in the previous theory.

Final Verification
The final check of the Astronomical theory of climate change consisted in the comparison of its results with paleoclimate. While analyzing data gained by geologists, geographers, and other specialists in the field of paleoclimate, we have found [22,31,32] four extremes of summer insolation 65N s Q having occurred over the past 50 thousand years; namely, 4.16, 15.88, 31.28, and 46.44 thousand years ago; those dates correspond to the middle of Holocene, to the middle of the last ice age, to the middle of a warm period, and to the middle of the penultimate ice maximum, respectively. Those events are called differently in different regions of the world; yet, all of them have left their traces in Siberia, Europe, and Northern America.
The whole complex of performed studies and their tests outlined here give us grounds to assert that, in the present article, results of Astronomical theory of climate change obtained by taking into account all studies performed during past centuries, are reported.
Changes of obliquity, ε and angle of the perihelion, φ pγ , as well as the eccentricity of the orbit, e are available at http://www.ikz.ru/~smulski//Data/Insol/ in the file OrAl1c_8.prn for ±100 years, in the file OrAl-5kyr.prn for 5 thousand years ago (ka), in the file OrAl-200ky.prn for 200 ka, in the file OrAl0-5My.prn for 5 Ma. The program, Ins12bdEn.mcd for analyzing these data, calculating the insolation of the Earth and the building of graphs is available on this site.

CONCLUSION AND FURTHER DEVELOPMENT OF THE ASTRONOMICAL THEORY OF CLIMATE CHANGE
As a result of the interaction of Solar-system bodies, evolution of the Earth's orbital and rotational motions proceeds; this evolution in turn gives rise to insolation oscillations being the cause of climate changes observed over time intervals of tens of thousands of years. The same interactions lead to the evolution of the Sun's motion about the center of mass of the Solar system [20,33] and also, to a change in the rotational motion of the Sun. Studies show [34], that the change in motions presents the cause of varia-Natural Science tions in solar activity. The radiation fluxes due to the Sun and its substance act on the Earth's upper shells, thus leading to circulation changes arising in its atmosphere and ocean. Very likely, these factors act as the cause of short-term climate variations occurring over periods of several ten and hundred years. Further development of the Astronomical theory of climate change will be related to the determination of those oscillations.