nbody.go

package nbody

import (
	"encoding/csv"
	"math"
	"os"
	"strconv"
	"strings"
)

type Bodies []Body

type Body struct {
	Pos [3]float64
	Vel [3]float64
	GM float64
}

func EOM (y []float64, dy []float64, t float64, bodies Bodies) {
	var dr = [3]float64{}
	var rm2, rm3 float64

	for bod1 := range bodies {
		for i := 0; i < 3; i++ {
			dy[bod1*6+i] = y[bod1*6+i+3] // dx = v
			dy[bod1*6+i+3] = 0
			if math.IsNaN(dy[bod1*6+i]){
				panic("boom2")
			}
		}
		for bod2 := range bodies {
			if bod1 == bod2 {
				continue
			}
			rm2 = 0
			for i := 0; i < 3; i++ {
				dr[i] = y[bod1*6+i] - y[bod2*6+i]
				rm2 += dr[i] * dr[i]
			}
			rm3 = math.Sqrt(rm2)
			rm3 *= rm2 // r.magnitude^3
			if rm3 == 0 {
				panic("boom")
			}
			for i := 0; i < 3; i++ {
				dy[bod1*6+i+3] -= bodies[bod2].GM * dr[i] / rm3
				if math.IsNaN(dy[bod1*6+i+3]){
					panic("boom3")
				}
			}
		}
	}
}

func BuildInitial(bodies Bodies) []float64 {
	y0 := make([]float64, len(bodies) * 6)
	for bod := range bodies {
		for i := 0; i < 3; i++ {
			y0[bod*6+i] = bodies[bod].Pos[i]
			y0[bod*6+i+3] = bodies[bod].Vel[i]
		}
	}
	return  y0
}

func BodiesFromCSV(csvfile string) (Bodies, error) {
	file, err := os.Open(csvfile)
	if err != nil {
		return nil, err
	}
	reader := csv.NewReader(file)
	records, _ := reader.ReadAll()
	var bodies Bodies
	for i := range records {
		var gm float64
		var pos [3]float64
		var vel [3]float64

		if gm, err = strconv.ParseFloat(strings.TrimSpace(records[i][1]), 64); err != nil {
			return nil, err
		}

		for j := 0; j < 3; j++ {
			var p, v float64
			if p, err = strconv.ParseFloat(strings.TrimSpace(records[i][4+j]), 64); err != nil {
				return nil, err
			}
			pos[j] = p
			if v, err = strconv.ParseFloat(strings.TrimSpace(records[i][7+j]), 64); err != nil {
				return nil, err
			}
			vel[j] = v
		}

		body := Body { Pos: pos, Vel: vel, GM: gm}
		bodies = append(bodies, body)
	}
	return bodies, nil
}