nbody.go

Documentation: golang.org/x/exp/shootout

     1  //go:build ignore
     2  // +build ignore
     3  
     4  /*
     5  Redistribution and use in source and binary forms, with or without
     6  modification, are permitted provided that the following conditions are met:
     7  
     8      * Redistributions of source code must retain the above copyright
     9      notice, this list of conditions and the following disclaimer.
    10  
    11      * Redistributions in binary form must reproduce the above copyright
    12      notice, this list of conditions and the following disclaimer in the
    13      documentation and/or other materials provided with the distribution.
    14  
    15      * Neither the name of "The Computer Language Benchmarks Game" nor the
    16      name of "The Computer Language Shootout Benchmarks" nor the names of
    17      its contributors may be used to endorse or promote products derived
    18      from this software without specific prior written permission.
    19  
    20  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
    21  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
    22  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
    23  ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
    24  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
    25  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
    26  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
    27  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
    28  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    29  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
    30  POSSIBILITY OF SUCH DAMAGE.
    31  */
    32  
    33  /* The Computer Language Benchmarks Game
    34   * http://shootout.alioth.debian.org/
    35   *
    36   * contributed by The Go Authors.
    37   * based on C program by Christoph Bauer
    38   */
    39  
    40  package main
    41  
    42  import (
    43  	"flag"
    44  	"fmt"
    45  	"math"
    46  )
    47  
    48  var n = flag.Int("n", 1000, "number of iterations")
    49  
    50  type Body struct {
    51  	x, y, z, vx, vy, vz, mass float64
    52  }
    53  
    54  const (
    55  	solarMass   = 4 * math.Pi * math.Pi
    56  	daysPerYear = 365.24
    57  )
    58  
    59  func (b *Body) offsetMomentum(px, py, pz float64) {
    60  	b.vx = -px / solarMass
    61  	b.vy = -py / solarMass
    62  	b.vz = -pz / solarMass
    63  }
    64  
    65  type System []*Body
    66  
    67  func NewSystem(body []Body) System {
    68  	n := make(System, len(body))
    69  	for i := 0; i < len(body); i++ {
    70  		n[i] = new(Body) // copy to avoid overwriting the inputs
    71  		*n[i] = body[i]
    72  	}
    73  	var px, py, pz float64
    74  	for _, body := range n {
    75  		px += body.vx * body.mass
    76  		py += body.vy * body.mass
    77  		pz += body.vz * body.mass
    78  	}
    79  	n[0].offsetMomentum(px, py, pz)
    80  	return n
    81  }
    82  
    83  func (sys System) energy() float64 {
    84  	var e float64
    85  	for i, body := range sys {
    86  		e += 0.5 * body.mass *
    87  			(body.vx*body.vx + body.vy*body.vy + body.vz*body.vz)
    88  		for j := i + 1; j < len(sys); j++ {
    89  			body2 := sys[j]
    90  			dx := body.x - body2.x
    91  			dy := body.y - body2.y
    92  			dz := body.z - body2.z
    93  			distance := math.Sqrt(dx*dx + dy*dy + dz*dz)
    94  			e -= (body.mass * body2.mass) / distance
    95  		}
    96  	}
    97  	return e
    98  }
    99  
   100  func (sys System) advance(dt float64) {
   101  	for i, body := range sys {
   102  		for j := i + 1; j < len(sys); j++ {
   103  			body2 := sys[j]
   104  			dx := body.x - body2.x
   105  			dy := body.y - body2.y
   106  			dz := body.z - body2.z
   107  
   108  			dSquared := dx*dx + dy*dy + dz*dz
   109  			distance := math.Sqrt(dSquared)
   110  			mag := dt / (dSquared * distance)
   111  
   112  			body.vx -= dx * body2.mass * mag
   113  			body.vy -= dy * body2.mass * mag
   114  			body.vz -= dz * body2.mass * mag
   115  
   116  			body2.vx += dx * body.mass * mag
   117  			body2.vy += dy * body.mass * mag
   118  			body2.vz += dz * body.mass * mag
   119  		}
   120  	}
   121  
   122  	for _, body := range sys {
   123  		body.x += dt * body.vx
   124  		body.y += dt * body.vy
   125  		body.z += dt * body.vz
   126  	}
   127  }
   128  
   129  var (
   130  	jupiter = Body{
   131  		x:    4.84143144246472090e+00,
   132  		y:    -1.16032004402742839e+00,
   133  		z:    -1.03622044471123109e-01,
   134  		vx:   1.66007664274403694e-03 * daysPerYear,
   135  		vy:   7.69901118419740425e-03 * daysPerYear,
   136  		vz:   -6.90460016972063023e-05 * daysPerYear,
   137  		mass: 9.54791938424326609e-04 * solarMass,
   138  	}
   139  	saturn = Body{
   140  		x:    8.34336671824457987e+00,
   141  		y:    4.12479856412430479e+00,
   142  		z:    -4.03523417114321381e-01,
   143  		vx:   -2.76742510726862411e-03 * daysPerYear,
   144  		vy:   4.99852801234917238e-03 * daysPerYear,
   145  		vz:   2.30417297573763929e-05 * daysPerYear,
   146  		mass: 2.85885980666130812e-04 * solarMass,
   147  	}
   148  	uranus = Body{
   149  		x:    1.28943695621391310e+01,
   150  		y:    -1.51111514016986312e+01,
   151  		z:    -2.23307578892655734e-01,
   152  		vx:   2.96460137564761618e-03 * daysPerYear,
   153  		vy:   2.37847173959480950e-03 * daysPerYear,
   154  		vz:   -2.96589568540237556e-05 * daysPerYear,
   155  		mass: 4.36624404335156298e-05 * solarMass,
   156  	}
   157  	neptune = Body{
   158  		x:    1.53796971148509165e+01,
   159  		y:    -2.59193146099879641e+01,
   160  		z:    1.79258772950371181e-01,
   161  		vx:   2.68067772490389322e-03 * daysPerYear,
   162  		vy:   1.62824170038242295e-03 * daysPerYear,
   163  		vz:   -9.51592254519715870e-05 * daysPerYear,
   164  		mass: 5.15138902046611451e-05 * solarMass,
   165  	}
   166  	sun = Body{
   167  		mass: solarMass,
   168  	}
   169  )
   170  
   171  func main() {
   172  	flag.Parse()
   173  
   174  	system := NewSystem([]Body{sun, jupiter, saturn, uranus, neptune})
   175  	fmt.Printf("%.9f\n", system.energy())
   176  	for i := 0; i < *n; i++ {
   177  		system.advance(0.01)
   178  	}
   179  	fmt.Printf("%.9f\n", system.energy())
   180  }
   181
View as plain text