nbodys.odin

package main

import "core:fmt";
import "core:math";

FloatType :: f64;
Vec3 :: distinct [3]FloatType;

Body :: struct {
	pos: Vec3,
	v: Vec3,
	mass: FloatType,
}

System :: []^Body;

solarMass :: 4 * math.PI * math.PI;
daysPerYear :: 365.24;

offsetMomentum :: proc(b: ^Body, p: Vec3) {
	b.v = -p / solarMass;
};

NewSystem :: proc(body : []Body) -> System {
	n := make(System, len(body));
	for v, i in body {
		n[i] = &body[i];
	}
	p: Vec3;
	for body, _ in n {
		p += body.v * body.mass;
	}
	offsetMomentum(n[0], p);
	return n;
}

energy :: proc(sys : System) -> FloatType {
	e : FloatType;
	for body, i in sys {
		e += 0.5 * body.mass *
			(body.v.x*body.v.x + body.v.y*body.v.y + body.v.z*body.v.z);
		for j in i+1..len(sys)-1 {
			body2 := sys[j];
			d := body.pos - body2.pos;
			distance := math.sqrt(d.x*d.x + d.y*d.y + d.z*d.z);
			e -= (body.mass * body2.mass) / distance;
		}
	}
	return e;
}

advance :: proc(sys: ^System, dt: FloatType) {
	l := len(sys) - 1;

	for body, i in sys {
		for j in i+1..l {
			body2 := sys[j];
			d := body.pos - body2.pos;

			dSquared := d.x*d.x + d.y*d.y + d.z*d.z;
			distance := math.sqrt(dSquared);
			mag := dt / (dSquared * distance);

			body.v -= d * body2.mass * mag;

			body2.v += d * body.mass * mag;
		}
	}

	for body, _ in sys {
		body.pos += dt * body.v;
	}
}

main :: proc() {
	n := 50000000;

	Jupiter := Body{
		pos = {4.84143144246472090e+00,
		-1.16032004402742839e+00,
		-1.03622044471123109e-01},
		v = {1.66007664274403694e-03 * daysPerYear,
		7.69901118419740425e-03 * daysPerYear,
		-6.90460016972063023e-05 * daysPerYear},
		mass=9.54791938424326609e-04 * solarMass,
	};

	Saturn := Body{
		pos={8.34336671824457987e+00,
		4.12479856412430479e+00,
		-4.03523417114321381e-01,},
		v={-2.76742510726862411e-03 * daysPerYear,
		4.99852801234917238e-03 * daysPerYear,
		2.30417297573763929e-05 * daysPerYear,},
		mass=2.85885980666130812e-04 * solarMass,
	};

	Uranus := Body{
		pos={1.28943695621391310e+01,
		-1.51111514016986312e+01,
		-2.23307578892655734e-01,},
		v={2.96460137564761618e-03 * daysPerYear,
		2.37847173959480950e-03 * daysPerYear,
		-2.96589568540237556e-05 * daysPerYear,},
		mass=4.36624404335156298e-05 * solarMass,
	};

	Neptune := Body{
		pos={1.53796971148509165e+01,
		-2.59193146099879641e+01,
		1.79258772950371181e-01,},
		v={2.68067772490389322e-03 * daysPerYear,
		1.62824170038242295e-03 * daysPerYear,
		-9.51592254519715870e-05 * daysPerYear,},
		mass=5.15138902046611451e-05 * solarMass,
	};

	Sun := Body{
		mass=solarMass
	};

	system := NewSystem([]Body{Sun, Jupiter, Saturn, Uranus, Neptune});
	fmt.printf("%.9f\n", energy(system));
	for i in 1..n {
		advance(&system, 0.01);
	}
	fmt.printf("%.9f\n", energy(system));
}