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

     1  n =0 
     2  pi = 3.14159265358 
     3  solarMass = 4 *pi*pi 
     4  daysPerYear =365.24 
     5  N = 5 
     6  
     7   sysV = [ 
     8      [
     9          0.0,
    10          0.0, 
    11          0.0, 
    12          solarMass
    13      ],
    14      [
    15          1.66007664274403694e-03 *daysPerYear,
    16          7.69901118419740425e-03 *daysPerYear,
    17          -6.90460016972063023e-05* daysPerYear,
    18          9.54791938424326609e-04 *solarMass
    19      ],
    20      [
    21          -2.76742510726862411e-03*daysPerYear,
    22          4.99852801234917238e-03 *daysPerYear,
    23          2.30417297573763929e-05 *daysPerYear,
    24          2.85885980666130812e-04 *solarMass
    25      ],
    26      [
    27          2.96460137564761618e-03 *daysPerYear,
    28          2.37847173959480950e-03 *daysPerYear,
    29          -2.96589568540237556e-05* daysPerYear,
    30          4.36624404335156298e-05 *solarMass
    31      ],
    32      [
    33          2.68067772490389322e-03 *daysPerYear,
    34          1.62824170038242295e-03 *daysPerYear,
    35          -9.51592254519715870e-05* daysPerYear,
    36          5.15138902046611451e-05 *solarMass
    37      ]
    38  ]
    39  sysS = [
    40      [ 0.0,0.0,0.0 ],
    41      [ 4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01 ],
    42      [ 8.34336671824457987e+00, 4.12479856412430479e+00 , -4.03523417114321381e-01 ],
    43      [ 1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01 ],
    44      [ 1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01  ]
    45  ] 
    46  
    47  function offsetMomentum()
    48      px=0
    49      py=0
    50      pz = 0 
    51      for i = 0, N do
    52           m = sysV[i][3] 
    53          px = px + sysV[i][0] * m 
    54          py = py + sysV[i][1] * m 
    55          pz = pz + sysV[i][2] * m 
    56      end
    57  
    58      sysV[0][0] = (-px) / solarMass 
    59      sysV[0][1] = (-py) / solarMass 
    60      sysV[0][2] = (-pz) / solarMass 
    61  end
    62  
    63  sqrt = math.sqrt 
    64  function energy()
    65      e = 0
    66      for i=0, N do
    67          x = sysV[i][0]
    68          y = sysV[i][1]
    69          z = sysV[i][2] 
    70  
    71          e = e + sysV[i][3] * 0.5 * (x * x + y * y + z * z) 
    72  
    73          for j = i+1,N do
    74              dx = sysS[i][0] - sysS[j][0] 
    75              dy = sysS[i][1] - sysS[j][1] 
    76              dz = sysS[i][2] - sysS[j][2] 
    77  
    78              distance = sqrt(dx * dx + dy * dy + dz * dz) 
    79              e = e - sysV[i][3] * sysV[j][3] / distance 
    80          end
    81      end
    82      return e 
    83  end
    84  
    85  function advance(dt)
    86      for i=0, N - 1 do
    87          sysVi = sysV[i]
    88          sysSi = sysS[i]  
    89  
    90          _vx = sysVi[0] 
    91          _vy = sysVi[1] 
    92          _vz = sysVi[2] 
    93          for j = i+1, N do
    94              sysSj = sysS[j] 
    95              sysVj = sysV[j] 
    96  
    97              dx = sysSi[0] - sysSj[0] 
    98              dy = sysSi[1] - sysSj[1] 
    99              dz = sysSi[2] - sysSj[2] 
   100  
   101              dSquared = dx * dx + dy * dy + dz * dz 
   102              distance = sqrt(dSquared) 
   103              mag = dt / (dSquared * distance) 
   104  
   105              mi = sysVi[3] * mag
   106              m = -sysVj[3] * mag
   107  
   108              _vx = _vx + dx * m
   109              _vy = _vy + dy * m
   110              _vz = _vz + dz * m
   111  
   112              sysVj[0] = sysVj[0] + dx * mi
   113              sysVj[1] = sysVj[1] + dy * mi
   114              sysVj[2] = sysVj[2] + dz * mi
   115          end
   116  
   117          sysVi[0] = _vx
   118          sysVi[1] = _vy
   119          sysVi[2] = _vz
   120      end
   121  
   122      for i = 0, N do
   123          sysSi = sysS[i]
   124          sysVi = sysV[i]
   125          sysSi[0] = sysSi[0] + dt * sysVi[0]
   126          sysSi[1] = sysSi[1] + dt * sysVi[1]
   127          sysSi[2] = sysSi[2] + dt * sysVi[2]
   128      end
   129  end
   130  
   131  offsetMomentum() 
   132  println(energy()) 
   133  
   134  for i=0, 500000 do
   135      advance(0.01) 
   136  end
   137  
   138  println(energy())