Initial Commit with "Backwards" Orbit class.
Initial Commit with "Backwards" Orbit class.

file:b/Body.py (new)
--- /dev/null
+++ b/Body.py
@@ -1,1 +1,15 @@
-
+#!/usr/bin/python3

+

+# Gravitational Constant in m³/(kg×s²)

+G = 6.67384e-11

+

+class Body():

+	# Mass in kg, Radius in m.

+	def __init__(self, Name, Mass = None, Radius = None):

+		self.Name = Name

+		self.Mass = Mass

+		self.Radius = Radius

+	

+	# Gravitational parameter µ in m³/s²

+	def GravParm(self):

+		return self.Mass * G

file:b/Orbit.py (new)
--- /dev/null
+++ b/Orbit.py
@@ -1,1 +1,172 @@
-
+#!/usr/bin/python3

+

+from Body import Body

+from Vector3D import *

+from math import acos, atan2, cos, pi, sin

+

+class Orbit():

+	def __init__(self, PosVector = None, VelVector = None, Primary = None):

+		'''

+		@param PosVector:

+		@type Vector3D:

+		@param VelVector:

+		@type Vector3D:

+		@param Primary:

+		@type Body:

+		'''

+		self.PosVector = PosVector

+		self.VelVector = VelVector

+		self.Primary = Primary

+	

+	def VisVivaEnergy(self):

+		return self.VelVector.Magnitude() ** 2. / 2 - self.Primary.GravParm() / self.PosVector.Magnitude()

+	

+	def SemiMajorAxis(self):

+		return -self.Primary.GravParm() / (2. * self.VisVivaEnergy())

+	

+	def SpecAngMmntm(self):

+		return self.PosVector * self.VelVector

+	

+	def EccentricityVector(self):

+		return self.VelVector * self.SpecAngMmntm() / self.Primary.GravParm() - self.PosVector.UnitVector()

+	

+	def Eccentricity(self):

+		return self.EccentricityVector().Magnitude()

+	

+	def Inclination(self):

+		h = self.SpecAngMmntm()

+		hZ = h._k

+		return acos(hZ / h.Magnitude())

+	

+	def AscNodeVector(self):

+		return kVector * self.SpecAngMmntm()

+	

+	def LongitudeAscNode(self):

+		n = self.AscNodeVector()

+		

+		try:

+			Omega = acos(n._i / n.Magnitude())

+		except ZeroDivisionError:

+			return None

+		if n._j >= 0:

+			return Omega

+		else:

+			return 2 * pi - Omega

+	

+	def ArgPeriapsis(self):

+		e = self.EccentricityVector()

+		try:

+			n = self.AscNodeVector().UnitVector()

+		except ZeroDivisionError:

+			arg = atan2(e._j / e.Magnitude(), e._i / e.Magnitude())

+			if self.SpecAngMmntm()._k < 0:

+				arg = 2. * pi - arg

+			return arg

+		

+		arg = acos(n.Dot(e) / (n.Magnitude() * e.Magnitude()))

+		

+		if e._k >= 0:

+			return arg

+		else:

+			return 2 * pi - arg

+	

+	def TrueAnomaly(self):

+		e = self.EccentricityVector()

+		r = self.PosVector

+		f = acos(e.Dot(r) / (e.Magnitude() * r.Magnitude()))

+		if r.Dot(self.VelVector) >= 0:

+			return f

+		else:

+			return 2 * pi - f

+	

+	def EccentricAnomaly(self):

+		e = self.Eccentricity()

+		f = self.TrueAnomaly()

+		E = acos((e + cos(f)) / (1 + e * cos(f)))

+		if f <= pi:

+			return E

+		else:

+			return 2 * pi - E

+	

+	def MeanAnomaly(self):

+		E = self.EccentricAnomaly()

+		e = self.Eccentricity()

+		M = E - e * sin(E)

+		return M

+	

+	def PeriapsisVector(self):

+		dummy = Orbit(self.PosVector, self.VelVector, self.Primary)

+		dummy.UpdateTrueAnomaly(0)

+		

+		return dummy.PosVector

+	def Periapsis(self):

+		return self.PeriapsisVector().Magnitude()

+	

+	def ApoapsisVector(self):

+		dummy = Orbit(self.PosVector, self.VelVector, self.Primary)

+		dummy.UpdateTrueAnomaly(pi)

+		

+		return dummy.PosVector

+	def Apoapsis(self):

+		return self.ApoapsisVector().Magnitude()

+	

+	def UpdateFromElements(self, a, e, i, Omega, w, f):

+		p = (a * (1 - e ** 2.) / (1 + e * cos(f)))

+		Rx = p * (cos(Omega) * cos(w + f) - sin(Omega) * cos(i) * sin(w + f))

+		Ry = p * (sin(Omega) * cos(w + f) + cos(Omega) * cos(i) * sin(w + f))

+		Rz = p * sin(i) * sin(w + f)

+		PosVector = Vector3D(Rx, Ry, Rz)

+		

+		Vx = -(self.Primary.GravParm() / p) ** .5 * (cos(Omega) * (sin(w + f) + e * sin(w)) + sin(Omega) * cos(i) * (cos(w + f) + e * cos(w)))

+		Vy = -(self.Primary.GravParm() / p) ** .5 * (sin(Omega) * (sin(w + f) + e * sin(w)) - cos(Omega) * cos(i) * (cos(w + f) + e * cos(w)))

+		Vz = -(self.Primary.GravParm() / p) ** .5 * (sin(i) * (cos(w + f) + e * cos(w)))

+		VelVector = Vector3D(Vx, Vy, Vz)

+		

+		self.PosVector = PosVector

+		self.VelVector = VelVector

+	

+	def UpdateFromShapeAndPos(self, a, e, i, Omega, w, PosVector):

+		self.UpdateFromElements(a, e, i, Omega, w, 0)

+		self.PosVector = PosVector

+		

+# 		f = self.TrueAnomaly()

+# 		self.UpdateFromElements(a, e, i, Omega, w, f)

+# 		

+# 		if (self.PosVector - PosVector).Magnitude() > .001 * self.SemiMajorAxis():

+# 			s  = "Specified PosVector is not on specified orbit!\n"

+# 			s += "Specified PosVector: {0:.2}\n".format(PosVector)

+# 			s += "Calculated PosVector: {0:.2}\n".format(self.PosVector)

+# 			raise ValueError(s)

+	

+	def UpdateTrueAnomaly(self, f):

+		return self.UpdateFromElements(self.SemiMajorAxis(), self.Eccentricity(), self.Inclination(), self.LongitudeAscNode(), self.ArgPeriapsis(), f)

+	

+	def __str__(self):

+		s  = "Orbiting {0}\n".format(self.Primary.Name)

+		s += "a = {0} [m]\n".format(self.SemiMajorAxis())

+		s += "e = {0}\n".format(self.Eccentricity())

+		s += "i = {0} [radians]\n".format(self.Inclination())

+		s += "Ω = {0} [radians]\n".format(self.LongitudeAscNode())

+		s += "ω = {0} [radians]\n".format(self.ArgPeriapsis())

+		s += "υ = {0} [radians]\n".format(self.TrueAnomaly())

+		s += "\n"

+		s += "Current Position: {0} (radius: {1})\n".format(self.PosVector, self.PosVector.Magnitude())

+		s += "Current Velocity: {0} (magnitude: {1})\n".format(self.VelVector, self.VelVector.Magnitude())

+		

+		return s

+

+if __name__ == "__main__":

+	Kerbin = Body("Kerbin", Mass = 5.2915793e22, Radius = 600e3)

+	Orbit1 = Orbit(Primary = Kerbin)

+	Orbit1.UpdateFromElements(4.35e6, 0.28, pi / 6., 0, 3. * pi / 2., pi / 2.)

+	print(Orbit1)

+	

+	Orbit2 = Orbit(Primary = Kerbin)

+	Orbit2.UpdateFromShapeAndPos(Orbit1.PosVector.Magnitude(), 0., 0, 0, 0, Orbit1.PosVector)

+	print(Orbit2)

+	

+	deltaV = Orbit2.VelVector - Orbit1.VelVector

+	print("delta-V to transfer: {0} (magnitude: {1})".format(deltaV, deltaV.Magnitude()))

+	

+	dVPRN = deltaV.ConvertToFrame(Orbit1.VelVector, Orbit1.PosVector, Orbit1.VelVector * Orbit1.PosVector)

+	print("delta-V in Prograde, Radial, Normal: {0} (magnitude {1})".format(dVPRN, dVPRN.Magnitude()))

file:b/Vector3D.py (new)
--- /dev/null
+++ b/Vector3D.py
@@ -1,1 +1,167 @@
-
+#!/usr/bin/python3

+

+from numbers import Number

+from math import acos

+

+__all__ = ['Vector3D', 'iVector', 'jVector', 'kVector']

+

+class Vector3D:

+	def __init__(self, i=0, j=0, k=0):

+		self._i = i

+		self._j = j

+		self._k = k

+	

+	@classmethod

+	def CrossProduct(cls, uVector, vVector):

+		i = uVector._j * vVector._k - uVector._k * vVector._j

+		j = uVector._k * vVector._i - uVector._i * vVector._k

+		k = uVector._i * vVector._j - uVector._j * vVector._i 

+		return cls(i, j, k)

+	

+	def LeftCross(self, vVector):

+		return self.CrossProduct(self, vVector)

+	

+	def RightCross(self, uVector):

+		return self.CrossProduct(uVector, self)

+	

+	@classmethod

+	def DotProduct(cls, uVector, vVector):

+		return uVector._i * vVector._i + uVector._j * vVector._j + uVector._k * vVector._k

+	

+	def Dot(self, vVector):

+		return self.DotProduct(self, vVector)

+	

+	def ScalarMultiply(self, Scalar):

+		return Vector3D(self._i * Scalar, self._j * Scalar, self._k * Scalar)

+	

+	def DivideByScalar(self, Scalar):

+		return self.ScalarMultiply(1. / Scalar)

+	

+	def Magnitude(self):

+		return (self._i ** 2. + self._j ** 2. + self._k ** 2.) ** .5

+	

+	def UnitVector(self):

+		return self / self.Magnitude()

+	

+	def AngleTo(self, vVector):

+		Theta = acos(self.Dot(vVector) / (self.Magnitude() * vVector.Magnitude()))

+		return Theta

+	

+	def AddVector(self, vVector):

+		return Vector3D(self._i + vVector._i, self._j + vVector._j, self._k + vVector._k)

+	

+	def SubtractVector(self, vVector):

+		return self.AddVector(-vVector)

+	

+	def Negate(self):

+		return Vector3D(-self._i, -self._j, -self._k)

+	

+	def ConvertToFrame(self, iVector1, jVector1, kVector1):

+		if abs(iVector1.Magnitude() - 1.) > .0001:

+			iVector1 = iVector1.UnitVector()

+		if abs(jVector1.Magnitude() - 1.) > .0001:

+			jVector1 = jVector1.UnitVector()

+		if abs(kVector1.Magnitude() - 1.) > .0001:

+			kVector1 = kVector1.UnitVector()

+		

+		iPart = self.Dot(iVector1)

+		jPart = self.Dot(jVector1)

+		kPart = self.Dot(kVector1)

+		

+		return Vector3D(iPart, jPart, kPart)

+	

+	def __neg__(self):

+		return self.Negate()

+	

+	def __mul__(self, vVector):

+		if isinstance(vVector, Vector3D):

+			return self.LeftCross(vVector)

+		if isinstance(vVector, Number):

+			return self.ScalarMultiply(vVector)

+		

+		return NotImplemented

+	

+	def __rmul__(self, uVector):

+		if isinstance(uVector, Vector3D):

+			return self.RightCross(uVector)

+		if isinstance(uVector, Number):

+			return self.ScalarMultiply(uVector)

+		

+		return NotImplemented

+	

+	def __add__(self, vVector):

+		if isinstance(vVector, Vector3D):

+			return self.AddVector(vVector)

+		

+		return NotImplemented

+	

+	def __sub__(self, vVector):

+		if isinstance(vVector, Vector3D):

+			return self.SubtractVector(vVector)

+		

+		return NotImplemented

+	

+	def __truediv__(self, Scalar):

+		if isinstance(Scalar, Number):

+			return self.DivideByScalar(Scalar)

+		

+		return NotImplemented

+	

+	def __eq__(self, vVector):

+		if isinstance(vVector, Vector3D):

+			return (self._i == vVector._i) and (self._j == vVector._j) and (self._k == vVector._k)

+		

+		return NotImplemented

+	

+	def __neq__(self, vVector):

+		return not self.__eq__(vVector) 

+	

+	def __format__(self, FormatSpec):

+		return("[{0}i, {1}j, {2}k]".format(self._i.__format__(FormatSpec), self._j.__format__(FormatSpec), self._k.__format__(FormatSpec)))

+	

+	def __str__(self):

+		return self.__format__("")

+	

+	def __repr__(self):

+		return("Vector3D({0}, {1}, {2})".format(self._i, self._j, self._k))

+

+iVector = Vector3D(1, 0, 0)

+jVector = Vector3D(0, 1, 0)

+kVector = Vector3D(0, 0, 1)

+

+if __name__ == "__main__":

+	from random import random

+	

+	uVector = Vector3D(random() * 10. - 5., random() * 10. - 5., random() * 10. - 5.)

+	print("u: {0}".format(uVector))

+	vVector = Vector3D(random() * 10. - 5., random() * 10. - 5., random() * 10. - 5.)

+	print("v: {0}".format(vVector))

+	

+	print("u = repr(u) = {0}".format(repr(uVector)))

+	

+	print("u × v = {0}".format(uVector * vVector))

+	print("v × u = {0}".format(vVector * uVector))

+	

+	print("u · v = {0}".format(uVector.Dot(vVector)))

+	print("v · u = {0}".format(vVector.Dot(uVector)))

+	

+	print("|u| = {0}".format(uVector.Magnitude()))

+	print("|v| = {0}".format(vVector.Magnitude()))

+	

+	print("û = {0}".format(uVector.UnitVector()))

+	print("|û| = {0}".format(uVector.UnitVector().Magnitude()))

+	

+	print("Angle between u and v: {0}".format(uVector.AngleTo(vVector)))

+	print("Angle between v and u: {0}".format(vVector.AngleTo(uVector)))

+	

+	print("u + v = {0}".format(uVector + vVector))

+	print("v + u = {0}".format(vVector + uVector))

+	

+	print("-u = {0}".format(-uVector))

+	print("u - v = {0}".format(uVector - vVector))

+	print("v - u = {0}".format(vVector - uVector))

+	

+	a = random() * 10. - 5.

+	print("a = {0}".format(a))

+	

+	print("u / a = {0}".format(uVector / a))