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

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