added comparison of Newton and Yukawa potential

Author: azimonti

README.md

@@ -28,6 +28,10 @@ To get started with these simulations:
    python solar_system.py
    ```
 
+## Yukawa Potential
+
+Yukawa potential (a potential which has an exponential decay in addition to the standard Newtonian gravitational force) is also simulated between the Sun and the Earth, but the change in orbit is minimal (no change up to 7 significant digits).
+
 ## Contributing
 
 Contributions to the Orbits Simulation project are welcome. Whether it's through submitting bug reports, proposing new features, or contributing to the code, your help is appreciated. For major changes, please open an issue first to discuss what you would like to change.

body_data_extraction.py

@@ -1,3 +1,11 @@
+#!/usr/bin/env python3
+'''
+/***************************/
+/* body_data_extraction.py */
+/*       Version 1.0       */
+/*      2024/08/06         */
+/***************************/
+'''
 import json
 import re
 import requests

data_input/bodies_2Dplane.json

@@ -37,7 +37,6 @@
     "position": [ 1.496e11, 0.0, 0.0 ],
     "velocity": [ 0.0, 2.978e4, 0.0 ]
   },
-
   "mars": {
     "label": "Mars",
     "mass": 6.4171e23,

data_input/bodies_3D.json

@@ -62,7 +62,6 @@
       1.184224839788195
     ]
   },
-
   "mars": {
     "label": "Mars",
     "mass": 6.417099999999999e+23,

data_input/bodies_yukawa.json

@@ -0,0 +1,22 @@
+{
+  "bodies" : ["sun", "earth"],
+  "reference_date": "2024-01-01 0:00:00",
+  "sun": {
+    "label": "Sun",
+    "mass": 1.989e+30,
+    "scale": 1.5,
+    "color": [1, 1, 0],
+    "radius": 6.96342e+8,
+    "position": [ 0.0, 0.0, 0.0 ],
+    "velocity": [ 0.0, 0.0, 0.0 ]
+  },
+  "earth": {
+    "label": "Earth",
+    "mass": 5.97219e+24,
+    "scale": 2.7,
+    "color": [0, 0, 1],
+    "radius": 6371010.0,
+    "position": [ -24810993259.6539, 144994861273.6719, -8215203.670851886 ],
+    "velocity": [ -29841.463655186788, -5126.262286859616, 1.184224839788195 ]
+  }
+}

solar_system.py

@@ -1,3 +1,11 @@
+#!/usr/bin/env python3
+'''
+/**********************/
+/*  solar_system.py   */
+/*    Version 1.0     */
+/*     2024/08/03     */
+/**********************/
+'''
 import argparse
 import json
 import matplotlib.pyplot as plt
@@ -18,6 +26,7 @@
     animation_format='mp4',   # animation format (mp4 or gif)
     fps=30,                   # frame per second
     useRK45=True,             # use RK45
+    step=86400,               # integration step 1 day
     verbose=True,             # verbose output
     use_yukawa=False,         # use Yukawa potential
     yukawa_coeff=1.5e15,      # Yukawa potential coefficient
@@ -28,8 +37,7 @@
 
 
 cst = SimpleNamespace(
-    G=6.67430e-11,
-    day_in_seconds=86400
+    G=6.67430e-11
 )
 
 
@@ -67,7 +75,7 @@ def runge_kutta4(f, y0, t0, tf, dt, args=()):
 
 
 def gravitational_force_newton(m1, m2, r_ij):
-    """Compute the gravitational force between two masses using Newtwon
+    """Compute the gravitational force between two masses using Newton
     gravitational law."""
     r = np.linalg.norm(r_ij)
     return cst.G * m1 * m2 * r_ij / r**3
@@ -77,8 +85,8 @@ def gravitational_force_yukawa(m1, m2, r_ij):
     """Compute the gravitational force between two masses using Yukawa
     gravitational law."""
     r = np.linalg.norm(r_ij)
-    return cst.G * m1 * m2 * r_ij / r**3 * (1 + r_ij / cfg.yukawa_coeff) * \
-        np.exp(-r_ij / cfg.yukawa_coeff)
+    return cst.G * m1 * m2 * r_ij / r**3 * (1 + r / cfg.yukawa_coeff) * \
+        np.exp(-r / cfg.yukawa_coeff)
 
 
 gravitational_force = gravitational_force_yukawa if cfg.use_yukawa else \
@@ -113,8 +121,8 @@ class SolarSystemSimulation:
     def __init__(self, outfile):
         self._bodies = []
         self._num_frames = 30
-        self.outfile = outfile
-        self.perc = None
+        self._outfile = outfile
+        self._perc = None
 
     @property
     def bodies(self):
@@ -132,7 +140,7 @@ def create_bodies(self):
                      position=b.position, velocity=b.velocity))
 
     def compute(self):
-        t_eval = np.arange(0, cfg.time_span[1], cst.day_in_seconds)
+        t_eval = np.arange(0, cfg.time_span[1], cfg.step)
         # extract masses
         masses = np.array([body.mass for body in self._bodies])
         # flatten initial conditions using the `position` #
@@ -146,15 +154,15 @@ def compute(self):
             solution = solve_ivp(n_body_problem, cfg.time_span,
                                  initial_conditions, args=(masses,),
                                  method='RK45', t_eval=t_eval,
-                                 max_step=cst.day_in_seconds)
+                                 max_step=cfg.step)
             self._positions = solution.y[:num_bodies * 3].reshape(
                 (num_bodies, 3, -1))
             self._velocities = solution.y[num_bodies * 3:].reshape(
                 (num_bodies, 3, -1))
         else:
             t, y = runge_kutta4(n_body_problem, initial_conditions,
                                 cfg.time_span[0], cfg.time_span[1],
-                                cst.day_in_seconds, args=(masses,))
+                                cfg.step, args=(masses,))
             self._positions = y[:num_bodies * 3].reshape((num_bodies, 3, -1))
             self._velocities = y[num_bodies * 3:].reshape((num_bodies, 3, -1))
 
@@ -202,7 +210,7 @@ def format_func(value, tick_number):
         plt.tight_layout()
 
     def animate(self):
-        self.perc = 0
+        self._perc = 0
         if cfg.save_anim:
             fig = plt.figure(figsize=(12, 8), dpi=300)
         else:
@@ -245,8 +253,8 @@ def animate(self):
         def animate_frame(frame):
             if cfg.verbose:
                 perc = (frame + 1) / self._num_frames * 100
-                if perc // 10 > self.perc // 10:
-                    self.perc = perc
+                if perc // 10 > self._perc // 10:
+                    self._perc = perc
                     print(f"completed {int(perc)}% of the animation")
             for i in range(len(self._bodies)):
                 # Extract the positions for body i at the current frame
@@ -264,7 +272,7 @@ def animate_frame(frame):
             fig, animate_frame, frames=self._num_frames,
             interval=1000 / cfg.fps, blit=True)
         if cfg.save_anim:
-            base, ext = self.outfile.rsplit('.', 1)
+            base, ext = self._outfile.rsplit('.', 1)
             animation_format = cfg.animation_format
             outfile_a = f"{base}.{animation_format}"
             if animation_format == 'mp4':

solar_system_body.py

@@ -1,3 +1,11 @@
+#!/usr/bin/env python3
+'''
+/************************/
+/* solar_system_body.py */
+/*    Version 1.0       */
+/*     2024/08/03       */
+/************************/
+'''
 import numpy as np
 import sys