📜 ⬆️ ⬇️

Modeling Dynamic Systems: How does the moon move?

The blessed memory of my teacher - the first dean of the Faculty of Mathematics and Physics of the Novocherkassk Polytechnic Institute, Head of the Department "Theoretical Mechanics" Alexander N. Kabelkov

Introduction


August, summer is coming to an end. People violently rushed to the sea, but it is not surprising - the season. And on Habré, meanwhile, with pompous color it dissolves and smells of pseudoscience . If we talk about the topic of this issue of "Modeling ...", then we will combine the pleasant with the useful in it - we will continue the promised cycle and just a bit more will fight with this very pseudoscience for the inquisitive minds of modern youth.


And the question is really not idle - from school years we used to think that our closest satellite in outer space - the Moon moves around the Earth with a period of 29.5 days, especially without going into the accompanying details. In fact, our neighbor is a peculiar and to some extent unique astronomical object, with the movement of which around the Earth is not so simple as it might have been desirable for some of my colleagues from the near abroad.
')
So, leaving the controversy aside, we will try from different sides, to the best of our competence, to consider this unconditionally beautiful, interesting and very revealing task.

1. The law of world wideness and what conclusions we can make of it


Discovered in the second half of the 17th century, by Sir Isaac Newton, the law of universal perception speaks of the fact that the moon is attracted to the Earth (and the Earth to the Moon!) With a force directed along the straight line connecting the centers of the considered celestial bodies and equal to

F 1 , 2 = G f r a c m 1m 2 r 2 1 , 2


where m 1 , m 2 are the masses, respectively, of the Moon and the Earth; G = 6.67e-11 m 3 / (kg * s 2 ) - gravitational constant; r 1,2 is the distance between the centers of the moon and the earth. If we take into account only this force, then, having solved the problem of the motion of the Moon as a satellite of the Earth and having learned to calculate the position of the Moon in the sky against the background of stars, we will soon be convinced by direct measurements of the equatorial coordinates of the Moon that in our conservatory not everything is as smooth as I would like to. And the point here is not in the law of universal perception (and in the early stages of the development of celestial mechanics such thoughts were expressed quite often), but in unaccounted perturbation of the motion of the moon by other bodies. What kind? We look at the sky and our gaze immediately rests on a hefty, weighing as much as 1.99e30 kilograms plasma ball right under our nose - the Sun. Is the moon attracted to the sun? More like, with a force equal in magnitude

F 1 , 3 = G f r a c m 1m 3 r 2 1 , 3


where m 3 is the mass of the Sun; r 1,3 is the distance from the moon to the sun. Compare this force with the previous one.

 f r a c F 1 , 3 F 1 , 2 = f r a c G  f r a c m 1m 3 r 2 1 , 3 G f r a c m 1m2r21,2= fracm3m2 left( fracr1,2r1,3 right)2


Take the position of the bodies, in which the attraction of the Moon to the Sun will be minimal: all three bodies are on one straight line and the Earth is located between the Moon and the Sun. In this case, our formula will look like:

 fracF1,3F1,2= fracm3m2 left( frac rhoa+ rho right)2


Where  rho=3,844 cdot108 , m - the average distance from the Earth to the Moon; a=1,496 cdot1011 , m - the average distance from the Earth to the Sun. Substitute the real parameters in this formula

\ frac {F_ {1,3}} {F_ {1,2}} = \ frac {1.99 \ cdot 10 ^ {30} {5.98 \ cdot10 ^ {24}} \, \ left (\ frac {3.844 \ cdot10 ^ {8}} {1.496 \ cdot10 ^ {11} + 3.844 \ cdot10 ^ {8}} \ right) ^ 2 = 2.19

\ frac {F_ {1,3}} {F_ {1,2}} = \ frac {1.99 \ cdot 10 ^ {30} {5.98 \ cdot10 ^ {24}} \, \ left (\ frac {3.844 \ cdot10 ^ {8}} {1.496 \ cdot10 ^ {11} + 3.844 \ cdot10 ^ {8}} \ right) ^ 2 = 2.19


Here is the number! It turns out that the Moon is attracted to the Sun by a force that is more than two times greater than the force of its attraction to the Earth.

Such a disturbance can no longer be ignored, and it will definitely affect the final trajectory of the Moon. Let us go further, taking into account the assumption that the Earth’s orbit is circular with radius a, we find the locus of points around the Earth, where the force of attraction of any object to the Earth is equal to the force of its attraction to the Sun. This will be a sphere with a radius

R= fraca sqrt gamma1 gamma


displaced along a straight line connecting the Earth and the Sun in the opposite direction to the Sun at a distance

l=R sqrt gamma


Where  gamma=m2/m3 - The ratio of the mass of the earth to the mass of the sun. Substituting the numerical values ​​of the parameters, we obtain the actual dimensions of this area: R = 259300 kilometers, and l = 450 kilometers. This sphere is called the sphere of the Earth relative to the Sun.

The known orbit of the moon lies outside this area. That is, at any point of the trajectory, the Moon suffers from the Sun significantly more attraction than from the Earth.

2. Satellite or planet? Gravitational scope


This information often gives rise to disputes that the Moon is not a satellite of the Earth, but an independent planet of the Solar System, whose orbit is outraged by the attraction of the near Earth.

Let us estimate the perturbation introduced by the Sun into the trajectory of the Moon relative to the Earth, as well as the disturbance introduced by the Earth into the trajectory of the Moon 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 us assume that the orbits of the Earth relative to the Sun and the Moon relative to the Earth 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 system is determined by the forces acting on it and is equal to:

 veca1= veca(3)1+ veca(2)1= frac1m1 vecF1,3+ frac1m1 vecF1,2


On the other hand, according to the Coriolis theorem, the absolute acceleration of the moon

 veca1= veca2+ veca1,2


Where  veca2 - portable acceleration equal to the acceleration of the Earth relative to the Sun;  veca1,2 - acceleration of the moon relative to the earth. Accelerations of Coriolis will not be here - the chosen coordinate system moves progressively. Hence we get the acceleration of the moon relative to the earth.

 veca1,2= frac1m1 vecF1,3+ frac1m1 vecF1,2 veca2


A portion of this acceleration equal to  veca(2)1= frac1m1 vecF1,2 due to the attraction of the Moon to the Earth and characterizes its undisturbed geocentric movement. Remaining part

 Delta veca1,3= frac1m1 vecF1,3 veca2


the acceleration of the moon caused by disturbance from the sun.

If we consider the motion of the moon in the heliocentric inertial reference frame, then everything is much simpler, the acceleration  veca(3)1= frac1m1 vecF1,3 characterizes the undisturbed heliocentric motion of the moon, and the acceleration  Delta veca1,2= frac1m1 vecF1,2 - outrage of this movement from the Earth.

With the parameters of the orbits of the Earth and the Moon existing in the current epoch, the inequality is true at each point of the trajectory of the Moon

 frac| Delta veca1,3|| veca(2)1|< frac| Delta veca1,2|| veca(3)1| quad quad(1)


which can be verified by direct calculation, but I refer to the source in order not to overly clutter the article.

What does inequality (1) mean? Yes, the fact that in relative terms the effect of perturbation of the Moon by the Sun (and very significantly) is less than the effect of the attraction of the Moon to the Earth. And vice versa, the Earth's perturbation of the geoliocentric trajectory of the Moon has a decisive influence on the nature of its movement. The influence of the earth's gravity in this case is more significant, which means the moon "belongs" to the earth by right and is its satellite.

Another interesting thing is that by turning inequality (1) into an equation, you can find the locus of points where the effects of perturbations of the Moon (and any other body) by the Earth and the Sun are the same. Unfortunately, this is not so easy as in the case of the realm of the sphere. Calculations show that this surface is described by an equation of crazy order, but close to an ellipsoid of revolution. All that we can do without unnecessary problems is to assess the overall dimensions of this surface relative to the center of the Earth. Numerically solving the equation

 frac| Delta veca1,3|| veca(2)1|= frac| Delta veca1,2|| veca(3)1| quad quad(2)


relative to the distance from the center of the Earth to the desired surface on a sufficient number of points, we obtain the cross section of the desired surface by the ecliptic plane


For clarity, the geocentric orbit of the moon is shown here, and the sphere we found above is relative to the sun. It can be seen from the figure that the sphere of influence, or the sphere of the gravitational action of the Earth relative to the Sun, is the surface of rotation about the X axis, flattened along the straight line connecting the Earth and the Sun (along the axis of eclipses). The moon's orbit is deep inside this imaginary surface.

For practical calculations, this surface is conveniently approximated by a sphere with a center in the center of the Earth and a radius equal to

r=a left( fracmM right) frac25 quad quad(3)


where m is the mass of a smaller celestial body; M is the mass of the larger body, in the field of which the smaller body moves; a - the distance between the centers of bodies. In our case

r=a left( fracm2m3 right) frac25=1.496 cdot1011 left( frac5.98 cdot10241.99 cdot1030 right) frac25=925000,km



This unfinished million kilometers is the theoretical limit beyond which the power of the old woman of the Earth does not extend - its influence on the trajectories of astronomical objects is so small that it 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) does not work, it is physically impossible.

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

For example, it is now clear that in order to have a theoretical opportunity to perform maneuvers for entering a circumlunar orbit, the spacecraft must fall within the sphere of the Moon’s action relative to the Earth. Its radius is easy to calculate by the 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, which means its trajectory is not a conic section.

3. The three-body problem in the classical formulation


So, we will consider a model problem in a general formulation, known in celestial mechanics as a three-body problem. Consider three bodies of arbitrary mass, arranged arbitrarily in space and moving exclusively under the action of forces of mutual gravitational attraction


The bodies are considered material points. The position of the bodies will be counted in an arbitrary basis, which is associated with the inertial reference system Oxyz . The position of each of the bodies is given by the radius vector, respectively  vecr1 ,  vecr2 and  vecr3 . Each body is affected by the force of gravitational attraction from the other two bodies, and in accordance with the third axiom of the dynamics of a point (Newton's 3rd law)

 vecFi,j= vecFj,i quad quad(4)



We write the differential equations of motion of each point in vector form.

\ begin {align} & m_1 \, \ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ vec F_ {1,2} + \ vec F_ {1,3} \\ & m_2 \, \ frac {d ^ 2 \ vec r_2} {dt ^ 2} = \ vec F_ {2,1} + \ vec F_ {2,3} \\ & m_3 \, \ frac {d ^ 2 \ vec r_3} {dt ^ 2} = \ vec F_ {3,1} + \ vec F_ {3,2} \ end {align}



or, taking into account (4)

\ begin {align} & m_1 \, \ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ vec F_ {1,2} + \ vec F_ {1,3} \\ & m_2 \, \ frac {d ^ 2 \ vec r_2} {dt ^ 2} = - \ vec F_ {1,2} + \ vec F_ {2,3} \\ & m_3 \, \ frac {d ^ 2 \ vec r_3} { dt ^ 2} = - \ vec F_ {1,3} - \ vec F_ {2,3} \ end {align}


In accordance with the law of the world, the interaction forces are directed along the vectors

\ begin {align} & \ vec r_ {1,2} = \ vec r_2 - \ vec r_1 \\ & \ vec r_ {1,3} = \ vec r_3 - \ vec r_1 \\ & \ vec r_ {2 , 3} = \ vec r_3 - \ vec r_2 \\ \ end {align}


Along each of these vectors we will release the corresponding ort

 vecei,j= frac1ri,j vecri,j


then each of the gravitational forces is calculated by the formula

 vecFi,j=G fracmimjr2i,j vecei,j


Given all this, the system of equations of motion takes the form

\ begin {align} & \ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ frac {G \, m_2} {r_ {1,2} ^ 3} \, \ vec r_ {1,2 } + \ frac {G \, m_3} {r_ {1,3} ^ 3} \, \ vec r_ {1,3} \\ & \ frac {d ^ 2 \ vec r_2} {dt ^ 2} = - \ frac {G \, m_1} {r_ {1,2} ^ 3} \, \ vec r_ {1,2} + \ frac {G \, m_3} {r_ {2,3} ^ 3} \, \ vec r_ {2,3} \\ & \ frac {d ^ 2 \ vec r_3} {dt ^ 2} = - \ frac {G \, m_1} {r_ {1,3} ^ 3} \, \ vec r_ {1,3} - \ frac {G \, m_2} {r_ {2,3} ^ 3} \, \ vec r_ {2,3} \ end {align}


We introduce the notation adopted in celestial mechanics.

 mui=Gmi


- gravitational parameter of the attracting center. Then the equations of motion will take the final vector view.

\ begin {align} & \ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ frac {\ mu_2} {r_ {1,2} ^ 3} \, \ vec r_ {1,2} + \ frac {\ mu_3} {r_ {1,3} ^ 3} \, \ vec r_ {1,3} \\ & \ frac {d ^ 2 \ vec r_2} {dt ^ 2} = - \ frac {\ mu_1} {r_ {1,2} ^ 3} \, \ vec r_ {1,2} + \ frac {\ mu_3} {r_ {2,3} ^ 3} \, \ vec r_ {2,3} \ \ & \ frac {d ^ 2 \ vec r_3} {dt ^ 2} = - \ frac {\ mu_1} {r_ {1,3} ^ 3} \, \ vec r_ {1,3} - \ frac {\ mu_2} {r_ {2,3} ^ 3} \, \ vec r_ {2,3} \ end {align}



4. Rationing of equations to dimensionless variables


A fairly popular technique in mathematical modeling is to reduce differential equations and other relations describing the process to dimensionless phase coordinates and dimensionless time. 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 tasks. I leave the question of how justified this is in each solvable problem, but I agree that in this case, this approach is quite fair.

So, we introduce some abstract celestial body with a gravitational parameter  mu such that the orbital period of the satellite in an elliptical orbit with the major semiaxis a around it is equal T . All these quantities, by virtue of the laws of mechanics, are related by

T=2 pi left( fraca3 mu right) frac12


We introduce the replacement of parameters. For the position of the points of our system

 vecri=a vec xii


Where  vec xii - dimensionless radius vector of the i-th point;
for gravity body parameters

 mui= varkappai mu


Where  varkappai - dimensionless gravitational parameter of the i-th point;
for time

t=T tau


Where  tau - dimensionless time.

Now recalculate the acceleration points of the system through these dimensionless parameters. Let's apply direct double differentiation in time. For speeds

 vecvi= fracd vecridt=a fracd vec xiidt= fracaT fracd vec xiid tau= frac12 pi sqrt frac mua fracd vec xiid tau.


For acceleration

 vecai= fracd vecvidt= frac12 pi sqrt frac mua frac1dt left( fracd vec xiid tau right)= frac14 pi2 frac mua2 fracd2 vec xiid tau2



When substituting the obtained relations into the equations of motion, everything elegantly collapses into beautiful equations:

\ begin {align} & \ frac {d ^ 2 \ vec \ xi_1} {d \ tau ^ 2} = 4 \, \ pi ^ 2 \, \ varkappa_2 \, \ frac {\ vec \ xi_2 - \ vec \ xi_1} {| \ vec \ xi_2 - \ vec \ xi_1 | ^ 3} + 4 \, \ pi ^ 2 \, \ varkappa_3 \, \ frac {\ vec \ xi_3 - \ vec \ xi_1} {| \ vec \ xi_3 - \ vec \ xi_1 | ^ 3} \\ & \ frac {d ^ 2 \ vec \ xi_2} {d \ tau ^ 2} = -4 \, \ pi ^ 2 \, \ varkappa_1 \, \ frac {\ vec \ xi_2 - \ vec \ xi_1} {| \ vec \ xi_2 - \ vec \ xi_1 | ^ 3} + 4 \, \ pi ^ 2 \, \ varkappa_3 \, \ frac {\ vec \ xi_3 - \ vec \ xi_2} {| \ vec \ xi_3 - \ vec \ xi_2 | ^ 3} \ quad \ quad (5) \\ & \ frac {d ^ 2 \ vec \ xi_3} {d \ tau ^ 2} = -4 \, \ pi ^ 2 \, \ varkappa_1 \, \ frac {\ vec \ xi_3 - \ vec \ xi_1} {| \ vec \ xi_3 - \ vec \ xi_1 | ^ 3} - 4 \, \ pi ^ 2 \, \ varkappa_2 \, \ frac {\ vec \ xi_3 - \ vec \ xi_2} {| \ vec \ xi_3 - \ vec \ xi_2 | ^ 3} \ end {align}



This system of equations is still considered non-integrable in analytic functions. Why is it considered and is not? Because the successes of the theory of the function of a complex variable led to the fact that the general solution of the three-body problem did appear in 1912 — Karl Sundman found an algorithm for finding coefficients for infinite series with respect to a complex parameter, which theoretically is a general solution of the three-body problem. But ... for the application of Zundman's series in practical calculations with the accuracy required for them, it requires obtaining such a number of members of these series that this task is much greater than the capabilities of computers even today.

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

5. Calculation of the initial conditions: extract the source data


As I wrote earlier , before starting the numerical integration, you should take care to calculate the initial conditions for the problem to be solved. In the problem under consideration, the search for initial conditions turns into an independent subtask, since system (5) gives us nine second-order scalar equations, which, when going to normal form, Cauchy increases the order of the system by another 2 times. That is, we need to calculate as many as 18 parameters — the initial positions and components of the initial velocity of all points of the system. Where do we get data on the position of the celestial bodies of interest to us? We live in a world where man walked on the moon - naturally, humanity must have information about how this moon moves and where it is located.

That is, you say, you, dude, offer us to take thick astronomical reference books from the shelves, blow dust from them ... They didn’t guess! I propose to go for these data to those who actually walked on the moon, to NASA, namely in the Jet Propulsion Laboratory, Pasadena, California. Here is the JPL Horizonts web interface .

Here, having spent some time studying the interface, we will get all the data we need. Choose a date, for example, yes we do not care, but let it be July 27, 2018 UT 20:21. Just at that moment there was a complete phase of the lunar eclipse. The program will give us a huge footwoman

Full output for the ephemeris of the Moon on 07/27/2018 20:21 (origin of coordinates in the center of the Earth)
******************************************************************************* Revised: Jul 31, 2013 Moon / (Earth) 301 GEOPHYSICAL DATA (updated 2018-Aug-13): Vol. Mean Radius, km = 1737.53+-0.03 Mass, x10^22 kg = 7.349 Radius (gravity), km = 1738.0 Surface emissivity = 0.92 Radius (IAU), km = 1737.4 GM, km^3/s^2 = 4902.800066 Density, g/cm^3 = 3.3437 GM 1-sigma, km^3/s^2 = +-0.0001 V(1,0) = +0.21 Surface accel., m/s^2 = 1.62 Earth/Moon mass ratio = 81.3005690769 Farside crust. thick. = ~80 - 90 km Mean crustal density = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 k2 = 0.024059 Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 Rot. Rate, rad/s = 0.0000026617 Geometric Albedo = 0.12 Mean angular diameter = 31'05.2" Orbit period = 27.321582 d Obliquity to orbit = 6.67 deg Eccentricity = 0.05490 Semi-major axis, a = 384400 km Inclination = 5.145 deg Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.393142 beta (CA/B), x10^-4 = 6.310213 gamma (BA/C), x10^-4 = 2.277317 Perihelion Aphelion Mean Solar Constant (W/m^2) 1414+-7 1323+-7 1368+-7 Maximum Planetary IR (W/m^2) 1314 1226 1268 Minimum Planetary IR (W/m^2) 5.2 5.2 5.2 ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Aug 15 20:45:05 2018 Pasadena, USA / Horizons ******************************************************************************* Target body name: Moon (301) {source: DE431mx} Center body name: Earth (399) {source: DE431mx} Center-site name: BODY CENTER ******************************************************************************* Start time : AD 2018-Jul-27 20:21:00.0003 TDB Stop time : AD 2018-Jul-28 20:21:00.0003 TDB Step-size : 0 steps ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {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/J2000.0 Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ******************************************************************************* JDTDB XYZ VX VY VZ LT RG RR ******************************************************************************* $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06 VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05 LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06 $$EOE ******************************************************************************* Coordinate system description: Ecliptic and Mean Equinox of Reference Epoch Reference epoch: J2000.0 XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76) X-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch 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 : Jon.D.Giorgini@jpl.nasa.gov ******************************************************************************* 


Br-rr, what is it? Without panic, for someone who taught astronomy well at school, there is nothing to fear for mechanics and mathematics. So, the most important final search coordinates and components of the velocity of the moon.

 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06 VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05 LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06 $$EOE 

Yes, yes, yes, they are Cartesian! If you carefully read the entire footcloth, then 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 (the plane of the ecliptic) for the epoch J2000. The X axis is directed along the line of intersection of the equatorial plane of the Earth and the ecliptic to the vernal equinox.The Z axis looks in the direction of the north pole of the Earth perpendicular to the plane of the ecliptic. Well, the Y axis complements all this happiness to the right three vectors. By default, the units of measurement of coordinates are astronomical units (the clever ones from NASA also give the value of the autronomic unit in kilometers). Speed ​​units: astronomical units per day, day is assumed to be 86400 seconds. Full stuffing!

Similar information can be obtained for the Earth.

Full output of the Earth's ephemeris on 07/27/2018 20:21 (the origin of coordinates in the center of mass of the Solar system)
 ******************************************************************************* Revised: Jul 31, 2013 Earth 399 GEOPHYSICAL PROPERTIES (revised Aug 13, 2018): Vol. Mean Radius (km) = 6371.01+-0.02 Mass x10^24 (kg)= 5.97219+-0.0006 Equ. radius, km = 6378.137 Mass layers: Polar axis, km = 6356.752 Atmos = 5.1 x 10^18 kg Flattening = 1/298.257223563 oceans = 1.4 x 10^21 kg Density, g/cm^3 = 5.51 crust = 2.6 x 10^22 kg J2 (IERS 2010) = 0.00108262545 mantle = 4.043 x 10^24 kg g_p, m/s^2 (polar) = 9.8321863685 outer core = 1.835 x 10^24 kg g_e, m/s^2 (equatorial) = 9.7803267715 inner core = 9.675 x 10^22 kg g_o, m/s^2 = 9.82022 Fluid core rad = 3480 km GM, km^3/s^2 = 398600.435436 Inner core rad = 1215 km GM 1-sigma, km^3/s^2 = 0.0014 Escape velocity = 11.186 km/s Rot. Rate (rad/s) = 0.00007292115 Surface Area: Mean sidereal day, hr = 23.9344695944 land = 1.48 x 10^8 km Mean solar day 2000.0, s = 86400.002 sea = 3.62 x 10^8 km Mean solar day 1820.0, s = 86400.0 Moment of inertia = 0.3308 Love no., k2 = 0.299 Mean Temperature, K = 270 Atm. pressure = 1.0 bar Vis. mag. V(1,0) = -3.86 Volume, km^3 = 1.08321 x 10^12 Geometric Albedo = 0.367 Magnetic moment = 0.61 gauss Rp^3 Solar Constant (W/m^2) = 1367.6 (mean), 1414 (perihelion), 1322 (aphelion) ORBIT CHARACTERISTICS: Obliquity to orbit, deg = 23.4392911 Sidereal orb period = 1.0000174 y Orbital speed, km/s = 29.79 Sidereal orb period = 365.25636 d Mean daily motion, deg/d = 0.9856474 Hill's sphere radius = 234.9 ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Aug 15 21:16:21 2018 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 : AD 2018-Jul-27 20:21:00.0003 TDB Stop time : AD 2018-Jul-28 20:21:00.0003 TDB Step-size : 0 steps ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {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/J2000.0 Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ******************************************************************************* JDTDB XYZ VX VY VZ LT RG RR ******************************************************************************* $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05 VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07 LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05 $$EOE ******************************************************************************* Coordinate system description: Ecliptic and Mean Equinox of Reference Epoch Reference epoch: J2000.0 XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76) X-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch 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 : Jon.D.Giorgini@jpl.nasa.gov ******************************************************************************* 


Here, the barycenter (center of mass) of the Solar system is chosen as the origin of coordinates. The data we are interested in

 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05 VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07 LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05 $$EOE 

For the moon, we will need the coordinates and speed relative to the barycenter of the solar system, we can calculate them, or we can ask NASA to give us such data

Full output of the Moon's ephemeris on 07/27/2018 20:21 (the origin of coordinates in the center of mass of the Solar System)
 ******************************************************************************* Revised: Jul 31, 2013 Moon / (Earth) 301 GEOPHYSICAL DATA (updated 2018-Aug-13): Vol. Mean Radius, km = 1737.53+-0.03 Mass, x10^22 kg = 7.349 Radius (gravity), km = 1738.0 Surface emissivity = 0.92 Radius (IAU), km = 1737.4 GM, km^3/s^2 = 4902.800066 Density, g/cm^3 = 3.3437 GM 1-sigma, km^3/s^2 = +-0.0001 V(1,0) = +0.21 Surface accel., m/s^2 = 1.62 Earth/Moon mass ratio = 81.3005690769 Farside crust. thick. = ~80 - 90 km Mean crustal density = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 k2 = 0.024059 Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 Rot. Rate, rad/s = 0.0000026617 Geometric Albedo = 0.12 Mean angular diameter = 31'05.2" Orbit period = 27.321582 d Obliquity to orbit = 6.67 deg Eccentricity = 0.05490 Semi-major axis, a = 384400 km Inclination = 5.145 deg Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.393142 beta (CA/B), x10^-4 = 6.310213 gamma (BA/C), x10^-4 = 2.277317 Perihelion Aphelion Mean Solar Constant (W/m^2) 1414+-7 1323+-7 1368+-7 Maximum Planetary IR (W/m^2) 1314 1226 1268 Minimum Planetary IR (W/m^2) 5.2 5.2 5.2 ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Aug 15 21:19:24 2018 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 : AD 2018-Jul-27 20:21:00.0003 TDB Stop time : AD 2018-Jul-28 20:21:00.0003 TDB Step-size : 0 steps ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {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/J2000.0 Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ******************************************************************************* JDTDB XYZ VX VY VZ LT RG RR ******************************************************************************* $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05 VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05 LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05 $$EOE ******************************************************************************* Coordinate system description: Ecliptic and Mean Equinox of Reference Epoch Reference epoch: J2000.0 XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76) X-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch 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 : Jon.D.Giorgini@jpl.nasa.gov ******************************************************************************* 


 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05 VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05 LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05 $$EOE 

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

6. 38 parrots and one parrot wing


To begin with, we will define the scale, because our equations of motion (5) are written in dimensionless form. The data provided by NASA themselves tell us that for the scale of the coordinates is to take one astronomical unit. Accordingly, as a reference body, to which we will ration the masses of other bodies, we take the Sun, and as the time scale, the period of the Earth’s orbit around the Sun.

All this is of course very good, but we did not set the initial conditions for the sun. “Why?” Some linguist would ask me. And I would say that the Sun is not motionless at all, 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 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 6.520050993518213E+04 Y = 1.049687363172734E+06 Z =-1.304404963058507E+04 VX=-1.265326939350981E-02 VY= 5.853475278436883E-03 VZ= 3.136673455633667E-04 LT= 3.508397935601254E+00 RG= 1.051791240756026E+06 RR= 5.053500842402456E-03 $$EOE 

Looking at the RG parameter, we will see that the Sun rotates around the solar system barycenter, and as of 07/27/2018, the center of the star is at a distance of one 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 star. Why?Yes, because all other bodies interacting with the Sun also inform him of the acceleration, mainly, of course, the heavy Jupiter. Accordingly, the Sun also has its own orbit.

Of course, we can choose these data as initial conditions, but no - we solve the three-body model problem, and Jupiter and other characters are not included in it. So to the detriment of realism, knowing the position and velocity of the Earth and the Moon, we recalculate 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

(m1+m2+m3)rC=m1r1+m2r2+m3r3



We place the center of mass at the origin, that is, let us set rC=0 then

m1r1+m2r2+m3r3=0


from where

m3r3=m1r1m2r2r3=m1m3r1m2m3r2


We turn to the dimensionless coordinates and parameters by selecting μ=μ3

ξ3=ϰ1ξ1ϰ2ξ2(6)


Differentiating (6) over time and moving to dimensionless time, we also obtain the ratio for the velocities

u3=ϰ1u1ϰ2u2


Where ui=dξidτ,i=¯1,3

Now we will write a program that will form the initial conditions in the “parrots” we have chosen. What will we write on? Of course on Python! After all, as you know, this is the best language for mathematical modeling.

However, if we get away from sarcasm, then we will really try a python for this purpose, and why not? I will definitely provide a link to all the code in my Github profile .

Calculation of the initial conditions for the system Moon - Earth - Sun
 # #    # #   G = 6.67e-11 #   (, , ) m = [7.349e22, 5.792e24, 1.989e30] #     mu = [] print("  ") for i, mass in enumerate(m): mu.append(G * mass) print("mu[" + str(i) + "] = " + str(mu[i])) #      kappa = [] print("  ") for i, gp in enumerate(mu): kappa.append(gp / mu[2]) print("xi[" + str(i) + "] = " + str(kappa[i])) print("\n") #   a = 1.495978707e11 import math #   , c T = 2 * math.pi * a * math.sqrt(a / mu[2]) print("  T = " + str(T) + "\n") #  NASA   xL = 5.771034756256845E-01 yL = -8.321193799697072E-01 zL = -4.855790760378579E-05 import numpy as np xi_10 = np.array([xL, yL, zL]) print("  , ..: " + str(xi_10)) #  NASA   xE = 5.755663665315949E-01 yE = -8.298818915224488E-01 zE = -5.366994499016168E-05 xi_20 = np.array([xE, yE, zE]) print("  , ..: " + str(xi_20)) #    ,     -      xi_30 = - kappa[0] * xi_10 - kappa[1] * xi_20 print("  , ..: " + str(xi_30)) #       Td = 86400.0 u = math.sqrt(mu[2] / a) / 2 / math.pi print("\n") #    vxL = 1.434571674368357E-02 vyL = 9.997686898668805E-03 vzL = -5.149408819470315E-05 vL0 = np.array([vxL, vyL, vzL]) uL0 = np.array([0.0, 0.0, 0.0]) for i, v in enumerate(vL0): vL0[i] = v * a / Td uL0[i] = vL0[i] / u print("  , /: " + str(vL0)) print(" -//- : " + str(uL0)) #    vxE = 1.388633512282171E-02 vyE = 9.678934168415631E-03 vzE = 3.429889230737491E-07 vE0 = np.array([vxE, vyE, vzE]) uE0 = np.array([0.0, 0.0, 0.0]) for i, v in enumerate(vE0): vE0[i] = v * a / Td uE0[i] = vE0[i] / u print("  , /: " + str(vE0)) print(" -//- : " + str(uE0)) #    vS0 = - kappa[0] * vL0 - kappa[1] * vE0 uS0 = - kappa[0] * uL0 - kappa[1] * uE0 print("  , /: " + str(vS0)) print(" -//- : " + str(uS0)) 


Program exhaust

    mu[0] = 4901783000000.0 mu[1] = 386326400000000.0 mu[2] = 1.326663e+20    xi[0] = 3.6948215183509304e-08 xi[1] = 2.912016088486677e-06 xi[2] = 1.0   T = 31563683.35432583   , ..: [ 5.77103476e-01 -8.32119380e-01 -4.85579076e-05]   , ..: [ 5.75566367e-01 -8.29881892e-01 -5.36699450e-05]   , ..: [-1.69738146e-06 2.44737475e-06 1.58081871e-10]   , /: [24838.98933473 17310.56333294 -89.15979106] -//- : [ 5.24078311 3.65235907 -0.01881184]   , /: [2.40435899e+04 1.67586567e+04 5.93870516e-01] -//- : [5.07296163e+00 3.53591219e+00 1.25300854e-04]   , /: [-7.09330769e-02 -4.94410725e-02 1.56493465e-06] -//- : [-1.49661835e-05 -1.04315813e-05 3.30185861e-10] 

7. Integrating the equations of motion and analyzing the results


Actually, the integration itself boils down to a more or less standard for SciPy procedure for preparing a system of equations: transforming an ODE system to the Cauchy form and calling the corresponding solver functions. To transform the system into a Cauchy form, we recall that

ui=dξidτ,i=¯1,3(7)


Then by entering the system state vector

y=[ξ1,ξ2,ξ1,u1,u2,u3]T


we reduce (7) and (5) to one vector equation

d yd τ =f (τ,y )( 8 )


To integrate (8) with the initial conditions, we write a little, quite a bit of code

Integrating the equations of motion in the three-body problem
 # #     # def calcAccels(xi): k = 4 * math.pi ** 2 xi12 = xi[1] - xi[0] xi13 = xi[2] - xi[0] xi23 = xi[2] - xi[1] s12 = math.sqrt(np.dot(xi12, xi12)) s13 = math.sqrt(np.dot(xi13, xi13)) s23 = math.sqrt(np.dot(xi23, xi23)) a1 = (k * kappa[1] / s12 ** 3) * xi12 + (k * kappa[2] / s13 ** 3) * xi13 a2 = -(k * kappa[0] / s12 ** 3) * xi12 + (k * kappa[2] / s23 ** 3) * xi23 a3 = -(k * kappa[0] / s13 ** 3) * xi13 - (k * kappa[1] / s23 ** 3) * xi23 return [a1, a2, a3] # #       # def f(t, y): n = 9 dydt = np.zeros((2 * n)) for i in range(0, 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 = n for accel in accels: for a in accel: dydt[i] = a i = i + 1 return dydt #     y0 = [xi_10[0], xi_10[1], xi_10[2], xi_20[0], xi_20[1], xi_20[2], xi_30[0], xi_30[1], xi_30[2], uL0[0], uL0[1], uL0[2], uE0[0], uE0[1], uE0[2], uS0[0], uS0[1], uS0[2]] # #    # #   t_begin = 0 #   t_end = 30.7 * Td / T; #      N_plots = 1000 #     step = (t_end - t_begin) / N_plots import scipy.integrate as spi solver = spi.ode(f) solver.set_integrator('vode', nsteps=50000, method='bdf', max_step=1e-6, rtol=1e-12) solver.set_initial_value(y0, t_begin) ts = [] ys = [] i = 0 while 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 did. The result was a spatial trajectory of the moon for the first 29 days from the chosen starting point.


as well as its projection into the ecliptic plane.


“Hey, uncle, what are you giving us?” It's a circle! ”

First, it’s not a circle - there is a noticeable shift in the projection of the trajectory from the origin to the right and down. Secondly - you notice nothing? No, well, really?


I promise to prepare a rationale for that (based on the analysis of the errors of the account and NASA data) that the resulting shift in the trajectory is not due to integration errors. While I suggest the reader to take my word for it - this shift is a consequence of the solar perturbation of the lunar trajectory. Spin one more turn



In both! And pay attention to the fact that, based on the initial data of the task, the Sun is just in that direction, where the trajectory of the Moon is shifted on each turn. Yes, this insolent Sun steals our favorite companion! Oh, this is the sun!

It can be concluded that solar gravity affects the moon's orbit quite significantly - the old woman does not walk twice through the sky in the same way. The picture for six months of movement allows (at least qualitatively) to make sure of this (the picture is clickable)

image

Interesting?Still would. Astronomy is generally entertaining science.

P.S


In the university where I studied and worked for almost seven years - the Novocherkassk Polytechnic Institute - an annual student Olympiad on the theoretical mechanics of universities in the North Caucasus was held annually. Three times we took 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 have achieved in my life and career have happened thanks to this science and my wonderful teachers. I respect the mechanics.

Therefore, I will never allow anyone to scoff at this science and brazenly exploit it for their own purposes to anyone, even if it is at least three times a PhD and a linguist four times, and has developed at least a million training programs. I sincerely believe that writing articles on a popular public resource should provide for their careful reading, normal design (LaTeX formulas are not the whim of resource developers!) And the absence of errors leading to results that violate the laws of nature. The latter is generally a “must have”.

I often tell my students: “The computer frees your hands, but this does not mean that you need to turn off the brain as well.”

I also urge you, my dear readers, to appreciate and respect the mechanics. I will gladly answer any questions, and the source code of the example of solving the three-body problem in Python, as promised,I post it on my Github profile .

Thanks for attention!

Source: https://habr.com/ru/post/420133/


All Articles