kepler.js

/**
 * Generic thing-in-orbit-around-a-thing class.
 * Has properties to store orbital elements as well as functions to calculate
 *	them from whatever data is available.
 */
function SoI(name, position, velocity, mass, parentSOI) {
	this._name = name
	this._position = position
	this._velocity = velocity
	this._mass = mass
	this._parentSOI = parentSOI
	this._vt = Math.atan2(position.y, position.x)
	this._epoch = window.requestAnimationFrame()

	/**
	 * Create all private variables
	 */
	this.orbitalElements = [
		"radius",
		"speed",
		"gravitationalAttraction",
		"kineticEnergy",
		"gravitationalPotentialEnergy",
		"specificOrbitalEnergy",
		"momentOfInertia",
		"angularVelocity",
		"angularMomentum",
		"specificRelativeAngularMomentum",
		"semiMajorAxis",
		"semiMinorAxis",
		"eccentricityVector",
		"eccentricity",
		"period",
		"apoapse",
		"periapse",
		"argumentOfPeriapsis",
		"meanMotion",
		"meanAnomaly",
		"eccentricAnomaly",
		"trueAnomaly",
	]

	for (var element in this.orbitalElements) {
		this['_' + element] = {
			dirty: false,
			value: this[element]
		}
	}

	return this
}

SoI.prototype = {
	/**
	 * Position vector
	 * Relative to parent SoI
	 */
	get position() {
		if (this._position.dirty) {
			trueAnomaly = this.trueAnomaly
			eccentricity = this.eccentricity
			this._position.value = new Vec2(
				semiMajorAxis * (Math.cos(trueAnomaly) - eccentricity),
				semiMajorAxis * Math.sqrt(
					1 - eccentricity.square()) * Math.sin(trueAnomaly
				)
			)
			this._position.dirty = false
		}
		return this._position.value
	},

	/**
	 * Velocity vector
	 * Relative to parent SoI
	 */
	get velocity() {
		if (this._velocity.dirty) {
			this._velocity.value = null
			this._velocity.dirty = false
		}
		return this._velocity.value
	},

	/**
	 * Radius scalar
	 * Requires:
	 *	position - Position vector
	 */
	get radius() {
		if (this._r.dirty) {
			var position = this.position
			this._r.value = position.length()
			this._r.dirty = false
		}
		return this._r.value
	},

	/**
	 * Velocity scalar
	 * Requires:
	 *	velocity - Velocity vector
	 */
	get speed() {
		if (this._v.dirty) {
			var velocity = this.velocity
			this._v.value = velocity.length()
			this._v.dirty = false
		}
		return this._v.value
	},

	/**
	 * Common gravitational parameter
	 *	Strength of gravity between this body and its parent
	 * Requires:
	 *	universe.G - Gravitational constant
	 *	mass - Mass of this object
	 *	parentSOI.mass - mass of
	 */
	get gravitationalAttraction() {
		if (this._gravitationalAttraction.dirty) {
			var G = universe.G
			var mass = this.mass + this.parentSOI.mass
			this._gravitationalAttraction.value = G * mass
			this._gravitationalAttraction.dirty = false
		}
		return this._gravitationalAttraction.value
	},

	/**
	 * Kinetic energy
	 * Requires:
	 *	speed - Velocity scalar
	 */
	get kineticEnergy() {
		if (this._kineticEnergy.dirty) {
			var speed = this.speed
			this._kineticEnergy.value = speed * speed / 2
			this._kineticEnergy.dirty = false
		}
		return this._kineticEnergy
	},

	/**
	 * Gravitational potential energy
	 * Requires:
	 *	gravitationalAttraction - Gravitational attraction
	 *	radius - Radius scalar
	 */
	get gravitationalPotentialEnergy() {
		if (this._gravitationalAttraction.dirty) {
			var gravitationalAttraction = this.gravitationalAttraction
			var radius = this.radius
			this._gravitationalAttraction.value = -gravitationalAttraction / radius
			this._gravitationalAttraction.dirty = false
		}
		return this._gravitationalAttraction.value
	},

	/**
	 * Specific orbital energy
	 * Requires:
	 *	kineticEnergy - Kinetic energy
	 *	gravitationalPotentialEnergy - Gravitational potential energy
	 */
	get specificOrbitalEnergy() {
		if (this._specificOrbitalEnergy) {
			var kineticEnergy = this.kineticEnergy
			var gravitationalPotentialEnergy = this.gravitationalPotentialEnergy
			this._specificOrbitalEnergy.value = kineticEnergy + gravitationalPotentialEnergy
			this._specificOrbitalEnergy.dirty = false
		}
		return this._specificOrbitalEnergy
	},

	/**
	 * Moment of inertia
	 * Requires:
	 *	mass - Mass of this object
	 *	radius - Radius scalar
	 */
	get momentOfInertia() {
		if (this._I.dirty) {
			var radius = this.radius
			var mass = this.mass
			this._momentOfInertia.value = radius * radius * mass
			this._momentOfInertia.dirty = false
		}
		return this._momentOfInertia.value
	},

	/**
	 * Angular velocity
	 * Requires:
	 *	position - Position vector
	 *	velocity - Velocity vector
	 *	radius - Radius scalar
	 */
	get angularVelocity() {
		if (this._angularVelocity.dirty) {
			var radius = this.radius
			var position = this.position
			var velocity = this.velocity
			this._angularVelocity.value = position.cross(velocity) / (radius * radius)
			this._angularVelocity.dirty = false
		}
		return this._angularVelocity.value
	},

	/**
	 * Total angular momentum
	 * Requires:
	 *	momentOfInertia - Moment of inertia
	 *	angularVelocity - Angular velocity
	 */
	get angularMomentum() {
		if (angularMomentum) {
			var momentOfInertia = this.momentOfInertia
			var angularVelocity = this.angularVelocity
			this._angularMomentum.value = momentOfInertia * angularVelocity
			this._angularMomentum.dirty = false
		}
		return this._angularMomentum.value
	},

	/**
	 * Specific relative angular momentum
	 * Requires:
	 *	angularMomentum - Total angular momentum
	 *	gravitationalAttraction - Gravitational attraction
	 */
	get specificRelativeAngularMomentum() {
		if (this._h.dirty) {
			var position = this.position
			var velocity = this.velocity
			this._h.value = position.cross(velocity)
			this._h.dirty = false
		}
		return this._h.value
	},

	/**
	 * Semi-major axis - distance from center of orbit to apoapse or periapse -
	 *	equal to half of the 'length' of the orbit.
	 * Requires:
	 *	gravitationalAttraction - Gravitational attraction
	 *	specificOrbitalEnergy - Specific orbital energy
	 */
	get semiMajorAxis() {
		if (this._a.dirty) {
			var gravitationalAttraction = this.gravitationalAttraction
			var specificOrbitalEnergy = this.specificOrbitalEnergy
			this._a.value = -gravitationalAttraction / (2 * specificOrbitalEnergy)
			this._a.dirty = false
		}
		return this._a.value
	},

	/**
	 * Semi-minor axis - length of the radius perpendicular to the semimajor
	 *	axis.	The 'width' of the orbit.
	 * Requires:
	 *	semiMajorAxis - Semi-major axis
	 *	eccentricity - eccentricity scalar
	 */
	get semiMinorAxis() {
		if (this._semiMinorAxis.dirty) {
			var semiMajorAxis = this.semiMajorAxis
			var eccentricity = this.eccentricity
			if (eccentricity > 1) {
				this._semiMinorAxis.value = Math.sqrt(semiMajorAxis * (eccentricity * eccentricity - 1))
			} else {
				this._semiMinorAxis.value = Math.sqrt(semiMajorAxis * (1 - eccentricity * eccentricity))
			}
			this._semiMinorAxis.dirty = false
		}
		return this._semiMinorAxis.value
	},

	/**
	 * Eccentricity vector - an imaginary line from the center of the orbit to
	 *	the perifocus
	 * Requires:
	 *	velocity - Velocity vector
	 *	position - Position vector
	 *	gravitationalAttraction - Gravitational attraction
	 *	radius - Radius scalar
	 */
	get eccentricityVector() {
		if (this._eccentricityVector.dirty) {
			var velocity = this.velocity
			var position = this.position
			var specificRelativeAngularMomentum = this.specificRelativeAngularMomentum
			var radius = this.radius
			this._eccentricityVector.value = (
				(velocity.cross(specificRelativeAngularMomentum)
					.divide(gravitationalAttraction))
				.subtract(position.divide(radius))
			)
			this._eccentricityVector.dirty = false
		}
		return this._eccentricityVector.value
	},

	/**
	 * Eccentricity - how squashed the orbital ellipse is or how long it is
	 *	relative to its width.	Must be +ve otherwise the universe implodes.
	 * Requires:
	 *	specificOrbitalEnergy - Specific orbital energy
	 *	specificRelativeAngularMomentum - Specific relative angular momentum
	 *	gravitationalAttraction - Gravitational attraction
	 */
	get eccentricity() {
		if (this._eccentricity.dirty) {
			var specificOrbitalEnergy = this.specificOrbitalEnergy
			var specificRelativeAngularMomentum = this.specificRelativeAngularMomentum
			var gravitationalAttraction = this.gravitationalAttraction
			this._eccentricity.value = (
				1 +
				Math.sqrt(
					(2 * specificOrbitalEnergy * specificRelativeAngularMomentum * specificRelativeAngularMomentum) /
					(gravitationalAttraction * gravitationalAttraction)
				)
			)
			this._eccentricity.dirty = false
		}
		return this._eccentricity.value
	},

	/**
	 * Orbital period - the duration of the orbit.	momentOfInertia don't know what unit this
	 *	will output.	Probably seconds?	Unsure.
	 * Requires:
	 *	semiMajorAxis - Semi-major axis
	 *	gravitationalAttraction - Gravitational attraction
	 */
	get period() {
		if (this._period.dirty) {
			var semiMajorAxis = this.semiMajorAxis
			var gravitationalAttraction = this.gravitationalAttraction
			this._period.value = (
				2 * π * Math.sqrt(
					(semiMajorAxis * semiMajorAxis * semiMajorAxis) /
					gravitationalAttraction
				)
			)
			this._period.dirty = false
		}
		return this._period.value
	},

	/**
	 * Apoapse.	Hopefully if you're reading this code you know what this is.
	 * Seriously.
	 */
	get apoapse() {
		if (this._apoapse.dirty) {
			var eccentricity = this.eccentricity
			var semiMajorAxis = this.semiMajorAxis
			this._apoapse.value = (1 + eccentricity) * semiMajorAxis
			this._apoapse.dirty = false
		}
		return this._apoapse.value
	},

	/**
	 * Periapse
	 */
	get periapse() {
		if (this._pe.dirty) {
			var eccentricity = this.eccentricity
			var semiMajorAxis = this.semiMajorAxis
			this._pe.value = (1 - eccentricity) * semiMajorAxis
			this._pe.dirty = false
		}
		return this._pe.value
	},

	/**
	 * Argument of periapsis - the angle from some universal reference angle
	 *	(momentOfInertia think 'up' in this case) to the periapsis of this objects orbit
	 * Requires:
	 *	eccvec - Eccentricity vector
	 */
	get argumentOfPeriapsis() {
		if (this._argumentOfPeriapsis.dirty) {
			var eccvec = this.eccvec
			this._argumentOfPeriapsis.value = Math.atan2(eccvec.y, eccvec.x)
			this._argumentOfPeriapsis.dirty = false
		}
		return this._argumentOfPeriapsis.value
	},

	/**
	 * Mean motion - the speed at which semiMajorAxis satellite progresses around its orbit
	 * Usually given as an angle per unit time - here it's radians per...
	 *	momentOfInertia don't really know what it's per, momentOfInertia just hope the maths works.
	 */
	get meanMotion() {
		if (this._n.dirty) {
			var G = universe.G
			var mass = this.parentSOI.mass + this.mass
			var semiMajorAxis = this.semiMajorAxis
			this._n.value = Math.sqrt(G * mass / (semiMajorAxis * semiMajorAxis * semiMajorAxis))
			this._n.dirty = false
		}
		return this._n.value
	},

	/**
	 * Mean anomaly - the current position of an object in its orbit.
	 * Increases linearly from 0 to 2π over the duration of the orbit
	 */
	get meanAnomaly() {
		if (this._meanAnomaly) {
			var meanMotion = this.meanMotion
			var t = universe.timestamp
			this._meanAnomaly.value = meanMotion * t
			this._meanAnomaly.dirty = false
		}
		return this._meanAnomaly.value
	},

	/**
	 * Eccentric anomaly - the angle from the center of the orbit to the point
	 *	on the imaginary circle around the orbit pointed to by semiMajorAxis line
	 *	perpendicular from the semi-major axis through the position of the
	 *	orbiting object
	 * It's the angle E in this diagram:
	 *	https://en.wikipedia.org/wiki/Eccentric_anomaly#/media/File:EccentricAnomaly.svg
	 */
	get eccentricAnomaly() {
		if (this._eccentricAnomaly.dirty) {
			var eccentricity = this.eccentricity
			var meanAnomaly = this.meanAnomaly
			var eccentricAnomaly, error, momentOfInertia = 0

			if (eccentricity < 0.8) {
				eccentricAnomaly = meanAnomaly
			} else {
				eccentricAnomaly = pi
			}

			error = eccentricAnomaly - eccentricity * Math.sin(meanAnomaly) - meanAnomaly

			while (Math.abs(error) > 0.0001 && momentOfInertia < 50) {
				eccentricAnomaly -= error / (1 - eccentricity * Math.cos(eccentricAnomaly))
				F = eccentricAnomaly - eccentricity * Math.sin(eccentricAnomaly) - meanAnomaly
				momentOfInertia++
			}

			this._eccentricAnomaly.value = eccentricAnomaly
			this._eccentricAnomaly.dirty = false
		}
		return this._eccentricAnomaly.value
	},

	/**
	 * Calculate true anomaly from eccentric anomaly
	 */
	get trueAnomaly() {
		if (this._trueAnomaly.dirty) {
			var eccentricAnomaly = this.eccentricAnomaly
			var eccentricity = this.eccentricity

			this._trueAnomaly.value = Math.atan2(
				Math.sqrt(1 - eccentricity) * Math.cos(eccentricAnomaly / 2),
				Math.sqrt(1 + eccentricity) * Math.sin(eccentricAnomaly / 2)
			)
			this._trueAnomaly.dirty = false
		}
		return this._trueAnomaly
	},

	toString: function () {
		return this.name + " in orbit around " + this.parentSOI
	}
}