New Computing Algorithm of the Earth ' s Insolation

M. Milankovitch developed the theory of the Earth’s insolation for research of paleoclimate. The calculations of insolation by the simplified methods are also used in meteorology, city planning, climatology and other fields, for example, for studying habitability and detectability of extrasolar planets. The new algorithm for computing insolation is based on the results of the exact solution of the two-body problem. It differs from the previous method by simplicity and greater adaptability to computing. The coherent reasoning is provided for a new algorithm and its realization is given in the MathCad environment. The distributions of insolation in latitude are computed for a year and for caloric half-year in the contemporary epoch. The change of insolation is computed in equivalent latitudes for 200 thousand years. The results of computing are compared with the results of other researchers, and the accuracy of the approach is proven. The dynamics of insolation in the contemporary epoch is computed. The prospects of application of its results are shown to determine the causes of changes of natural processes that depend on insolation. The developed approach and the program are applicable for calculating insolation of other planets. The paper points out some specifications which will be needed when computing insolation of Mercury, Venus, Jupiter and Saturn.


Introduction
The amount of heat reaching the Earth from the Sun, i.e. insolation of the Earth, is an important factor of processes taking place on our planet.One of the most likely causes of geological changes of the Earth's history may be a change of its insolation.The effect of insolation on the Earth's climate of past eras is studied for nearly two centuries.These studies are based on a method for insolation calculation, which M. Milankovitch has completely presented in his works, for example, in Milankovitch (1939).It is quite complicated.Therefore a series of operations such as the insolation calculation in equivalent latitudes are usually omitted.Climate is determined by insolation for all year or for half-year.However, in some papers (Borisenkov et al., 1983;Edvardsson et al., 2002) climate changes are studied by insolation change on the summer solstice day.In other papers (Berger & Loutre, 1991;Loutre & Berger, 2000), climate changes are considered in terms of insolation dynamics in July.In Gasperini and Chierici (1997) the Earth's insolation is considered, taking into account only the Earth's axis oscillation with the period of 18.6 years.
Recently, in the literature sources, there has been a lack of understanding of this theory due to its complexity.For example, it is claimed in Bolshakov and Kapitsa (2011) that Milankovitch ignored the direct contribution of the eccentricity changes to insolation curve in his calculations.In fact, M. Milankovitch developed the theory for insolation calculation of the Earth, taking into account all factors, including the eccentricity.Therefore, the development of a simple algorithm to compute insolation would facilitate its understanding and appropriate application.
It is worth noting that except Milankovitch, other researchers contributed to the modern theory of insolation.Hargreaeves (1896) has developed a theory of the mean annual insolation depending on astronomical factors.He was not interested in seasonal insolation and other components.Garfinkel (1953) in his paper laid the theory of R. Hargreaeves (1896).He used it for research of the dependence of the mean annual insolation on latitude and obliquity.A. Berger and his colleagues have greatly contributed to the development and the explanation of the insolation theory.In (Berger et al., 1993), there are the insolation theory and the definitions of various terms, including various kinds of insolation: insolation, irradiance, irradiation, daily irradiation etc.The paper (Berger at al., 2010) gives a clearer view of insolation.The history of its development is analyzed, hypotheses are underlined in derivation of equations, a foundation of the solutions of elliptic integrals and the definition of astronomical terms, used for insolation calculation, are given.Nevertheless, the basic principles of insolation calculation remain as is, therefore the insolation theory is named as the Milankovitch theory.
Change of insolation on the Earth's surface and its dynamics are of interest to specialists in various fields of human activity and in different fields of science.For example, the calculations of insolation and its regulation have an important role for justification of building cities (Baharev & Orlova, 2006).Insolation calculations are used in various power researches (Harvey, 1995) and for studying habitability and detectability of extrasolar planets (Dobrovolskis, 2009).Evolution of insolation is calculated not only for the Earth but also for Mars (Laskar, Correia, Gastineau, & Robutel, 2004;Edvardsson & Karlsson, 2008).The dynamics of insolation is more and more often used to find correlations while studying the causes of changes of various factors in biology, medicine and meteorology.In these works the researchers have to develop their own approaches to insolation calculation of the Earth's surface and influence of motion of Venus (Ivanov, 2002), the other planets (Bogdanov, Katruschenko, & Surkov, 2006) and the Moon (Bogdanov & Katruschenko, 2008) on it.These approaches, as a rule, take into account only the change in the distance between the Earth and the Sun from all factors that affect insolation.For example, Bogdanov and Katruschenko (2008) derived the oscillation amplitude of insolation of 84.9 mW/m 2 due to monthly changes of this distance at the motion of the Moon.From the examples of insolation application, it is clear that the creation of a simple algorithm for its computing, taking into account all the factors which it depends on, is urgent.In the algorithm developed by Milankovitch (1939) for insolation calculation of the Earth, the differential equation is integrated to determine the dependence of the longitude λ of the Sun on time t by an approximate analytical method.It comes from the Kepler second law.We have derived (Smulsky, 2007;Smulsky, 1999) the exact analytical dependence of time on polar angle ϕ o or longitude λ by solving the two-body problem.This relation greatly simplifies the expression for insolation.Therefore, the method becomes clearer and the calculations are simplified.

The Main Results of the Two-Body Problem
To calculate insolation of the Earth, it is necessary to determine positions of the Sun above a point of the Earth's surface.Since the Sun is moving relative to the Earth in an orbit that is identical to the Earth's orbit, the Earth's orbital parameters are used to describe the movement of the Sun throughout one year.Similarly, the parameters of rotational motion of the Earth are used to calculate the daily movement of the Sun.
Figure 1.Schematics of motion of the Earth (E) in its orbit around the Sun (a) and the Sun (S) relative to the Earth (b): γ is the vernal equinox; P E is the Earth's perihelion; P S is perigee of the Sun, ϕ о is the polar angle of the Earth; ϕ pγ is the angle of the perihelion of Earth; ν = ∠γEP S is the perigee angle of the Sun; λ is longitude of the Sun; the dotted line and the its points E′, P′ E and γ E are conditional images of the Earth's orbit on the schematic of the orbital motion of the Sun Figure 1 shows the schemes of motion of the Earth relative to the Sun and the Sun relative to the Earth.The Earth's orbital plane (Figure 1a) crosses the plane of the equator by line Sγ.The Sun passes the vernal equinox γ in its orbit in the spring, and the Earthin the autumn.In point γ the Sun is in the plane of the Earth's equator, so the day duration is equal to the night duration.The planes of equator and orbit of the Earth change in space, so point γ is moved on the Earth's orbit (Figure 1).Angle ϕ pγ is measured between these two points γ and P E .
Let's consider the motion of a planet with mass m at gravitational Newton's law interaction with the Sun, which has a mass M. The differential equations of its motion relative to the Sun in the dimensionless variables are given by Smulsky (1999): where is the dimensionless radius of planet position relative to the Sun; ( ) is the trajectory parameter; μ 1 = −G(M+m) is the interaction parameter; G is the gravitational constant; P R is the perihelion radius; v p is the planet velocity at perihelion.
The movement in dimensionless variables r , t , as follows from (1), is completely determined by the trajectory parameter α 1 .It is worth noting that this parameter is identical to the eccentricity e, and they are related by the following expression: e = -(1 + α 1 )/α 1 .In the differential equations of the electromagnetic interaction (Smulsky, 1999), parameter α 1 plays a similar role of the parameter of trajectories to which the notion of "eccentricity" is inapplicable.Therefore, it is preferable to use α 1 instead of eccentricity e.
As a result of integration of equation ( 1), the equation for trajectory and relations for motion time t(r) from the distance r between the bodies for various cases of motion are derived (Smulsky, 1999).The trajectory equation in polar coordinates (r, ϕ о ) has the following view in dimensional variables ( ) where ϕ о is the polar angle of the body on the orbit, measured from the perihelion radius r = SP E = P R (see Figure 1а).Equation (2) represents: a circle at α 1 = -1, an ellipse at -1 < α 1 < -0.5, a parabola at α 1 = -0.5, a hyperbola at -0.5 < α 1 < 0, a straight line at α 1 = 0.It should be noted that the trajectory parameter α 1 can be determined at a known orbit through the smallest P R and the biggest a R orbit radii as follows: ( ) This equation shows that P R = a R and parameter α 1 = -1 for the circle, and and α 1 = -0.5 for the parabola.
According to the solutions of the two-body problem (Smulsky, 2007), the velocity at perihelion can be calculated by the equation: where P sd is the orbital period.In the two-body problem, it is considered in relation to a non-accelerated coordinate system.It's the period of its revolution around the Sun relative to the fixed stars (sidereal period) for a planet.The average period P sd = 365.25636042days for the Earth.This value of P sd is derived as a result of generalization of the observation data for all the history of astronomical observations.While computing the interaction of the planets, the Sun and the Moon, it was received (Melnikov & Smulsky, 2009;Smulsky & Smulsky, 2012) that in different epochs the orbital periods change little about P sd .Therefore, this difference is neglected.
To find the dependence of time on a polar angle, the function r(ϕ o ) is substituted in the solution t(r) of the two-body problem (Smulsky, 2007) according to (2).Then the time of a body's movement on its orbit (see Figure 1a) from the perihelion point P E to point E with angle ϕ o is calculated by equations: Here, as follows from (4), the motion time from the perihelion to aphelion at After substituting v p from (3) to (5), it is derived that t a = 0.5⋅P sd .In the solutions of the full orbital task and in observations, the motion time from perihelion to aphelion t a varies little about 0.5⋅P sd .In this formulation we neglect these oscillations.Let's consider the neglect of oscillations P sd and t a as the first simplification.

Geometric Characteristics of Insolation
Let's consider solar irradiation of point M on the Earth's surface (see Figure 2).Let's use all basic notations accepted by Milankovitch (1939).The plane of the horizon for point М is displayed on the celestial sphere 1 by the horizontal circle HH′.A perpendicular to plane HH′ intersects the plane of the celestial sphere 1 in the point of zenith Z.The Sun S performs annual motion on its orbit around the Earth, which is projected onto the celestial sphere by the ecliptic circle EE′.The movement occurs counterclockwise with the origin of longitude λ of the Sun at the vernal equinox γ.At this point the Sun is in the equatorial plane, when it passes from the southern hemisphere to the northern hemisphere.
The daily movement of the Sun occurs as follows.The Earth together with observer M, horizon circle HH′ and meridian NZE′A′H′ rotates around rotation axis NM, which is called a world axis.The rotation occurs counterclockwise.Therefore, the Sun moves clockwise in circle SrMdSd in parallel to equator AA′ relative to the Earth and, in particular, relative to horizon circle HH′.It rises above horizon HH′ at point Sr, it is at Md at noon, and it goes down over the horizon at point Sd.The hour angle of the Sun ω will be measured from the noon point Md.The arcs from sunrise to noon SrMd and from noon to sunset MdSd have the same length, which is designated as ω 0 .Therefore, the hour angle of the day varies within the limits of -ω 0 ≤ ω ≤ ω 0 .The duration of the day is equal to ω d = 2ω 0 , and the duration of the night is ω n = 24 -ω d .Here we used the hour angle ω in hours.In addition, it will apply hereinafter in angular units.
Figure 2. The basic geometrical characteristics of the Sun S at irradiation of point M on the Earth's surface: 1 is the celestial sphere; НН′ is the plane of the horizon; N is North pole; AA′ is the plane of the mobile equator; ЕЕ′ is the plane of the mobile ecliptic, and ε is the angle between planes AA′ and ЕЕ′; Z is zenith of point M, and z = ∠ ZMS is the zenithal angle of the Sun; the arc HN = ϕ is the geographical latitude of point M; ω = ∠ ZNS is the hour angle of the Sun, measured from noon; δ = SB is the declination of the Sun; λ = γS is the longitude of the Sun At the position of observer M and the Sun S, as shown in Figure 2, the duration of the day ω d is longer than the duration of the night ω n .When the Sun S is at E′, the duration of the day will be the longest.This is the point of the summer solstice.When the Sun S is at point E, the duration of the night is the longest.This is the point of the winter solstice.And when the Sun S is at γ or γ′, its daily motion will occur by the equator's circle AA′.This circle intersects the horizon's circle HH′ in its diameter, so the passage of the Sun over the horizon and below it is the same, i.e. the duration of the day is equal to the duration of the night.
If the observer at point M (see Figure 2) is located at higher latitudes, i.e. arc HN will be higher, then circle SMd will not cross the horizon circle HH′.In this case observer M will have a polar day.When the Sun S is in the southern hemisphere close to point E, the circle of its daily motion will not cross the line of horizon HH′.In this example of the latitude, the observer at point M will have a polar night.
The plane of equator AA′ and the plane of ecliptic EE′ change in space, so that the vernal equinox γ moves clockwise per year for 50″.25641.As already noted, the annual motion of the Sun passes counter-clockwise on the circle of ecliptic EE′, that is reflected by change of longitude λ, beginning from the vernal equinox γ.Therefore, the time of the Sun's passage of two consecutive points of the vernal equinox, i.e. tropical year P tr = 365.24219879days, is the less than sidereal year P sd .To avoid shifting the seasons of year in dates, our calendar is based on the tropical year.Therefore, the seasonal and half-year insolation of the Earth is determined by the duration of the tropical year.

Flux of Solar Heat
By neglecting the absorption of solar radiation in interplanetary space we can accept that the same amount of heat per unit of time proceeds through heliocentric spherical surfaces with radii r and a where dW n (r)/dt is the flux of radiant energy per unit area at the distance r from the Sun per unit of time; a is the average radius of the Earth's orbit, i.e. its semi-major axis a = 0.5 (R p + R a ), where R a is the aphelion radius.
The flux of heat from the Sun at distance a is called the solar constant J 0 M. Milankovitch used its value of J 0 = 2 cal/(cm 2 ⋅min).This value can be written in other units as: J 0 = 83.736kJ/(m 2 ⋅min) = 1395.6W/m 2 .
At present, the solar constant is known with a greater accuracy.Since 1907, the radiation of the Sun has been investigated by the Davos Physical Meteorological Observatory in Switzerland (Fröhlich & Finsterle, 2005), which has been operating under the aegis of the World Meteorological Organization since 1971.In 1996, J = 1366.784W/m 2 was accepted as the World Radiometric Reference (WRR) on the basis of data analysis from different radiometers.Since 1978, the radiation of the Sun has been measured outside of the atmosphere, i.e. on the satellites.Both in ground-based measurements, and in satellite ones the flux of solar heat in each series changes within the tenth parts of a percent.And in all series the radiation of the Sun ranges from 1358 W/m 2 to 1375 W/m 2 (Willson and Mordvinov 2003).The value of the solar radiation J = 1366.22W/m 2 is accepted as the Space Absolute Radiometric Reference (SARR) (Crommelynck, Fichot, Lee, & Romero, 1995).In order to maintain continuity with the calculations of other researchers in all calculations below, we use J 0 = 83.736kJ/(m 2 ⋅min), accepted by M. Milankovitch.
Based on (7), the flux of solar radiation according to expression (6) can be written as dW n (r)/dt = J 0 ⋅a 2 /r 2 .Then the flux of solar heat at distance r from the Sun on the Earth's surface area, perpendicular to the rays of the Sun, would be: where ρ = r/a is the relative distance to the Sun.
The line of zenith MZ (see Figure 2) is perpendicular to the considered area.If the angle between the Sun and the zenith is z, then the heat flux at point M is The zenithal angle z of the Sun depends on the angular coordinates of the points M and S. In the spherical triangle NZS, angle ∠ N = ω and the two adjacent sides NZ = π/2 -ϕ and NS =π/2 -δ are known.The side ZS of the triangle is an arc, which is measured by the central angle z (see Figure 2), i.e. z = ZS.This angle and the arc are determined by the law of cosines (Duboshin, 1976) After substitution ( 8) and ( 10) in ( 9), we receive the radiation law of point M depending on its latitude ϕ and angles δ and ω of the Sun's position Further, it is accepted that points M of the Earth's surface have the same relative distance ρ to the Sun at all latitudes ϕ.In reality this distance varies within value of the Earth's radius R E .The relative error is accepted about R E /a = 4.3⋅10 -5 for the calculation of insolation.Let's consider it as simplification 2. The third simplification concerns the form of the Earth.The zenithal angle z of the Sun is determined in relation to the line of zenith ZM (see Figure 2), which is drawn perpendicular to the surface of the Earth and passes through its center.For an ellipsoidal surface of the Earth it is true for points M at the poles and the equator.And at intermediate latitudes, the zenithal angle for the ellipsoid will differ slightly from expression (10).The fourth simplification is due to the fact that the solar constant J 0 is accepted unchanged.As already noted, it has some minor variations in time.Perhaps in the future, its variation law will be revealed which can be put into the calculation of insolation.At this stage the change of insolation, caused by these four factors, are not taken into account.

Daily Insolation
According to (10), the zenithal angle z of the Sun, and, hence, flux of radiation ( 11) can take values greater than zero, zero and negative depending on latitude ϕ of point M and angles δ and ω of the Sun's position.The positive values of insolation dW/dt correspond to the day-time, zero onesto the time of sunrise and sunset, and negative onesto the night-time.
As noted earlier, the reading of the zenithal angle z of the Sun in Figure 2 begins at noon.Therefore its greatest negative value (-90°) will correspond to sunrise, and the largest positive value (90°)to sunset.From cos z = 0 in (10) or dW/dt = 0 in (11), we derive the dependence of the hour angles of sunrise and sunset on latitude ϕ of point M and on declination δ of the Sun ω 0 = arcos(-tg ϕ ⋅ tg δ).
So, arcs SrMd and MdSd are equal to ω 0 in Figure 2, and the day-time is defined by the hour angles of the Sun, -ω 0 ≤ ω ≤ ω 0 .Therefore, the specific amount of heat, coming in per light day, i.e. to point M on the Earth's surface per day, may be defined for the day-time by integration (11).After transition from time t to the hour angle ω the daily insolation is received as: Here W is the amount of solar radiation coming to point M on 1 m 2 of the Earth's surface per day.
As time is counted on the basis of the mean movement of the Sun, a complete turn ω = 2π of the Sun around the Earth corresponds to the time of the day in minutes τ = 24⋅60 = 1440 min, i.e. on the average dt/dω = τ/2π.After the substitution of dt/dω and flux of solar radiation according to (11) in ( 13), we receive the specific daily insolation in the following form In record ( 14) we have neglected the change of distance ρ from point M to the Sun during one turn of the Earth around its axis.During this time, inclination δ of the Sun (see Figure 2) also changes due to its motion along ecliptic ЕЕ′.Let's also neglect this change.Let's accept it as simplification 5. Its influence is greatly reduced if ρ and δ refer to the middle of day.Under these conditions in integrand ( 14) only cos ω depends on the hour angle ω.
As a result of integration ( 14), the specific daily insolation is derived as: where the hour angle of the boundary of day ω 0 is defined by expression (12).
Declination δ of the Sun can be expressed from the rectangular spherical triangle γSB (see Figure 2), where angle B is right, ∠γ = ε, and side γS specifies longitude λ of the Sun, i.e. λ = γS.According to the law of sines: sin δ/sin ε = sin λ/ sin π/2 we receive the declination angle of the Sun Expression ( 15) for the daily insolation includes the relative distance ρ from the center of the Earth to the Sun.According to (2), the distance to the Sun r = ρ⋅a is defined by angle ϕ o , which is counted from perigee.
The angle of the Earth's perihelion ϕ pγ is counted from the ascending node of its orbit γ.As perigee of the Sun P S (see Figure 1b) is located diametrically opposite to perihelion of the Earth P′ E , the angle of the Sun's perigee is Because angles λ and ν (See Figure 1b) are counted from the same point, the ascending node γ, the angular distance of the Sun from its perigee will be defined by the angle In the calculation of the polar angle ϕ o , it is necessary to take into account some issues related to the cycle of changes in the angles from 0 to 2π.In our solutions of the orbital and rotational tasks the perihelion angle ϕ pγ from the mobile node γ is represented as a series of increasing numbers.Therefore, it's necessary to change angle ν of the Sun's perigee into a series of numbers with cyclic change of angles from 0 up to 2π.After the expression of longitude λ in (18), angle ϕ o must also be reduced to the range of angles from 0 to 2π.It's necessary to perform these operations during the research of the evolution of insolation for the long intervals of time.
So, longitude λ is included into distance r according to (2) and ( 18).Therefore, according to (15), the daily insolation W depends on longitude λ of the Sun and latitude ϕ of point М on the Earth's surface.Longitude λ is a short-period parameter with the cycle of change of 1 year.Eccentricity e or the trajectory parameter α 1 , the angle of inclination ε and the angle of the perihelion ϕ pγ are long-period parameters, which are part of the daily insolation W according to (15).

Change of the Sun s Longitude for Days of Year
The Sun moves non-uniformly in the orbit, for example, faster in perigee, and slower in apogee.Therefore, to calculate the insolation for the particular days of the year, it is necessary to define the Sun's longitudes appropriate to these days, i.e. it is necessary to define the dependence of the longitude on time.The time of motion t fp from perigee to the point in an orbit with angle ϕ o is defined by (4).According to (18), the Sun's longitude λ occurs in angle ϕ o .Then, using equation (4) for time t fp of the orbital motion, it is possible to define a series of longitudes, which are appropriate to days of year.This task is solved by the method of successive approximations.Initially an equal series of longitudes is set where λ0 0 = 0; Δλ0 = 2π/365.
According to (18), the appropriate series of angles ϕ o,i of the Sun's position is defined from perigee.Then according to (4), a series of the time moments t fp,1,i is defined appropriate to longitudes (19).
The derived discrete values of t fp1,i /ed, where ed = 24⋅3600 is the number of seconds in a day, differ from integers of days in a year Td i = 1, 2, … 365.Therefore, the adjusted series of longitudes is calculated where λ1 0 = 0.
Angles ϕ o,1,i are calculated from the new values of λ1 i according to (18), and then the adjusted time moments t fp,2,i are defined by equation (4).The intervals of time t fp,2,i+1 -t fp,2,i between the adjacent longitudes must be equal to length ed of a day.Therefore, the relative error value of time between the adjacent longitudes λ1 i is calculated as and also the average value of the error 22) δt fp,1,m represents the average relative error of time, which is defined by the adjusted series of longitudes.
Similarly, starting with the calculation of the adjusted series of longitudes according to (20), these operations can be repeated for derivation of new refinements λ2 i , λ3 i , λ4 i etc.The computations show that after the first refinement, the error is δt fp,1,m = 7.09⋅10 -4 , after the second refinement δt fp,2,m = 3.43⋅10 -5 and after the third refinement δt fp,3,m = 7.15⋅10 -7 .The relative error of about 10 -5 ÷ 10 -7 is quite sufficient for further insolation computation.If it is necessary to reduce it further, it is possible to perform some more refinements according to the above-stated algorithm.
Next, the computation error for longitude λ 365 for the last 365-th day of the year is computed.As a complete turn of the Sun around the Earth occurs at angle 2•π for the sidereal year P sd , the rest of the day in excess of 365 days is Δd sd = P sd -365.According to the computing longitude of the 365-th day λ 365 the excess of the day will be 364 365 where 1d is 1 day, and Δd sdc is calculated in days.
The difference between the calculated excess of the day and the actual one Δd sd represents the accumulated error of longitude λ 365 in days, i.e.Δλ 365 = Δd sdc -Δd sd .And the relative error of the longitude will be δλ 365 = Δλ 365 /Δd sd .
The following values of these quantities are received after each of the three refinements: Δλ 365 = -0.2days; -161 s; -11 s; δλ 365 = -0.8;-0.0073; -4.97⋅10 -4 .As it is seen from these data, each refinement reduces the accumulated error more than by an order of magnitude.

Change of Daily Insolation for Days of Year
Let's rewrite equation ( 15) for the daily insolation, expressing inclination δ in it through longitude λ using ( 16).Further, longitude of the Sun specified according to algorithm (20) -( 22) is meant under λ.Then according to (15), the daily insolation is )) sin (sin cos(arcsin sin cos sin sin sin ( where the hour angle of sunrises and sunsets can be written according to (12) as: The relative distance to the Sun in (24) ρ = r/a, where r is defined by (2).As according to the two-body problem (Smulsky 2007;Smulsky 1999) R p /a = (2α 1 + 1)/α 1 , then using (18) for the polar angle ϕ o , the relative distance can be written as where, as previously mentioned, the trajectory parameter α 1 is unambiguously related to the eccentricity: So, the specific daily insolation of the Earth is completely defined at any of its latitude and in any day of the year appropriate to the longitude λ by ( 24) -( 26).In the series of the specified daily longitudes λ i , where i = 1 ... 365, longitudes count from the moment of the vernal equinox, i.e.March 22, and till March 21 of the next year.Let's note that there is some difference of the calculated longitudes from the ones corresponding to calendar days.This is due to the fact that the moment of the equinox does not coincide with the beginning of a day, a calendar year does not coincide with a tropical year P tr and a leap-year jumps up by one day.Therefore, further the computing of exact dates will be considered.As it is seen from the graphs, the highest daily insolation occurs at the poles during the polar day, for example, at the North Pole W max = 46.4MJ/m 2 , and at the South Pole W max = 49.6 MJ/m 2 .At the equator, the maximum daily insolation is much lower W max = 38.8MJ/m 2 .At the same time, the minimum insolation W min = 34.1 MJ/m 2 is a little different from the maximum.However, the minimum insolation at the equator, as shown in Figure 3, is the highest among the minima at all latitudes.
It is seen from comparison W max and W min , the summer insolations are higher in the Southern Hemisphere, and the winter ones are lower than the appropriate insolations in the Northern Hemisphere.It is caused by the fact that the perihelion falls on the winter time in the Northern Hemisphere.

Insolation for a Year
As already noted, the tropical year is more than 365 days.The excess of a day in tropical year is defined by Δd tr = P tr -365.Then the complete insolation Q T for a year can be calculated by adding up daily insolations this way: where W i is the daily insolation calculated according to (24) depending on longitude λ i of the Sun, and W 365 is the daily insolation of the last day.
The main difference between our method of insolation calculation from the one of M. Milankovitch consists in calculation of insolation for the intervals of time which include a series of days, for example, for seasons, half-year and a year.M. Milankovitch derives analytical expressions for insolation for the required time intervals.To do so, he has to perform a series of mathematical transformations and approximations.Their point is to derive an analytical dependence of time on the longitude and then integrate the flux of radiation dW/dt, represented by ( 11), over the required interval of time.
In contrast with M. Milankovitch, we derive only the daily insolation by integration the flux dW/dt.Then we determine the longitude for each day on the basis of solution (4) of the two-body problem.As the longitude and days are counted from the vernal equinox, it is possible to set any period of time by number of days.In addition, the daily insolation is defined for each end of the day by expression (24).Therefore, the insolation for any interval of time is defined by summation of the daily insolation.That's the way the annual insolation has been defined in ( 27).
It was possible to do so because we have used time solutions of the two-body problem (see Section 5.2 Smulsky, 1999).M. Milankovitch relates the Sun's longitude λ with time t on the basis of the differential equation (see (34) Milankovitch, 1939): M. Milankovitch integrates equation ( 28) by the approximate analytical methods.As already noted, in (Smulsky, 1999) the exact analytical solutions of the two-body problem are given for four possible cases of motion: ellipse, parabola, hyperbola and rectilinear motion.The dependence of the time on the polar angle ϕ o or longitude λ is represented for elliptical motion by equation ( 4).Thus, the basic mathematical difference of our method from the method of M. Milankovitch is expressed by ( 4) and ( 28).
The algorithmic difference should also be noted.M. Milankovitch was based on technology of manual calculations.Therefore he had to get analytical dependences.We are based on computing.Therefore we compute insolation for intervals of time by summation or sampling from an array of daily insolations.

Insolation for Caloric Half-Year
In the astronomical theory of paleoclimate, the insolation for summer and winter half-years plays an important role.However due to the motion of the Earth's perihelion relative to the ascending node, the length of the astronomical summer and winter half-years changes.Therefore the caloric half-year of the same length were introduced.According to Milankovitch (1939), the caloric summer half-year is defined so that insolation W for any day of the half-year was more than the insolation for any day of the winter half-year.
As already noted, year P tr is defined by the number of days between two Sun's passages of the vernal equinox point γ.The value of P tr exceeds 365 by almost a quarter of a day.To find the caloric half-years, the insolation is considered for the two equal intervals of 182 days (see Figure 4).Their sum is less than P tr , but this difference is not important to define the half-years.We have developed two methods of computation for the caloric half-years.In the first method the maximum daily insolation W max for a year (see Figure 4) and the appropriate day number K max are defined.The initial number is set for the summer insolation where ΔK stp = 2 at the latitude |ϕ| ≥ 45°; ΔK stp = 5 at 25° < |ϕ| < 45°; ΔK stp = 70 at |ϕ| ≤ 25°.
Displacement ΔK stp for the initial day number was determined as a result of the summer insolation computation at different latitudes.
To simplify computing of the summer insolation, the insolation array V i2 , i2 = 1, 2… 730 is created, which includes two arrays W i from the insolations for a year: The extended array of daily insolation V i2 is required for the selected interval for the half-year (see Figure 4) not to split up.
Some variants of the summer half-years with the initial indexes in the range of j S = =K stp …K stp + 2ΔK stp , are considered.The insolation sums are computed for these half-years 182 days long: where j 3 = 0, 1…181; here js = j s and j3 = j 3 , i.e. differently written indexes are identical.
The maximum value of array VS jS with index j p0 is accepted as the main summand of insolation for the summer caloric half-year.And the complete insolation is defined by expression: In ( 32), the additional summand 0.5⋅(1 + Δd tr )⋅V jp0+182 is caused by adding the insolation for half of difference between a tropical year and for the considered period of 2⋅182 = 364 days.
So, in this method the summer caloric half-year is defined by the maximum insolation for 182 days.
In the second computing method of the summer caloric half-year, j S variants of the summer half-years are considered, and the initial and final days of the summer half-year are selected so that their insolation was more than the insolation for the nearest days of the winter half-year.These two methods are equivalent, as the selected half-years coincide.
As the insolation for the summer half-year is known, the insolation is computed for the winter caloric half-year Figure 5.The distribution of the specific heat in GJ/m 2 on the Earth's latitude in the contemporary epoch ( 1950): Q S is for summer caloric half-year; Q W is for the winter caloric half-year; Q T is for all year: in the graph Q T is reduced twice; ϕ > 0 is the northern hemisphere; ϕ < 0 is the southern hemisphere; Mil is the calculations of M. Milankovitch for the 1800 epoch Figure 5 shows the change on latitude ϕ°of summer Q S , winter Q W and twice-reduced annual insolations Q T in the contemporary epoch.The insolation for a year Q T monotonically decreases from the equator to the poles.In the meantime, at the poles it is 2.4 times less than at the equator.The winter insolation Q W near the equator has a maximum value, and at the poles it tends to zero.Summer insolation Q S has the maximum value near the tropics (ϕ = ± ε°= ± 23.4°) and the lowest values at the poles.In the meantime in the contemporary epoch, summer insolation Q S , as it is seen from the graph in Figure 5, is lower at the North Pole than at the South Pole.Also in the northern tropics the summer insolation is 1.04 times less than in the southern tropics.

Insolation in Equivalent Latitudes
The climate condition in a given geographical location, except insolation, depends on other factors.To have the ability to define changes of climate in different epochs by value of insolation, Milankovitch (1939) has proposed to express the Earth's insolation in equivalent latitudes.If in epoch T, summer insolation at latitude ϕ was the same as in contemporary epoch T 0 in summer at latitude ϕ 0, then the insolation in equivalent latitudes will be I = ϕ 0 .
Through the value of such insolation I in a point with geographical latitude ϕ, it is possible to say how much warmer or cooler the climate has become.
As shown in Figure 5, summer insolation Q S in the 1950 epoch was taken as standard.Let's denote its values at different latitudes as Q S,n,i3 and ϕ n,i3 , where n is a sign of the standard latitude, i3 is an index of the specific value of the latitude.So, there is a functional relation for standard summer insolation Q S,n,i3 (ϕ n,i3 ).To define the standard latitude, at which the insolation ).However, as the insolation Q S,n , as seen from the graph in Figure 5, has several maxima that the relation ϕ n,i3 (Q S,n,i3 ) is ambiguous.Therefore it is necessary to divide it on the monotonous sections.For example, for the northern hemisphere the section is considered from the minimum value Q S,n,min at ϕ = 87.5°andwith the imin index to the maximum value Q S,n,max at ϕ = ε°and with the imax index.
As the relation ϕ n,i3 (Q S,n,i3 ) for the standard insolation is discrete, then the insolation in equivalent latitudes is defined at any value of specific insolation Q S by parabolic interpolation where iQ is the index of origin of a parabolic section of approximation.
Coefficients A i3 , B i3 and C i3 included into (34) are calculated for any index i3 so: where Index iQ in expression ( 34) is selected so that the value of the insolation Q S considered in epoch T was between the values of Q S,n,iQ-1 and Q S,n,iQ .
So, the summer insolation I in the equivalent latitudes is defined depending on the standard summer insolation Q S,n (ϕ n ) in epoch 1950 by the known specific summer insolation Q S in kJ/m 2 in epoch T from equation ( 34) -( 37).
The insolation graph Q S is shown in Figure 5. Table 1 shows the distributions of the standard insolation in the latitude and coefficients A, B, C for the summer insolation in the northern hemisphere at the monotonous section 25°≤ϕ≤87.5°.It is similarly possible to compute the annual and winter insolations in equivalent latitudes.
The considered algorithm for computing daily W, annual Q T , summer Q S , winter Q W insolations and also of the insolation in equivalent latitudes I was realized in the MathCad environment as the program Insl2bdEn.mcd(see the Appendix).For computing of the insolation evolution, the data on evolution of the orbital and rotational motion of the Earth is additionally set as a series of the variables: T, e, ϕ pγ and ε.They are read out from the file, for example (see the Appendix), INSO_LA2004.txtor OrAl1c_8.prn.For computing summer insolation I in equivalent latitudes, a series of the variables is set in file InsCvSNJ.prn:i, Q S,n , ϕ n , A, B and C from table 1.In the program, a series of conditions, not mentioned here, is also used which provide working efficiency of the algorithm at possible combinations of used parameters.

Checking Reliability of the Algorithm
In table 14, Milankovitch (1939) represents the distribution of insolation on latitude for the caloric half-years in epoch 1800.These results are given in the canonical units used by M. Milankovitch.They can be calculated in kJ/m 2 with coefficient K Kn = 10 -5 ⋅P tr ⋅ed⋅J 0 /60 = 440.The summer and winter insolations we have computed for the northern and southern hemispheres in epoch 1800 using the above-considered algorithm.They practically coincide with the calculations by M. Milankovitch.The differences do not exceed 0.1% and, as a rule, are caused by the last digit of the numbers, which are given in table 14 (Milankovitch, 1939).
On the graphs of Figure 5, computed for epoch 1950 using the new algorithm for the summer Q S and winter insolations Q W , the calculation results of M. Milankovitch are plotted by diamonds.Although these results refer to different epochs 1950 and 1800, they are practically do not differ in the graphs.The relative difference between these results for different epochs is about hundredths of a percent and reaches the maximum value of 0.6% for latitude ϕ = 80°.
It is necessary to note that the insolation calculations of M. Milankovitch in the equatorial area are missing.It is an illustration of problems of the former method in this area.M. Milankovitch had to divide the whole of the interval of latitudes into ranges and to derive different expressions for insolation, for example, in polar and middle latitudes.
Computations were performed to verify the algorithm for insolation computation in equivalent latitudes.In the work by Sharaf and Budnikova (1969) the insolation in equivalent latitudes is given in a graphic view.However, in this work the numeric data about the evolution of the orbital and rotational motion of the Earth are missing.These data are given in the work by Laskar, Robutel, Joutel et al (2004).In this case in the program (see the Appendix), file INSO_LA2004.txtwith the parameters of the evolution of the orbital and rotational Earth's motion for 21 million years was used.These data are available at Laskar's website (2004).Based on the data in this file, we computed the insolation I in equivalent latitudes for 200 thousand years in the future.In Figure 6 our computing is compared with the calculations by Sharaf and Budnikova.As it is seen, they practically coincide.The small differences in the extreme points can be caused by two circumstances: the difference between the input data from Laskar, Robutel, Joutel, and Robutel (2004) and from Sharaf and Budnikova (1969), and also graphic character of the results of the latter.
Figure 6.The comparison of the insolations in equivalent latitudes, computed by two methods: the new method (1) and the method of Milankovitch (2): I is insolation in degrees of 15°longitude; T is time in thousands years.The insolation 1 computed by the parameters of the evolution of the orbital and rotational motion from the works by Laskar, Robutel, Joutel et al (2004) and Laskar (2004), 2the graphical results of computing by Sharaf and Budnikova (1969) So, all the checks have confirmed identity of the results of the new algorithm for insolation computation to the results of the algorithm of M. Milankovitch.It should be noted that computing insolation evolution is applicable for any models of the orbital and rotational tasks.The output of these tasks is the parameters e, ε and ϕ pγ , by which the program Insl2bdEn.mcdcomputes the insolation.Parameters of other authors can be inputted similar to inputting the Laskar parameters from file INSO_LA2004.txt(see the Appendix).In this case, it is need to comment the parameters input from file OrAl1c_8.prn.

Dynamics of Insolation in the Contemporary Epoch
As noted in the beginning, the calculations of the Earth's insolation in the contemporary epoch are used in different areas.According to the above-mentioned algorithm, using the program Ins12bd.mcd in the Appendix, insolation can be computed by setting the parameters of the Earth's orbit and the plane of the equator for epoch 2000.0 Simon, Bretagnon, Chapront et al. (1994).
e 0 = 0.01670863; ε 0 = 0.4090926; ϕ pγ0 = 1.7965956;JD p0 = 2451548.8, where JD p0 is the Julian day of passage of the middle Sun through perihelion, which is defined according to our data; JD p0 corresponds to the calendar date 05.01.2000 and the time 7:12:00 GMT.
Parameters e 0 , ε 0 and ϕ pγ0 are the averaged elements of the Earth's orbit and the positions of its axis, which motion of the "middle" Sun and the "middle" axis of the Earth's rotation are characterized by in astronomy.The Julian Day JD 0 = 2451545, the calendar date 01.01.2000 and the Greenwich time 12:00:00 correspond to epoch 2000.0.
The distribution of daily insolation W by latitude and days of a year is computed using the algorithm in item 6 of the program, described in (15).The days of year are counted from the point of a spring equinox, which day number j = 0 corresponds to.And day number j p , corresponding to the moment of passing the perihelion by the Sun, is defined by perihelion longitude λ p .The Julian date of each day JD d is defined, if it is necessary, by the number of days k d after the Julian perihelion date JD p In case a day precedes the perihelion date, k d < 0. A calendar day of the year is defined by the equivalent rendition tables of Julian days JD d to a certain date (for example, see IAA Transactions ( 2004)).The time interval of the Sun's movement from the perihelion to the equinox moment is defined by (4) at longitude λ = 0.In item 4 of the program, the amount of time from the perihelion to the day of the equinox is expressed by parameter t pγd , i.e. k d = t pγd = 76.84.Then according to (39), the Julian day of the equinox is JD γ0 = 2451625.6.The calendar date 22 March 2000 and the Greenwich time 02:24 correspond to it.
According to the series of longitudes λ i , which are considered in item 5 of the program, the days of the seasons of the year are computed, with the seasonal limits mentioned earlier: spring 0°< λ < 90°, summer 90°< λ < 180°, autumn 180°< λ < 270°and winter 270°< λ < 360°; or the days are computed for the astronomical summer: 0°< λ < 180°and the astronomical winter: 180°< λ < 360°.The insolation for these periods of time is defined by summation of the daily insolation.
According to item 8 of the program in the Appendix, the insolations for a year and the caloric half-years are computed.As already noted earlier, the distribution of insolation by latitude for 1950 and 1800 differ mainly in hundredths of a percent.Therefore, if the short-period oscillations of insolation are not important, then the insolation, computed with the changeless parameters e 0 , ε 0 , ϕ pγ0 and JD p0 , can be used for an interval of about 100 years.
For the research of reasons for short-period changes in meteorology, medicine, biology and other fields, which are caused by insolation, more exact calculations are necessary.Figure 7 shows the dynamics of insolation for 100 years.It is based on the results of the solutions for the orbital problem using the system Galactica (Smulsky, 2012) and the problem of the Earth's rotation (Smulsky, 2011) at simultaneous action of the planets, the Sun and the Moon on it.The dynamics of orbital motion is considered relative to the fixed Earth's equator, and the rotational motion of the Earth is studied relative to the fixed orbit of the Earth.Then the received parameters are recomputed into the relative variables ε and ϕ pγ , determining a position of the moving orbit relative to the moving equator.Therefore all oscillations of the Earth's orbit and the plane of the equator are included in oscillations of angles ε, ϕ pγ and eccentricity e.
The orbital parameters of the planets were defined with the interval of 1 year, therefore these computations included their short-period changes with the periods of some years.The solutions for rotational motion include daily, half-monthly, half-yearly oscillations and also the oscillations with the period of 18.6 years.If the period is longer, then the amplitude of these oscillations is higher.
Using parameters ε, ϕ pγ and e the dynamics of insolation was computed with help of the program in the Appendix.These parameters are given in table 2 with the interval of 5 years.Besides in table 2 there are given the Julian days of epoch JD and the Julian days JD p at the moments of perihelion passage by the Earth.In the program Insl2bdEn.mcdall these parameters are read out from file OrAl1c_8.prn.It is necessary to note that files OrAl1c_8.prn,InsCvSNJ.prnand program Insl2bdEn.mcdare free accessible at Smulsky (2013).
In the graph of the annual insolation Q T for latitude 80°of the northern hemisphere (see Figure 7), the oscillations with the period of 18.6 years are well visible.They are caused by the action of the precessional Moon's motion upon the rotational motion of the Earth: the Earth's axis oscillates with this period.In this case, the oscillation amplitude of the insolation is 532 J/m 2 .In the annual insolation Q T there are still half-yearly oscillations with the amplitude of 31 J/m 2 , which are not visible because of their small size on the graph.
Figure 7.The dynamics of insolation at different latitudes of the northern hemisphere for 100 years, beginning with 12.30.1949:Q S is summer insolation for caloric half-year; Q W is winter insolation for caloric half-year; Q T is total insolation for all year In insolation Q S for the summer caloric half-year at latitude 80°the oscillations of shorter period are visible together with the oscillations with the period of 18.6 years.To a greater extent they are seen in winter insolation Q W .These insolation oscillations are caused by oscillations of the perihelion angle ϕ p of the Earth's orbit with the periods of 2.75 and 11.86 years, and also its eccentricity e with the periods of 3.98 years and 11.86 years.If for the summer insolation the oscillations of perihelion are negligible, then for the winter insolation they even surpass the oscillations with the period of 18.6 years.
All three insolation: Q T , Q S and Q W have a common falling trend at 80°latitude.It is caused by the long-period decrease in the angle of inclination ε between the planes of the equator and the orbit.
As can be seen from Figure 7, the dynamics of the insolation qualitatively changes with the change of latitude ϕ = 65°, 45°, 25°and 0°.In summer insolation Q S the oscillations of the perihelion and the eccentricity are seen more clearly at 65°latitude.At 45°latitude these oscillations are more prominent for summer Q S and winter Q W insolations.At the same time, the winter and summer insolations change in antiphase to each other.Besides, in the summer and winter, insolations the falling trend, which was at the previous latitudes, is absent.
Table 2.The dynamics of the parameters of the Earth's orbit and its rotation axis for 100 years from 30.12.1949 with JD 49 = 2433280.5;T is the time in sidereal centuries; JD is the Julian day of epoch; e is the eccentricity of the orbit; ε is the angle between the moving orbit and the moving equator, ϕ pγ is the angle between the perihelion of orbit and its moving ascending node in the plane of the moving equator.JD P is the Julian day of perihelion passage by the Earth T, centuries JD e ε ϕ pγ JD P At 25°and 0°there is an ascending trend for the annual insolation Q T , and also for the summer one at 25°and the winter one at 0°.At these latitudes the winter and summer insolations also oscillate in antiphase.In their oscillations, the periods of change in the perihelion of 2.75 years and the eccentricity of 3.98 years are shown, and the 18.6-year oscillation period of the angle of inclination ε is recognized only in the change in the annual insolation Q T .
It is worth noting that in Borisenkov et al. (1983), the insolation dynamics is calculated for the summer solstice day for 200 years, beginning from 1800.The periods, amplitudes and trends of its change are similar to the ones we calculated and presented in Figure 7.In Stanhill and Moreshet (1994), the changes in annual insolation are analyzed from 1958 to 1990, which is measured by seven stations of observation.The measured insolation is approximately two times less than the calculated one due to the processes of reflection and absorption by atmosphere.However, the decreasing trend of the annual insolation in Figure 7 is observed at many stations.
The dynamics of insolation shown in Figure 7 represent a wide range of its change in time and over the surface of the Earth.These changes were received within the framework of the strictly deterministic analysis of the interaction of the Solar System bodies and the geometric characteristics of the Earth's lighting by the Sun.The deterministic approach provides different frequencies, different trends, which are differently manifested for insolations Q T , Q S and Q W .If we address any process occurring in nature, for example, to change in weather or climate, we also find a wide range of changes.As some reasons for change are unknown, and action mechanism of other reasons is not defined, the researcher is compelled to use statistical methods for a definition of probabilistic relations.Therefore the ascertained knowledge about this process is also probabilistic, which assumes unforeseen situations are possible.Nevertheless, time helps to establish some reasons, and the processes, which previously seemed random, get clear explanation.
As already mentioned, the measurements of insolation give its oscillations in a certain range, and they do not repeat from one series to another.Researchers tend to ascribe them to unstable activity of the Sun.In light of the results, presented in Figure 7, these measurements should be the way they are because the exact changes in insolation of the Earth do not repeat neither over the Earth's surface, nor in time.Apparently, it is necessary to subtract these changes from the results of observation.Then the remainder can be explained by activity of the Sun.
As already mentioned earlier, in the work by Bogdanov and Katruschenko (2008) the influence of the Moon's rotation on the insolation changes was investigated, which was calculated by the researchers taking into account only the change in the distance between the Earth and the Sun.As our computations of insolation implies, the precession of the Moon's orbit has the greatest influence on the insolation.It causes the oscillations of the annual insolation with the period of 18.6 years and the amplitude of 532 J/m 2 at the latitude of 80°, while the oscillations of the insolation with the period of 0.5 years have the amplitude of almost 20 times lower.In case of the Earth's annual orbital motion, the Earth-Sun distance changes by several orders of magnitude more than in case of the monthly orbiting of the Moon.Therefore, when performing the exact computation of insolation, the influence of the Moon's orbital motion on the insolation is practically absent.Bogdanov and Katruschenko (2008) had not known the actual mechanism of the insolation dependence on the motions of the Earth and the Moon, so they set up the hypothesis about the influence of the Moon's orbital motion on the insolation.Then they create the theory of insolation calculation depending on the Earth-Sun distance and calculate the insolation array for 300 years.As a result of the spectral analysis of the array, these researchers derived the probabilistic harmonic with the period of sidereal month and the amplitude of 84.9 mW/m 2 (milliwatts per m 2 ).This is a small value compared with the solar constant J 0 = 1395.6W/m 2 .As we noted, when using the exact method of insolation computation, the results do not seem to have the influence of the Earth's orbital motion which is several orders of magnitude higher.The Moon's orbital motion causes the precession of the Earth's rotation axis and they are in the form of the half-monthly rather than monthly oscillations.These are half-monthly oscillations that influence the insolation of the Earth, but not the oscillations of the Earth-Sun distance in the course of the Moon's orbital motion.At the same time, as we noted earlier, the Moon actually causes the greatest influence among other bodies on insolation with the period of 18.6 years and the amplitude of 532 J/m 2 in the annual insolation.However, this influence of the Moon is caused by the precession of the Moon's orbit and by the influence of this motion on the dynamics of the Earth's rotation axis.
We have considered this example in detail because it is a typical case of indeterministic approach to the study of natural processes.As noted in Smul'skii (2013), such approaches are not intended to ascertain all reasons determining the considered phenomenon.It seems sufficient to take into account some hypothetical reasons in such approaches, and then to confirm their probability of realization with a huge amount of statistical material.As shown in this example, the statistical evidence of such reasons may be available, but in reality the phenomenon is defined by quite different circumstances.
Both of these examples show that exact computations of insolation will contribute to a better quality analysis of those natural processes, which are defined by radiation of the Sun.

Conclusion
This paper considers a new mathematical theory of insolation computation.The difference from the algorithm of the Milankovitch theory consists, firstly, in a different method for computing the dependence of the Sun's longitude on time.In the new method, this relation is based on the exact solution of the two-body problem.Secondly, the new method is intended for computations.It allows the researcher to define on his own the type of processing for insolation data, and not to be limited only to those that were developed when the method was created.The description of the insolation computation method contains the assumptions used.The first simplification is related to the distinction between the Earth's orbital motion and the results of the two-body problem.The work by Smulsky (2007) shows that these distinctions are insignificant if the parameters of the orbit change according to its evolution.
The second simplification consists in the fact that the distance of point M on the Earth's surface to the Sun is equal to the distance between the centers of the Earth and the Sun at all latitudes.For insolation of the Earth, this simplification gives an insignificant error.Perhaps, for a large planet, such as Jupiter, in some cases it is necessary to consider the distinction between the distances at different latitudes.
The third simplification neglects the non-sphericity of the Earth.If insolation will be widely used in the future, it will be required to take into account the Earth's form in its computations.Perhaps, considering such fast-rotating planets, as Jupiter and Saturn, the planets' form, when accounted for, will result in significant corrections of insolation.
The fourth simplification is when we neglect changes in distance ρ and the angular position δ and λ of the Sun in the computation of integral (15) for the daily insolation.If we relate values of these parameters to the middle of the integration interval, i.e. to noon, the error will be insignificant.However, for the slow-rotating planets, such as Venus and Mercury, integral (15) must be computed with account for the changes in ρ, λ and δ.
The fifth simplification consists in the fact that the Sun's radiation is deemed to be constant.According to the researches of the last decades, small oscillations of solar radiation are revealed.If the insolation oscillations, caused by the Earth, will be deducted from the results of the measurements, then they will represent the dynamics of the Sun's insolation.It will be possible to replace J 0 in (24) with the daily average value of this quantity.Then the insolation computations will include the change of the Sun's activity.The new method of insolation computation was tested in all modes.The results, computed by it, have coincided with the results of calculations according to the Milankovitch method.On the one hand, this testifies to reliability of the new method, on the other hand, the coincidence of the results confirms the reliability of the Milankovitch method.
The approach and the program developed are applicable for insolation computation both for the Earth and for other planets.To use this approach for Jupiter and Saturn, it will probably require to take into account non-sphericity of a planet and the dependence of distance ρ on latitude φ.To determine it for Mercury and Venus, it is necessary to use the dependence of distance ρ on angles λ and δ of the Sun's position when computing the daily insolation.
When computing planet's insolation by program Ins12bdEn.mcd, it is necessary to use solar constant J 0 of the planet and also the parameters of its orbital and rotary motion.

Figure 3 .
Figure 3. Change of daily insolation W in MJ/m 2 for days of year at different latitudes ϕ from + 90°to -90°.Time Td in days is counted from the moment of the vernal equinox on March 22 (the 1950 epoch).The numerics in the graphs show the latitudes in degrees, and the maximum W max and the minimum W min daily insolation at each latitude in MJ/m 2 : NP is the North Pole; Eq is the equator; SP is the South Pole

Figure 4 .
Figure 4. To definition of the caloric half-years: 1 -the summer caloric half-year, 2 -the winter caloric half; W -the daily insolation, T d -the time in days

Table 1 .
The distribution of the standard summer insolation Q S,n in kJ/m 2 in the latitude of the northern hemisphere (epoch 1950) and coefficients A, B, C of an interpolating parabola at J 0 = 83.736kJ/(m 2 ⋅min)