Coverage for pySDC/implementations/problem_classes/FullSolarSystem.py: 100%

1import numpy as np 

2 

3from pySDC.implementations.datatype_classes.particles import particles, acceleration 

4from pySDC.implementations.problem_classes.OuterSolarSystem import outer_solar_system 

5 

6 

7# noinspection PyUnusedLocal 

8class full_solar_system(outer_solar_system): 

9 r""" 

10 Gravitational N-body problem of the sun, the eight planets (Earth and Moon as one body) and Pluto. 

11 

12 The :math:`N`-body problem describes the mutual influence of the motion of :math:`N` bodies. Formulation of the problem is 

13 based on Newton's second law. Therefore, the :math:`N`-body problem is formulated as 

14 

15 .. math:: 

16 m_i \frac{d^2 {\bf r}_i}{d t^2} = \sum_{j=1, i\neq j}^N G \frac{m_i m_j}{|{\bf r}_i - {\bf r}_j|^3}({\bf r}_j - {\bf r}_i), 

17 

18 where :math:`m_i` is the :math:`i`-th mass point with position described by the vector :math:`{\bf r}_i`, and :math:`G` 

19 is the gravitational constant. If only the sun (body 1) influences the motion of the bodies gravitationally, the 

20 equations become 

21 

22 .. math:: 

23 \frac{d^2 {\bf r}_i}{d t^2} = G \frac{m_1}{|{\bf r}_i - {\bf r}_1|^3}({\bf r}_1 - {\bf r}_i), \quad i \geq 2, 

24 

25 and the sun is not accelerated. 

26 

27 This class implements the full solar system: the sun, the eight planets, with Earth and its moon combined into one 

28 body, and Pluto, i.e., :math:`N=10`. Initial conditions are taken from [1]_, and masses relative to the sun taken 

29 from [2]_. 

30 

31 Parameters 

32 ---------- 

33 sun_only : bool, optional 

34 If True, only the gravity of the sun acts on the other bodies, see above. If False (default), all bodies 

35 attract each other. 

36 

37 Attributes 

38 ---------- 

39 G : float 

40 Gravitational constant. 

41 

42 References 

43 ---------- 

44 .. [1] https://www.aanda.org/articles/aa/full/2002/08/aa1405/aa1405.right.html 

45 .. [2] https://en.wikipedia.org/wiki/Planetary_mass#Values_from_the_DE405_ephemeris 

46 """ 

47 

48 dtype_u = particles 

49 dtype_f = acceleration 

50 

51 def __init__(self, sun_only=False): 

52 """Initialization routine""" 

53 

54 super().__init__(sun_only) 

55 self.init = ((3, 10), None, np.dtype('float64')) 

56 

57 def u_exact(self, t): 

58 r""" 

59 Routine to compute the exact/initial trajectory at time :math:`t`. Values here are taken from [1]_, [2]_. 

60 

61 Parameters 

62 ---------- 

63 t : float 

64 Time of the exact/initial trajectory. 

65 

66 Returns 

67 ------- 

68 me : dtype_u 

69 Exact/initial position and velocity. 

70 """ 

71 assert t == 0.0, 'error, u_exact only works for the initial time t0=0' 

72 me = self.dtype_u(self.init) 

73 

74 me.pos[:, 0] = [0.0, 0.0, 0.0] 

75 me.pos[:, 1] = [-2.503321047836e-01, +1.873217481656e-01, +1.260230112145e-01] 

76 me.pos[:, 2] = [+1.747780055994e-02, -6.624210296743e-01, -2.991203277122e-01] 

77 me.pos[:, 3] = [-9.091916173950e-01, +3.592925969244e-01, +1.557729610506e-01] 

78 me.pos[:, 4] = [+1.203018828754e00, +7.270712989688e-01, +3.009561427569e-01] 

79 me.pos[:, 5] = [+3.733076999471e00, +3.052424824299e00, +1.217426663570e00] 

80 me.pos[:, 6] = [+6.164433062913e00, +6.366775402981e00, +2.364531109847e00] 

81 me.pos[:, 7] = [+1.457964661868e01, -1.236891078519e01, -5.623617280033e00] 

82 me.pos[:, 8] = [+1.695491139909e01, -2.288713988623e01, -9.789921035251e00] 

83 me.pos[:, 9] = [-9.707098450131e00, -2.804098175319e01, -5.823808919246e00] 

84 

85 me.vel[:, 0] = [0.0, 0.0, 0.0] 

86 me.vel[:, 1] = [-2.438808424736e-02, -1.850224608274e-02, -7.353811537540e-03] 

87 me.vel[:, 2] = [+2.008547034175e-02, +8.365454832702e-04, -8.947888514893e-04] 

88 me.vel[:, 3] = [-7.085843239142e-03, -1.455634327653e-02, -6.310912842359e-03] 

89 me.vel[:, 4] = [-7.124453943885e-03, +1.166307407692e-02, +5.542098698449e-03] 

90 me.vel[:, 5] = [-5.086540617947e-03, +5.493643783389e-03, +2.478685100749e-03] 

91 me.vel[:, 6] = [-4.426823593779e-03, +3.394060157503e-03, +1.592261423092e-03] 

92 me.vel[:, 7] = [+2.647505630327e-03, +2.487457379099e-03, +1.052000252243e-03] 

93 me.vel[:, 8] = [+2.568651772461e-03, +1.681832388267e-03, +6.245613982833e-04] 

94 me.vel[:, 9] = [+3.034112963576e-03, -1.111317562971e-03, -1.261841468083e-03] 

95 

96 me.m[0] = 1.0 # Sun 

97 me.m[1] = 0.1660100 * 1e-06 # Mercury 

98 me.m[2] = 2.4478383 * 1e-06 # Venus 

99 me.m[3] = 3.0404326 * 1e-06 # Earth+Moon 

100 me.m[4] = 0.3227151 * 1e-06 # Mars 

101 me.m[5] = 954.79194 * 1e-06 # Jupiter 

102 me.m[6] = 285.88600 * 1e-06 # Saturn 

103 me.m[7] = 43.662440 * 1e-06 # Uranus 

104 me.m[8] = 51.513890 * 1e-06 # Neptune 

105 me.m[9] = 0.0073960 * 1e-06 # Pluto 

106 

107 return me