github.com/gopherd/gonum@v0.0.4/spatial/barneshut/galaxy_example_test.go (about)

     1  // Copyright ©2019 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package barneshut_test
     6  
     7  import (
     8  	"log"
     9  
    10  	"math/rand"
    11  
    12  	"github.com/gopherd/gonum/spatial/barneshut"
    13  	"github.com/gopherd/gonum/spatial/r2"
    14  )
    15  
    16  type mass struct {
    17  	d r2.Vec
    18  	v r2.Vec
    19  	m float64
    20  }
    21  
    22  func (m *mass) Coord2() r2.Vec { return m.d }
    23  func (m *mass) Mass() float64  { return m.m }
    24  func (m *mass) move(f r2.Vec) {
    25  	m.v = m.v.Add(f.Scale(1 / m.m))
    26  	m.d = m.d.Add(m.v)
    27  }
    28  
    29  func Example_galaxy() {
    30  	rnd := rand.New(rand.NewSource(1))
    31  
    32  	// Make 1000 stars in random locations.
    33  	stars := make([]*mass, 1000)
    34  	p := make([]barneshut.Particle2, len(stars))
    35  	for i := range stars {
    36  		s := &mass{
    37  			d: r2.Vec{
    38  				X: 100 * rnd.Float64(),
    39  				Y: 100 * rnd.Float64(),
    40  			},
    41  			v: r2.Vec{
    42  				X: rnd.NormFloat64(),
    43  				Y: rnd.NormFloat64(),
    44  			},
    45  			m: 10 * rnd.Float64(),
    46  		}
    47  		stars[i] = s
    48  		p[i] = s
    49  	}
    50  	vectors := make([]r2.Vec, len(stars))
    51  
    52  	// Make a plane to calculate approximate forces
    53  	plane := barneshut.Plane{Particles: p}
    54  
    55  	// Run a simulation for 100 updates.
    56  	for i := 0; i < 1000; i++ {
    57  		// Build the data structure. For small systems
    58  		// this step may be omitted and ForceOn will
    59  		// perform the naive quadratic calculation
    60  		// without building the data structure.
    61  		err := plane.Reset()
    62  		if err != nil {
    63  			log.Fatal(err)
    64  		}
    65  
    66  		// Calculate the force vectors using the theta
    67  		// parameter...
    68  		const theta = 0.5
    69  		// and an imaginary gravitational constant.
    70  		const G = 10
    71  		for j, s := range stars {
    72  			vectors[j] = plane.ForceOn(s, theta, barneshut.Gravity2).Scale(G)
    73  		}
    74  
    75  		// Update positions.
    76  		for j, s := range stars {
    77  			s.move(vectors[j])
    78  		}
    79  
    80  		// Rendering stars is left as an exercise for
    81  		// the reader.
    82  	}
    83  }