F 1 , 2 = G f r a c m 1m 2 r 2 1 , 2
F 1 , 3 = G f r a c m 1m 3 r 2 1 , 3
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
fracF1,3F1,2= fracm3m2 left( frac rhoa+ rho right)2
\ 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
R= fraca sqrt gamma1− gamma
l=R sqrt gamma
veca1= veca(3)1+ veca(2)1= frac1m1 vecF1,3+ frac1m1 vecF1,2
veca1= veca2+ veca1,2
veca1,2= frac1m1 vecF1,3+ frac1m1 vecF1,2− veca2
Delta veca1,3= frac1m1 vecF1,3− veca2
frac| Delta veca1,3|| veca(2)1|< frac| Delta veca1,2|| veca(3)1| quad quad(1)
frac| Delta veca1,3|| veca(2)1|= frac| Delta veca1,2|| veca(3)1| quad quad(2)
r=a left( fracmM right) frac25 quad quad(3)
r=a left( fracm2m3 right) frac25=1.496 cdot1011 left( frac5.98 cdot10241.99 cdot1030 right) frac25=925000,km
vecFi,j=− vecFj,i quad quad(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_ {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}
\ 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}
\ 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}
vecei,j= frac1ri,j vecri,j
vecFi,j=G fracmimjr2i,j vecei,j
\ 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}
mui=Gmi
\ 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}
T=2 pi left( fraca3 mu right) frac12
vecri=a vec xii
mui= varkappai mu
t=T tau
vecvi= fracd vecridt=a fracd vec xiidt= fracaT fracd vec xiid tau= frac12 pi sqrt frac mua fracd vec xiid tau.
vecai= fracd vecvidt= frac12 pi sqrt frac mua frac1dt left( fracd vec xiid tau right)= frac14 pi2 frac mua2 fracd2 vec xiid tau2
\ 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}
******************************************************************************* 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 *******************************************************************************
$$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
******************************************************************************* 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 *******************************************************************************
$$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
******************************************************************************* 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
$$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
(m1+m2+m3)→rC=m1→r1+m2→r2+m3→r3
m1→r1+m2→r2+m3→r3=0
m3→r3=−m1→r1−m2→r2→r3=−m1m3→r1−m2m3→r2
→ξ3=−ϰ1→ξ1−ϰ2→ξ2(6)
→u3=−ϰ1→u1−ϰ2→u2
# # # # 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))
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]
→ui=d→ξidτ,∀i=¯1,3(7)
→y=[→ξ1,→ξ2,→ξ1,→u1,→u2,→u3]T
d → yd τ = → f (τ, → y )( 8 )
# # # 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
Source: https://habr.com/ru/post/420133/
All Articles