# Modeling Dynamic Systems: How does the Moon move?

The bright memory of my teacher - the first dean of the Physics and Mathematics Faculty of the Novocherkassk Polytechnic Institute, the head of the department "Theoretical Mechanics" Alexander N. Cabelkova

# Introduction

August, summer is coming to an end. The people furiously jerked to the seas, but it is not surprising - the season itself. And on the habra, meanwhile, In a violent color, pseudoscience disintegrates and smells. . If we talk about the topic of this issue of "Modeling ", then we will combine business with pleasure - we will continue the promised cycle and quite a bit with this very pseudoscience for the inquisitive minds of modern youth. And the question is not really idle since school days, we used to think that our closest companion in outer space - the Moon is moving around the Earth with a period of 29.5 days, especially without going into the accompanying details. In fact, our neighbor is a unique and to some extent unique astronomical object, with the movement of which around the Earth is not so simple as, perhaps, I would like some of my colleagues from the near abroad.

So, leaving the polemic aside, we will try to consider this unconditionally beautiful, interesting and very revealing task from different sides, to the best of our competence.

often generates spores , that the moon is not a satellite of the Earth, but an independent planet of the solar system, whose orbit is disturbed by the attraction of a nearby Earth.

Let us estimate the perturbation introduced by the Sun into the Moon's trajectory relative to the Earth, as well as the perturbation introduced by the Earth into the Moon's trajectory relative to the Sun, using the criterion proposed by P. Laplace. Consider three bodies: the Sun (S), the Earth (E) and the Moon (M).

Let's assume that the Earth's orbits relative to the Earth and the Sun and the Moon are circular. Consider the motion of the Moon in a geocentric inertial frame of reference. The absolute acceleration of the Moon in the heliocentric reference frame is determined by the gravitational forces acting on it and is equal to: On the other hand, in accordance with the Coriolis theorem, the absolute acceleration of the Moon is where - a portable acceleration equal to the acceleration of the Earth relative to the Sun; - acceleration of the Moon relative to the Earth. There will be no acceleration of Coriolis here - the coordinate system chosen by us moves forward. From here we get the acceleration of the Moon relative to the Earth Part of this acceleration is is due to the attraction of the Moon to the Earth and characterizes its unperturbed geocentric movement. The rest of acceleration of the moon caused by perturbation from the sun.

If we consider the motion of the moon in the heliocentric inertial reference frame, then everything is much simpler, the acceleration is characterizes the unperturbed heliocentric motion of the moon, and the acceleration is - disturbance of this movement from the Earth.

With the parameters of the orbits of the Earth and the Moon existing at the current epoch, the inequality
is valid at each point of the trajectory of the moon. which can be verified by direct calculation, but I will refer to to the source , so as not to unnecessarily clutter up the article.

What does inequality (1) mean? Yes, that in relative terms, the effect of the Moon's perturbation of the Sun (and very significantly) is less than the effect of the attraction of the Moon to the Earth. Conversely, the disturbance of the Earth by the Earth's geocentric trajectory has a decisive influence on the nature of its motion. The influence of terrestrial gravity in this case is more significant, which means that the Moon "belongs" to the Earth by right and is its companion.

Another interesting fact is that by turning inequality (1) into an equation, one can find the locus of points where the effects of the perturbation of the Moon (and any other body) by the Earth and the Sun are the same. Unfortunately, this is not as simple as in the case of gravity. Calculations show that this surface is described by an equation of a crazy order, but close to an ellipsoid of rotation. All that we can do without unnecessary problems is to estimate the overall dimensions of this surface relative to the center of the Earth. Solving numerically the equation relative to the distance from the center of the Earth to the desired surface on a sufficient number of points, we obtain a section of the sought-for surface by the ecliptic plane For clarity, the geocentric orbit of the moon is also shown here, and the sphere of gravity of the Earth, which we found above relative to the Sun, found above. It can be seen from the figure that the sphere of influence, or the sphere of gravitational action of the Earth relative to the Sun, is a surface of revolution with respect to the X axis, oblate along a straight line connecting the Earth and the Sun (along the axis of eclipses). The Moon's orbit lies deep within this imaginary surface.

For practical calculations, this surface can be conveniently approximated by a sphere with the center at the center of the Earth and a radius equal to where m is the mass of the smaller celestial body; M is the mass of a larger body, in the gravitational field of which a smaller body moves; a is the distance between the centers of the bodies. In our case, This unfinished million kilometers is the theoretical limit for which the power of the old woman of the Earth does not spread - her influence on the trajectory of astronomical objects is so small that they can be neglected. So, to launch the Moon in a circular orbit at a distance of 38.4 million kilometers from the Earth (as some linguists do) will not work, it is physically impossible.

This sphere, for comparison, is shown in the figure with a blue dotted line. At estimated calculations it is considered that the body inside this sphere will experience gravitation exclusively from the side of the Earth. If the body is located outside this sphere - we believe that the body moves in the gravitational field of the Sun. In practical cosmonautics, the method of conjugation of conic sections is known, which makes it possible to approximate the trajectory of a spacecraft using the solution of the two-body problem. In this case, all the space that overcomes the apparatus is divided into similar spheres of influence.

For example, it is now clear that in order to have a theoretical opportunity to make maneuvers to enter the near-moon orbit, the spacecraft must get inside the sphere of action of the Moon relative to the Earth. Its radius can easily be calculated by formula (3) and it is equal to 66 thousand kilometers.

Thus, the Moon can rightly be considered a satellite of the Earth. However, due to the significant influence of the gravitational field of the Sun, it does not move in the central gravitational field, and hence its trajectory is not a conical section.

# 3. The three-body problem in the classical formulation of

So, let us consider a model problem in the general formulation known in celestial mechanics as a three-body problem. Consider three bodies of arbitrary mass located arbitrarily in space and moving solely under the action of the forces of mutual gravitational attraction We consider bodies material points. The position of the bodies will be counted in an arbitrary basis, with which the inertial frame of reference is connected. Oxyz . The position of each of the bodies is given by the radius vector, respectively, , and . Each body is affected by gravitational attraction from the other two bodies, and in accordance with the third axiom of the point dynamics (Newton's third law) Let us write the differential equations of motion of each point in the vector form or, taking into account (4) In accordance with the law of universal gravitation, the forces of interaction are directed along the vectors Along each of these vectors we will release the corresponding unit vector then each of the gravitational forces is calculated by the formula With all this taken into account, the system of equations of motion takes the form We introduce the notation adopted in celestial mechanics Is the gravitational parameter of the attracting center. Then the equations of motion will take the final vector form # 4. Rationing of equations to dimensionless variables

A fairly popular technique for mathematical modeling is to reduce the differential equations and other relations that describe the process to dimensionless phase coordinates and dimensionless time. The other parameters are also normalized. This allows us to consider, albeit with the use of numerical modeling, but in a fairly general form a whole class of typical problems. The question of how much this is justified in each problem being solved is left open, but I agree that in this case this approach is quite just.

So, let us introduce some abstract celestial body with the gravitational parameter , such that the period of revolution of the satellite in an elliptical orbit with the semimajor axis around it is . All these quantities, by virtue of the laws of mechanics, are related by the relation We introduce a change of parameters. For the position of the points of our system where - dimensionless radius vector of the i-th point;

for the gravitational parameters of bodies where - dimensionless gravitational parameter of the i-th point;

for the time where - dimensionless time.

Now we recount the acceleration of system points through these dimensionless parameters. We apply a direct two-fold differentiation with respect to time. For the speeds For the acceleration When you substitute the resulting relationships in the equations of motion, everything elegantly collapses into beautiful equations: This system of equations is still considered not to be integrable in analytic functions. Why is it considered and not is? Because the success of the theory of a function of a complex variable led to the fact that the general solution of the three-body problem appeared in 1912 - Karl Sundman found an algorithm for finding coefficients for infinite series with respect to a complex parameter, theoretically being a general solution of the three-body problem. But for the application of the Sundman series in practical calculations with the accuracy required for them, it is necessary to obtain the number of terms of these series that this task much exceeds the capabilities of computers even to the present day.

Therefore, numerical integration is the only way to analyze the solution of equation (5)

# 5. Calculation of the initial conditions: we extract the initial data

As I wrote earlier , before starting numerical integration, one should be concerned with calculating the initial conditions for the problem being solved. In the problem under consideration, the search for the initial conditions turns into an independent subtask, since the system (5) gives us nine second-order scalar equations, which increases the order of the system by a factor of 2 in the transition to the normal Cauchy form. That is, we need to calculate as many as 18 parameters - the initial positions and the components of the initial velocity of all points of the system. Where do we get data on the location of the heavenly bodies of interest to us? We live in a world where a person walked on the Moon - naturally, humanity must have information about how this Moon moves and where it is.

That is, you say, you, dude, suggest that we take from the shelves thick astronomical directories, blow off dust from them They did not guess! I propose to go for these data to those who actually walked on the Moon, to NASA, namely to the Jet Propulsion Laboratory, Pasadena, California. Here it is JPL Horizonts web interface .

Here, after spending a little time studying the interface, we will get all the data we need. Let's choose a date, for example, yes, we do not care, but let it be July 2? 2018 UT 20:21. Just at this moment, the full phase of the lunar eclipse was observed. The program will give us a huge portess

The full conclusion for the ephemeris of the Moon on ??? 20:21 (origin in the center of the Earth) [/b]
` ************************************************** ***************************** Revised: Jul 3? 2013 Moon /(Earth) 301  GEOPHYSICAL DATA (updated 2018-Aug-13): Vol. Mean Radius, km = ??? + -??? Mass, x10 ^ 22 kg = ???r3r31865.  Radius (gravity), km = 1738.0 Surface emissivity = ???r3r31865.  Radius (IAU), km = 1737.4 GM, km ^ 3 /s ^ 2 = ???r3r31865.  Density, g /cm ^ 3 = ??? GM 1-sigma, km ^ 3 /s ^ 2 = + -???r3r31865.  V (?0) = +??? Surface accel., M /s ^ 2 = ???r3r31865.  Earth /Moon mass ratio = ??? Farside crust. thick. = ~ 80 - 90 km Mean crustal density = ??? + -. 07 g /cm ^ 3 Nearside crust. thick. = 58 + -8 km Heat flow, Apollo 15 = 3.1 + -. 6 mW /m ^ 2 k2 = ???r3r31865.  Heat flow, Apollo 17 = 2.2 + -. 5 mW /m ^ 2 Rot. Rate, rad /s = ???r3r31865.  Geometric Albedo = ???r3r31865.   Mean angular diameter = 31'05.2 "Orbit period = ??? d Obliquity to orbit = ??? deg Eccentricity = ???r3r31865.  Semi-major axis, a = 384400 km Inclination = ??? deg Mean motion, rad /s = ???x10 ^ -6 Nodal period = ??? d Apsidal period = ??? d Mom. of inertia C /MR ^ 2 = ???r3r31865.  beta (C-A /B), x10 ^ -4 = ??? gamma (B-A /C), x10 ^ -4 = ???r3r31865.   Perihelion Aphelion Mean Solar Constant (W /m ^ 2) 1414 + -??? + -??? + -7 Maximum Planetary IR (W /m ^ 2) ??? Minimum Planetary IR (W /m ^ 2) ???.??? ************************************************** *****************************   ************************************************** ***************************** Ephemeris /WWW_USER Wed Aug ???:45:??? Pasadena, USA /Horizons ************************************************** ***************************** Target body name: Moon (301) {source: DE431mx} Center body name: Earth (399) {source: DE431mx} Center-site name: BODY CENTER ************************************************** ***************************** Start time: A.D. 2018-Jul-???: 21: ??? TDB Stop time: A.D. 2018-Jul-???: 21: ??? TDB Step-size: 0 steps ************************************************** ***************************** Center geodetic: ???.??? {E-lon (deg), Lat (deg), Alt (km)} Center cylindric: ???.??? {E-lon (deg), Dxy (km), Dz (km)} Center radii: 6378.1 x 6378.1 x 6356.8 km {Equator, meridian, pole} Output units: AU-D Output type: GEOMETRIC cartesian states Output format: 3 (position, velocity, LT, range, range-rate) Reference frame: ICRF /J???r3r31865.  Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ************************************************** ***************************** JDTDB X Y Z VX VY VZ LT RG RR ************************************************** ***************************** \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E-03 Y = -???E-03 Z = ???E-06 VX = ???E-04 VY = ???E-04 VZ = -???E-05 LT = ???E-05 RG = ???E-03 RR = -???E-06 \$\$ EOE ************************************************** ***************************** Coordinate system description:  Ecliptic and Mean Equinox of Reference Epoch  Reference epoch: J???r3r31865.  XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of ??? arcseconds wrt ICRF equator (IAU76) X-axis: out along the ascending node of the instantaneous plane of the Earth's orbit and the Earth's mean at the epoch reference Z-axis: perpendicular to the xy plane in the directional (+ or -) sense of Earth's north pole at the reference epoch.  Symbol meaning[1 au= 149597870.700 km, 1 day= 86400.0 s]:  JDTDB Julian Day Number, Barycentric Dynamical Time X X-component of position vector (au) Y Y-component of position vector (au) Z Z-component of position vector (au) VX X-component of velocity vector (au /day) VY Y-component of velocity vector (au /day) VZ Z-component of velocity vector (au /day) LT One-way down-leg Newtonian light-time (day) RG Range; distance from coordinate center (au) RR Range-rate; radial velocity wrt coord. center (au /day)  Geometric states /elements have no aberrations applied.  Computations by Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect: telnet: //ssd.jpl.nasa.gov: 6775 (via browser) http://ssd.jpl.nasa.gov/?horizons telnet ssd.jpl.nasa.gov 6775 (via command-line) Author: [email protected] ************************************************** ***************************** `

Br-rr, what is it? Without panic, for someone who taught astronomy well in school, mechanics and mathematics are scared of nothing. So, the most important final sought-for coordinates and components of the Moon's speed are

` \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E-03 Y = -???E-03 Z = ???E-06 VX = ???E-04 VY = ???E-04 VZ = -???E-05 LT = ???E-05 RG = ???E-03 RR = -???E-06 \$\$ EOE `

Yes, yes, they are Cartesian! If you carefully read the whole of the cloth, we learn that the origin of this coordinate system coincides with the center of the Earth. The XY plane lies in the plane of the Earth's orbit (ecliptic plane) for the J2000 era. The X axis is directed along the line of intersection of the plane of the equator of the Earth and the ecliptic to the point of the vernal equinox. The Z axis looks in the direction of the north pole of the Earth perpendicular to the ecliptic plane. Well, the Y axis complements all this happiness to the right three vectors. By default, the unit of measurement of coordinates: astronomical units (umnichki from NASA lead and the amount of the autonomous unit in kilometers). Units of measurement of speed: astronomical units per day, the day is assumed equal to 86400 seconds. Full stuffing!

Similar information we can get for the Earth

The complete derivation of the Earth's ephemerides as of ??? 20:21 (origin in the center of mass of the solar system) [/b]
` ************************************************** ***************************** Revised: Jul 3? 2013 Earth 399  GEOPHYSICAL PROPERTIES (revised Aug 1? 2018): Vol. Mean Radius (km) = ??? + -??? Mass x10 ^ 24 (kg) = ??? + -???r3r31865.  Equ. radius, km = ??? Mass layers: Polar axis, km = ??? Atmos = 5.1 x 10 ^ 18 kg Flattening = 1 /??? oceans = 1.4 x 10 ^ 21 kg Density, g /cm ^ 3 = ??? crust = 2.6 x 10 ^ 22 kg J2 (IERS 2010) = ??? mantle = ??? x 10 ^ 24 kg g_p, m /s ^ 2 (polar) = ??? outer core = ??? x 10 ^ 24 kg g_e, m /s ^ 2 (equatorial) = ??? inner core = ??? x 10 ^ 22 kg g_o, m /s ^ 2 = ??? Fluid core rad = 3480 km GM, km ^ 3 /s ^ 2 = ??? Inner core rad = 1215 km GM 1-sigma, km ^ 3 /s ^ 2 = ??? Escape velocity = ??? km /s Rot. Rate (rad /s) = ??? Surface Area: Mean sidereal day, hr = ??? land = ??? x 10 ^ 8 km Mean solar day 2000.? s = ??? sea = ??? x 10 ^ 8 km Mean solar day 1820.? s = ???r3r31865.  Moment of inertia = ??? Love no., K2 = ???r3r31865.  Mean Temperature, K = 270 Atm. pressure = 1.0 bar Vis. mag. V (?0) = -??? Volume, km ^ 3 = ??? x 10 ^ 12 Geometric Albedo = ??? Magnetic moment = ??? gauss Rp ^ 3 Solar Constant (W /m ^ 2) = 1367.6 (mean), 1414 (perihelion), 1322 (aphelion)  ORBIT CHARACTERISTICS: Obliquity to orbit, deg = ??? Sidereal orb period = ??? y Orbital speed, km /s = ??? Sidereal orb period = ??? d Mean daily motion, deg /d = ??? Hill's sphere radius = ???r3r31865.  ************************************************** *****************************   ************************************************** ***************************** Ephemeris /WWW_USER Wed Aug ???:16:??? Pasadena, USA /Horizons ************************************************** ***************************** Target body name: Earth (399) {source: DE431mx} Center body name: Solar System Barycenter (0) {source: DE431mx} Center-site name: BODY CENTER ************************************************** ***************************** Start time: A.D. 2018-Jul-???: 21: ??? TDB Stop time: A.D. 2018-Jul-???: 21: ??? TDB Step-size: 0 steps ************************************************** ***************************** Center geodetic: ???.??? {E-lon (deg), Lat (deg), Alt (km)} Center cylindric: ???.??? {E-lon (deg), Dxy (km), Dz (km)} Center radii: (undefined) Output units: AU-D Output type: GEOMETRIC cartesian states Output format: 3 (position, velocity, LT, range, range-rate) Reference frame: ICRF /J???r3r31865.  Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ************************************************** ***************************** JDTDB X Y Z VX VY VZ LT RG RR ************************************************** ***************************** \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E-01 Y = -???E-01 Z = -???E-05 VX = ???E-02 VY = ???E-03 VZ = ???E-07 LT = ???E-03 RG = ???E + 00 RR = -???E-05 \$\$ EOE ************************************************** ***************************** Coordinate system description:  Ecliptic and Mean Equinox of Reference Epoch  Reference epoch: J???r3r31865.  XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of ??? arcseconds wrt ICRF equator (IAU76) X-axis: out along the ascending node of the instantaneous plane of the Earth's orbit and the Earth's mean at the epoch reference Z-axis: perpendicular to the xy plane in the directional (+ or -) sense of Earth's north pole at the reference epoch.  Symbol meaning[1 au= 149597870.700 km, 1 day= 86400.0 s]:  JDTDB Julian Day Number, Barycentric Dynamical Time X X-component of position vector (au) Y Y-component of position vector (au) Z Z-component of position vector (au) VX X-component of velocity vector (au /day) VY Y-component of velocity vector (au /day) VZ Z-component of velocity vector (au /day) LT One-way down-leg Newtonian light-time (day) RG Range; distance from coordinate center (au) RR Range-rate; radial velocity wrt coord. center (au /day)  Geometric states /elements have no aberrations applied.  Computations by Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect: telnet: //ssd.jpl.nasa.gov: 6775 (via browser) http://ssd.jpl.nasa.gov/?horizons telnet ssd.jpl.nasa.gov 6775 (via command-line) Author: [email protected] ************************************************** ***************************** `

Here the barycenter (center of mass) of the solar system is chosen as the origin of coordinates. The data of interest to us is

` \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E-01 Y = -???E-01 Z = -???E-05 VX = ???E-02 VY = ???E-03 VZ = ???E-07 LT = ???E-03 RG = ???E + 00 RR = -???E-05 \$\$ EOE `

For the Moon, we need coordinates and speed relative to the barycenter of the solar system, we can count them, and we can ask NASA to give us such data

Full withdrawal of the ephemeris of the Moon on ??? 20:21 (the origin of coordinates in the center of mass of the solar system) [/b]
` ************************************************** ***************************** Revised: Jul 3? 2013 Moon /(Earth) 301  GEOPHYSICAL DATA (updated 2018-Aug-13): Vol. Mean Radius, km = ??? + -??? Mass, x10 ^ 22 kg = ???r3r31865.  Radius (gravity), km = 1738.0 Surface emissivity = ???r3r31865.  Radius (IAU), km = 1737.4 GM, km ^ 3 /s ^ 2 = ???r3r31865.  Density, g /cm ^ 3 = ??? GM 1-sigma, km ^ 3 /s ^ 2 = + -???r3r31865.  V (?0) = +??? Surface accel., M /s ^ 2 = ???r3r31865.  Earth /Moon mass ratio = ??? Farside crust. thick. = ~ 80 - 90 km Mean crustal density = ??? + -. 07 g /cm ^ 3 Nearside crust. thick. = 58 + -8 km Heat flow, Apollo 15 = 3.1 + -. 6 mW /m ^ 2 k2 = ???r3r31865.  Heat flow, Apollo 17 = 2.2 + -. 5 mW /m ^ 2 Rot. Rate, rad /s = ???r3r31865.  Geometric Albedo = ???r3r31865.   Mean angular diameter = 31'05.2 "Orbit period = ??? d Obliquity to orbit = ??? deg Eccentricity = ???r3r31865.  Semi-major axis, a = 384400 km Inclination = ??? deg Mean motion, rad /s = ???x10 ^ -6 Nodal period = ??? d Apsidal period = ??? d Mom. of inertia C /MR ^ 2 = ???r3r31865.  beta (C-A /B), x10 ^ -4 = ??? gamma (B-A /C), x10 ^ -4 = ???r3r31865.   Perihelion Aphelion Mean Solar Constant (W /m ^ 2) 1414 + -??? + -??? + -7 Maximum Planetary IR (W /m ^ 2) ??? Minimum Planetary IR (W /m ^ 2) ???.??? ************************************************** *****************************   ************************************************** ***************************** Ephemeris /WWW_USER Wed Aug ???:19:??? Pasadena, USA /Horizons ************************************************** ***************************** Target body name: Moon (301) {source: DE431mx} Center body name: Solar System Barycenter (0) {source: DE431mx} Center-site name: BODY CENTER ************************************************** ***************************** Start time: A.D. 2018-Jul-???: 21: ??? TDB Stop time: A.D. 2018-Jul-???: 21: ??? TDB Step-size: 0 steps ************************************************** ***************************** Center geodetic: ???.??? {E-lon (deg), Lat (deg), Alt (km)} Center cylindric: ???.??? {E-lon (deg), Dxy (km), Dz (km)} Center radii: (undefined) Output units: AU-D Output type: GEOMETRIC cartesian states Output format: 3 (position, velocity, LT, range, range-rate) Reference frame: ICRF /J???r3r31865.  Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ************************************************** ***************************** JDTDB X Y Z VX VY VZ LT RG RR ************************************************** ***************************** \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E-01 Y = -???E-01 Z = -???E-05 VX = ???E-02 VY = ???E-03 VZ = -???E-05 LT = ???E-03 RG = ???E + 00 RR = -???E-05 \$\$ EOE ************************************************** ***************************** Coordinate system description:  Ecliptic and Mean Equinox of Reference Epoch  Reference epoch: J???r3r31865.  XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of ??? arcseconds wrt ICRF equator (IAU76) X-axis: out along the ascending node of the instantaneous plane of the Earth's orbit and the Earth's mean at the epoch reference Z-axis: perpendicular to the xy plane in the directional (+ or -) sense of Earth's north pole at the reference epoch.  Symbol meaning[1 au= 149597870.700 km, 1 day= 86400.0 s]:  JDTDB Julian Day Number, Barycentric Dynamical Time X X-component of position vector (au) Y Y-component of position vector (au) Z Z-component of position vector (au) VX X-component of velocity vector (au /day) VY Y-component of velocity vector (au /day) VZ Z-component of velocity vector (au /day) LT One-way down-leg Newtonian light-time (day) RG Range; distance from coordinate center (au) RR Range-rate; radial velocity wrt coord. center (au /day)  Geometric states /elements have no aberrations applied.  Computations by Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect: telnet: //ssd.jpl.nasa.gov: 6775 (via browser) http://ssd.jpl.nasa.gov/?horizons telnet ssd.jpl.nasa.gov 6775 (via command-line) Author: [email protected] ************************************************** ***************************** `

` \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E-01 Y = -???E-01 Z = -???E-05 VX = ???E-02 VY = ???E-03 VZ = -???E-05 LT = ???E-03 RG = ???E + 00 RR = -???E-05 \$\$ EOE `

Wonderful! Now it is necessary to process the received data lightly with a file.

# 6. 38 parrots and one parrot's wings

To begin with, we determine the scale, because our equations of motion (5) are written in dimensionless form. Data provided by NASA themselves tell us that for the scale of the coordinates it is worth taking one astronomical unit. Accordingly, as a reference body, to which we will normalize the masses of other bodies, we take the Sun, and as the time scale - the period of revolution of the Earth around the Sun.

All this is of course very good, but we did not set the initial conditions for the Sun. "Why?" - I would be asked by some linguist. And I would say that the Sun is not at all motionless, but also rotates in its orbit around the center of mass of the solar system. This can be seen by looking at the NASA data for the Sun

` \$\$ SOE ??? = A.D. 2018-Jul-???: 21: ??? TDB X = ???E + 04 Y = ???E + 06 Z = -???E + 04 VX = -???E-02 VY = ???E-03 VZ = ???E-04 LT = ???E + 00 RG = ???E + 06RR = ???E-03 \$\$ EOE `

Looking at the parameter RG, we see that the Sun rotates around the barycenter of the solar system, and on ??? the center of the star is located from it at a distance of a million kilometers. The radius of the Sun, for reference - 696 thousand kilometers. That is, the barycenter of the solar system lies half a million kilometers from the surface of the luminary. Why? Because all other bodies interacting with the Sun also report acceleration to it, mainly, of course, the heavier Jupiter. Accordingly, the Sun also has its own orbit.

Of course, we can choose these data as initial conditions, but no - we solve the model problem of three bodies, and Jupiter and other characters do not enter it. So to the detriment of realism, knowing the position and speeds of the Earth and the Moon, we recount the initial conditions for the Sun, so that the center of mass of the Sun-Earth-Moon system is at the origin. For the center of mass of our mechanical system, the equation
is valid. We place the center of mass at the origin, that is, we set , then whence We proceed to dimensionless coordinates and parameters by choosing  Differentiating (6) with respect to time and passing to the dimensionless time, we obtain the relation for the velocities where Now write a program that will form the initial conditions in the "parrots" we have chosen. On what will we write? Of course, on Python! After all, as you know, this is the best language for mathematical modeling.

However, if we get away from sarcasm, we really will try for this purpose a python, and why not? I'm sure to provide a link to all the code in My profile is Github .

Calculation of initial conditions for the Moon-Earth-Sun system [/b]
`  `  ## The initial data of the task## Gravitational constantG = ???e-11# Masses of bodies (the Moon, the Earth, the Sun)m =[7.349e22, 5.792e24, 1.989e30]# Calculate the gravity parameters of the bodiesmu =[]print ("Gravitational parameters of bodies")for i, mass in enumerate (m):mu.append (G * mass)print ("mu[" + str(i) + "]=" + str (mu[i]))# We normalize the gravitational parameters to the Sunkappa =[]print ("Normalized gravitational parameters")for i, gp in enumerate (mu):kappa.append (gp /mu)print ("xi[" + str(i) + "]=" + str (kappa[i]))print ("n")# The astronomical unit isa = ???e11import math# Scale of the dimensionless time, cT = 2 * math.pi * a * math.sqrt (a /mu)print ("Time scale T =" + str (T) + "n")# NASA coordinates for LuxL = ???E-01yL = -???E-01zL = -???E-05import numpy as npxi_10 = np.array ([xL, yL, zL])print ("The initial position of the moon, as:" + str (xi_10))# NASA coordinates for the EarthxE = ???E-01yE = -???E-01zE = -???E-05xi_20 = np.array ([xE, yE, zE])print ("The initial position of the Earth, as .:" + str (xi_20))# Calculate the initial position of the Sun, assuming that the origin is at the center of mass of the entire systemxi_30 = - kappa* xi_10 - kappa* xi_20print ("The initial position of the Sun, as .:" + str (xi_30))# We introduce constants for the calculation of the dimensionless velocitiesTd = ???r3r31872. u = math.sqrt (mu/a) /2 /math.piprint ("n")The initial velocity of the Moon isvxL = ???E-02vyL = ???E-03vzL = -???E-05vL0 = np.array ([vxL, vyL, vzL])uL0 = np.array ([0.0, 0.0, 0.0])for i, v in enumerate (vL0):vL0[i]= v * a /TduL0[i]= vL0[i]/uprint ("Initial velocity of the Moon, m /s:" + str (vL0))print ("- //- dimensionless:" + str (uL0))The initial velocity of the Earth isvxE = ???E-02vyE = ???E-03vzE = ???E-07vE0 = np.array ([vxE, vyE, vzE])uE0 = np.array ([0.0, 0.0, 0.0])for i, v in enumerate (vE0):vE0[i]= v * a /TduE0[i]= vE0[i]/uprint ("Earth's initial velocity, m /s:" + str (vE0))print ("- //- dimensionless:" + str (uE0))The initial velocity of the Sun isvS0 = - kappa* vL0 - kappa* vE0uS0 = - kappa* uL0 - kappa* uE0print ("Initial velocity of the Sun, m /s:" + str (vS0))print ("- //- dimensionless:" + str (uS0))`  `

The exhaust of the program is

` Gravitational parameters of bodies mu= ???r3r31865.  mu= ???r3r31865.  mu= ???e + 20 Normalized gravity parameters xi= ???e-08 xi= ???e-06 xi= ???r3r31865.   The time scale is T = ???r3r31865.   The initial position of the moon, as:[5.77103476e-01 -8.32119380e-01 -4.85579076e-05] The initial position of the Earth, as:[5.75566367e-01 -8.29881892e-01 -5.36699450e-05] The initial position of the Sun, as:[-1.69738146e-06 2.44737475e-06 1.58081871e-10]  Initial velocity of the Moon, m /s:[24838.98933473 17310.56333294 -89.15979106] - //- dimensionless:[5.24078311 3.65235907 -0.01881184] The initial velocity of the Earth, m /s:[2.40435899e+04 1.67586567e+04 5.93870516e-01] - //- dimensionless:[5.07296163e+00 3.53591219e+00 1.25300854e-04] Initial speed of the Sun, m /s:[-7.09330769e-02 -4.94410725e-02 1.56493465e-06] - //- dimensionless:[-1.49661835e-05 -1.04315813e-05 3.30185861e-10] `

# 7. Integration of the equations of motion and analysis of the results

Actually integration itself is reduced to a more or less standard procedure for SciPy preparation of a system of equations: the transformation of the ODE system to the Cauchy form and the invocation of the corresponding solver functions. To transform the system to the Cauchy form, we recall that Then, introducing the state vector of the system we reduce (7) and (5) to one vector equation To integrate (8) with the existing initial conditions, we write a little, very little code

Integration of the equations of motion in the three-body problem [/b]
`  `  #Computation of generalized acceleration vectors#def calcAccels (xi):k = 4 * math.pi ** 2xi12 = xi- xixi13 = xi- xixi23 = xi- xis12 = math.sqrt (np.dot (xi1? xi12))s13 = math.sqrt (np.dot (xi1? xi13))s23 = math.sqrt (np.dot (xi2? xi23))a1 = (k * kappa/s12 ** 3) * xi12 + (k * kappa/s13 ** 3) * xi13a2 = - (k * kappa/s12 ** 3) * xi12 + (k * kappa/s23 ** 3) * xi23a3 = - (k * kappa/s13 ** 3) * xi13 - (k * kappa/s23 ** 3) * xi23return[a1, a2, a3]#The system of equations in the normal form of Cauchy#def f (t, y):n = 9dydt = np.zeros ((2 * n))for i in range (? n):dydt[i]= y[i + n]xi1 = np.array (y[0:3])xi2 = np.array (y[3:6])xi3 = np.array (y[6:9])accels = calcAccels ([xi1, xi2, xi3])i = nfor accel in accels:for a in accel:dydt[i]= ai = i + 1return dydtInitial conditions of the Cauchy problemy0 =[xi_10, xi_10, xi_10,xi_20, xi_20, xi_20,xi_30, xi_30, xi_30,uL0, uL0, uL0,uE0, uE0, uE0,uS0, uS0, uS0]#We integrate the equations of motion## The starting time ist_begin = 0# The final time ist_end = 30.7 * Td /T;The number of trajectory points we are interested in isN_plots = 1000# The time step between pointsstep = (t_end - t_begin) /N_plotsimport scipy.integrate as spisolver = spi.ode (f)solver.set_integrator ('vode', nsteps = 5000? method = 'bdf', max_step = 1e-? rtol = 1e-12)solver.set_initial_value (y? t_begin)ts =[]ys =[]i = 0while solver.successful () and solver.t  <= t_end:solver.integrate (solver.t + step)ts.append (solver.t)ys.append (solver.y)print (ts[i], ys[i])i = i + 1`   Let's see what we got. The spatial trajectory of the Moon was obtained for the first 29 days from the initial pointchosen by us.   ` as well as its projection into the plane of the ecliptic. "Hey, uncle, what are you vparivesh? This is a circle! ".

Firstly, it's not a circle - the displacement of the projection of the trajectory from the origin to the right and down is noticeable. Secondly - do not notice anything? No, well, the truth? I promise to prepare a rationale for this (based on the analysis of account errors and NASA data) that the resulting trajectory offset is not a consequence of integration errors. While I suggest to the reader to take my word for it - this displacement is a consequence of the solar perturbation of the lunar trajectory. We are spinning another turn  In how! And pay attention to the fact that, based on the initial data of the problem, the Sun is just in the direction to which the trajectory of the Moon is shifted on each turn. Yes, this arrogant Sun steals from us our beloved companion! Oh, this is the sun!

It can be concluded that solar gravity affects the Moon's orbit quite significantly - the old lady does not walk across the sky twice in the same way. The picture for six months of the movement allows (at least qualitatively) to see this (picture is clickable) Interesting? Still would. Astronomy in general is an interesting science.

# Postscript

In the university, where I studied and worked for almost seven years - the Novocherkassk Polytechnic - the zone Olympiad of students on theoretical mechanics of higher educational institutions of the North Caucasus was held annually. Three times we also took part in the All-Russian Olympiad. At the opening, our main "Olympian", Professor Kondratenko AI, always said: "Academician Krylov called mechanics the poetry of the exact sciences."

I love mechanics. All the good things that I achieved in my life and career came about thanks to this science and my wonderful teachers. I respect mechanics.

Therefore, I will never allow to scoff at this science and brazenly exploit it for my own purposes to anyone, be it at least three times Doctor of Science and four times a linguist, and have developed at least a million curricula. I sincerely believe that writing articles on a popular public resource should provide for their thorough proofreading, normal design (LaTeX formulas are not the whim of resource developers!) And the absence of errors leading to results violating the laws of nature. The latter in general "mast hev."

I often tell my students: "The computer frees up your hands, but that does not mean that you need to turn off the brain."

To appreciate and respect mechanics, I urge you, my dear readers. I will gladly answer any questions, and the original text of the example of solving the three-body problem in Python, as promised, is I spread in my profile Github .

Thank you for attention!
+ 0 -

• • • 