github.com/coyove/nj@v0.0.0-20221110084952-c7f8db1065c3/tests/bench/n-body.py (about)

     1  # The Computer Language Benchmarks Game
     2  # https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
     3  #
     4  # originally by Kevin Carson
     5  # modified by Tupteq, Fredrik Johansson, and Daniel Nanz
     6  # modified by Maciej Fijalkowski
     7  # 2to3
     8  
     9  import sys 
    10  
    11  def combinations(l):
    12      result = []
    13      for x in range(len(l) - 1):
    14          ls = l[x+1:]
    15          for y in ls:
    16              result.append((l[x],y))
    17      return result
    18  
    19  PI = 3.14159265358979323
    20  SOLAR_MASS = 4 * PI * PI
    21  DAYS_PER_YEAR = 365.24
    22  
    23  BODIES = {
    24      'sun': ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], SOLAR_MASS),
    25  
    26      'jupiter': ([4.84143144246472090e+00,
    27                   -1.16032004402742839e+00,
    28                   -1.03622044471123109e-01],
    29                  [1.66007664274403694e-03 * DAYS_PER_YEAR,
    30                   7.69901118419740425e-03 * DAYS_PER_YEAR,
    31                   -6.90460016972063023e-05 * DAYS_PER_YEAR],
    32                  9.54791938424326609e-04 * SOLAR_MASS),
    33  
    34      'saturn': ([8.34336671824457987e+00,
    35                  4.12479856412430479e+00,
    36                  -4.03523417114321381e-01],
    37                 [-2.76742510726862411e-03 * DAYS_PER_YEAR,
    38                  4.99852801234917238e-03 * DAYS_PER_YEAR,
    39                  2.30417297573763929e-05 * DAYS_PER_YEAR],
    40                 2.85885980666130812e-04 * SOLAR_MASS),
    41  
    42      'uranus': ([1.28943695621391310e+01,
    43                  -1.51111514016986312e+01,
    44                  -2.23307578892655734e-01],
    45                 [2.96460137564761618e-03 * DAYS_PER_YEAR,
    46                  2.37847173959480950e-03 * DAYS_PER_YEAR,
    47                  -2.96589568540237556e-05 * DAYS_PER_YEAR],
    48                 4.36624404335156298e-05 * SOLAR_MASS),
    49  
    50      'neptune': ([1.53796971148509165e+01,
    51                   -2.59193146099879641e+01,
    52                   1.79258772950371181e-01],
    53                  [2.68067772490389322e-03 * DAYS_PER_YEAR,
    54                   1.62824170038242295e-03 * DAYS_PER_YEAR,
    55                   -9.51592254519715870e-05 * DAYS_PER_YEAR],
    56                  5.15138902046611451e-05 * SOLAR_MASS) }
    57  
    58  
    59  SYSTEM = list(BODIES.values())
    60  PAIRS = combinations(SYSTEM)
    61  
    62  
    63  def advance(dt, n, bodies=SYSTEM, pairs=PAIRS):
    64  
    65      for i in range(n):
    66          for (([x1, y1, z1], v1, m1),
    67               ([x2, y2, z2], v2, m2)) in pairs:
    68              dx = x1 - x2
    69              dy = y1 - y2
    70              dz = z1 - z2
    71              mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))
    72              b1m = m1 * mag
    73              b2m = m2 * mag
    74              v1[0] -= dx * b2m
    75              v1[1] -= dy * b2m
    76              v1[2] -= dz * b2m
    77              v2[0] += dx * b1m
    78              v2[1] += dy * b1m
    79              v2[2] += dz * b1m
    80          for (r, [vx, vy, vz], m) in bodies:
    81              r[0] += dt * vx
    82              r[1] += dt * vy
    83              r[2] += dt * vz
    84  
    85  
    86  def report_energy(bodies=SYSTEM, pairs=PAIRS, e=0.0):
    87  
    88      for (((x1, y1, z1), v1, m1),
    89           ((x2, y2, z2), v2, m2)) in pairs:
    90          dx = x1 - x2
    91          dy = y1 - y2
    92          dz = z1 - z2
    93          e -= (m1 * m2) / ((dx * dx + dy * dy + dz * dz) ** 0.5)
    94      for (r, [vx, vy, vz], m) in bodies:
    95          e += m * (vx * vx + vy * vy + vz * vz) / 2.
    96      print("%.9f" % e)
    97  
    98  def offset_momentum(ref, bodies=SYSTEM, px=0.0, py=0.0, pz=0.0):
    99  
   100      for (r, [vx, vy, vz], m) in bodies:
   101          px -= vx * m
   102          py -= vy * m
   103          pz -= vz * m
   104      (r, v, m) = ref
   105      v[0] = px / m
   106      v[1] = py / m
   107      v[2] = pz / m
   108  
   109  def main(n, ref='sun'):
   110      offset_momentum(BODIES[ref])
   111      report_energy()
   112      advance(0.01, n)
   113      report_energy()
   114  
   115  if __name__ == '__main__':
   116      main(int(sys.argv[1]))