Skip to content

Instantly share code, notes, and snippets.

@grom358
Created June 7, 2014 07:39
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save grom358/367295eeb2c8efa10e69 to your computer and use it in GitHub Desktop.
Save grom358/367295eeb2c8efa10e69 to your computer and use it in GitHub Desktop.
"use strict";
/**
* High performance Vector library.
*
* Constructing vectors is expensive and therefore these functions take
* the required operands and a parameter to output the result into.
*/
var Vector = {};
Vector.create = function(x, y) {
return new Float64Array([x, y]);
};
Vector.copy = function(a, b) {
b[0] = a[0];
b[1] = a[1];
};
Vector.add = function(a, b, out) {
out[0] = a[0] + b[0];
out[1] = a[1] + b[1];
};
Vector.scale = function(a, scale, out) {
out[0] = a[0] * scale;
out[1] = a[1] * scale;
};
Vector.length = function(a) {
return Math.sqrt(Math.pow(a[0], 2) + Math.pow(a[1], 2));
};
Vector.lengthSquared = function(a) {
return Math.pow(a[0], 2) + Math.pow(a[1], 2);
};
Vector.angle = function(a) {
return Math.atan(a[1] / a[0]);
};
Vector.normalize = function(a, out) {
var mag = 1 / Vector.length(a);
out[0] = a[0] * mag;
out[1] = a[1] * mag;
};
Vector.lengthTo = function(a, scale, out) {
var mag = scale / Vector.length(a);
out[0] = a[0] * mag;
out[1] = a[1] * mag;
};
Vector.fromPolar = function(angle, length, out) {
out[0] = length * Math.cos(angle);
out[1] = length * Math.sin(angle);
};
// Gravitational constant.
var G = 6.67384e-11;
/**
* RK4 Integrator.
*/
function Simulator(dt) {
this.t = 0;
this.dt = dt;
// Allocate memory for immediate calculations.
this.pos = Vector.create(0, 0);
this.a = {velocity: Vector.create(0, 0), acceleration: Vector.create(0, 0)};
this.b = {velocity: Vector.create(0, 0), acceleration: Vector.create(0, 0)};
this.c = {velocity: Vector.create(0, 0), acceleration: Vector.create(0, 0)};
this.d = {velocity: Vector.create(0, 0), acceleration: Vector.create(0, 0)};
this.result = Vector.create(0, 0);
}
Simulator.prototype.evaluate = function(initialState, t, dt, derivative, output) {
Vector.scale(derivative.velocity, dt, this.result);
Vector.add(initialState.position, this.result, this.pos);
Vector.scale(derivative.acceleration, dt, this.result);
Vector.add(initialState.velocity, this.result, output.velocity);
this.calculateAcceleration(initialState, this.pos, t + dt, output.acceleration);
};
Simulator.prototype.integrate = function(state, t, dt) {
Vector.copy(state.velocity, this.a.velocity);
this.calculateAcceleration(state, state.position, t, this.a.acceleration);
this.evaluate(state, t, dt * 0.5, this.a, this.b);
this.evaluate(state, t, dt * 0.5, this.b, this.c);
this.evaluate(state, t, dt, this.c, this.d);
// position = position + 1/6 * dt * (a.velocity + 2 * (b.velocity + c.velocity) + d.velocity)
Vector.add(this.b.velocity, this.c.velocity, this.result);
Vector.scale(this.result, 2.0, this.result);
Vector.add(this.result, this.d.velocity, this.result);
Vector.add(this.result, this.a.velocity, this.result);
Vector.scale(this.result, dt / 6.0, this.result);
Vector.add(state.position, this.result, state.position);
// velocity = velocity + 1/6 * dt * (a.acceleration + 2 * (b.acceleration + c.acceleration) + d.acceleration)
Vector.add(this.b.acceleration, this.c.acceleration, this.result);
Vector.scale(this.result, 2.0, this.result);
Vector.add(this.result, this.d.acceleration, this.result);
Vector.add(this.result, this.a.acceleration, this.result);
Vector.scale(this.result, dt / 6.0, this.result);
Vector.add(state.velocity, this.result, state.velocity);
};
Simulator.prototype.update = function(state) {
this.integrate(state, this.t, this.dt);
this.t += this.dt;
};
Simulator.prototype.calculateAcceleration = function(state, pos, t, out) {
var rSquared = Vector.lengthSquared(pos);
var gravityForce = G * state.parentBody.mass / rSquared;
Vector.lengthTo(pos, -gravityForce, out);
};
function toRadians(x) {
return x * Math.PI / 180;
}
function toDegrees(x) {
return x * 180 / Math.PI;
}
function testCannon(h, angle, velocity, timeCap) {
var planet = {
radius: 6378100,
mass: 5.97219e24
};
// Useful information about orbit from h
console.log('Firing cannon from h:' + h);
var vo = Math.sqrt(G * planet.mass / (planet.radius + h));
var circumference = 2 * Math.PI * (planet.radius + h);
var orbitPeriod = circumference / vo;
console.log('=Horizontal velocity required=');
console.log('Orbit velocity: ' + vo);
console.log('Orbit period: ' + orbitPeriod);
console.log('------------------------------');
// Cannonball.
var state = {
parentBody: planet,
position: Vector.create(0, planet.radius + h),
velocity: Vector.create(0, 0)
};
Vector.fromPolar(angle, velocity, state.velocity);
// Setup simulator for cannon.
var orbitLoops = 0;
var lastOrbit = 0;
var orbitPeriod = 0;
var dt = 0.01;
var simulator = new Simulator(dt);
var previous = {
position: Vector.create(0, 0)
};
var ap = h;
var pe = h;
var alt = h;
// Run simulator until cannonball hits the group or achieves orbit.
while (simulator.t < timeCap && alt > 0) {
Vector.copy(state.position, previous.position);
simulator.update(state);
alt = Vector.length(state.position) - planet.radius;
pe = Math.max(0, Math.min(pe, alt));
ap = Math.max(ap, alt);
// Each time we pass the positive y axis the cannon has completed an orbit.
if (previous.position[0] < 0 && previous.position[1] > 0 && state.position[0] >= 0 && state.position[1] > 0) {
orbitLoops++;
orbitPeriod = simulator.t - lastOrbit;
}
}
if (alt <= 0) {
console.log('Impacted ground :(');
console.log('Impact t=' + simulator.t);
}
if (orbitLoops > 0) {
console.log('Orbit reached x' + orbitLoops);
console.log('Orbit period: ' + orbitPeriod);
}
console.log('pe=' + pe + ' ap=' + ap);
}
var h = 100000; // Height launching cannon ball from.
var angle = toRadians(30);
var velocity = 3000;
var timeCap = 20000;
testCannon(h, angle, velocity, timeCap);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment