Completely rebuilt Orbit. Its state vector calcs are still a bit wrong.
Completely rebuilt Orbit. Its state vector calcs are still a bit wrong.

file:a/Orbit.py -> file:b/Orbit.py
--- a/Orbit.py
+++ b/Orbit.py
@@ -5,168 +5,184 @@
 from math import acos, atan2, cos, pi, sin

 

 class Orbit():

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

+	def __init__(self, Primary, a = None, e = None, i = None, Omega = None, w = None, M = None):

 		'''

-		@param PosVector:

-		@type Vector3D:

-		@param VelVector:

-		@type Vector3D:

-		@param Primary:

-		@type Body:

+		Creates an Orbit object for determining the behavior of an object

+		with insignificant mass in orbit about a Primary body.

+		

+		@param Primary: The Primary body about which the object orbits.

+		@type Primary: Body

+		@param a: Semi-Major Axis [m]

+		@type a: Number

+		@param e: Eccentricity [radians]

+		@type e: Number

+		@param i: Inclination [radians]

+		@type i: Number

+		@param Omega: Longitude of the Ascending Node [radians]

+		@type Omega: Number

+		@param w: Argument of Periapsis [radians]

+		@type w: Number

+		@param M: Mean Anomaly [radians]

 		'''

-		self.PosVector = PosVector

-		self.VelVector = VelVector

 		self.Primary = Primary

+		self.a = a

+		self.e = e

+		self.i = i

+		self.Omega = Omega

+		self.w = w

+		self.M = M

 	

-	def VisVivaEnergy(self):

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

+	def EccentricAnomaly(self, M = None):

+		if M is None:

+			M = self.M

+		

+		if self.e > 0.8:

+			E = pi

+		else:

+			E = M

+		

+		while abs(E - self.e * sin(E) - M) > 1e-9:

+			E = E - (E - self.e * sin(E) - M) / (1. - self.e * cos(E))

+		

+		return E

 	

-	def SemiMajorAxis(self):

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

+	def TrueAnomaly(self, E = None):

+		if E is None:

+			E = self.EccentricAnomaly()

+		

+		return 2. * atan2((1. + self.e) ** .5 * sin(E / 2.), (1. - self.e) ** .5 * cos(E / 2.))

 	

-	def SpecAngMmntm(self):

-		return self.PosVector * self.VelVector

+	def SemiParameter(self):

+		return (self.a * (1. - self.e ** 2.))# / (1. + self.e * cos(self.TrueAnomaly()))

 	

-	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()

+	def PosVector(self, M = None):

+		i = self.i

+		Omega = self.Omega

+		w = self.w

+		f = self.TrueAnomaly(self.EccentricAnomaly(M))

+		p = self.SemiParameter() / (1. + self.e * cos(f))

 		

-		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)

+		

+		return Vector3D(Rx, Ry, Rz)

+	

+	def VelVector(self, M = None):

+		p = self.SemiParameter()

+		e = self.e

+		i = self.i

+		Omega = self.Omega

+		w = self.w

+		f = self.TrueAnomaly(self.EccentricAnomaly(M))

 		

 		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

+		return Vector3D(Vx, Vy, Vz)

 	

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

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

-		self.PosVector = PosVector

+	def UpdateFromVectors(self, PosVector, VelVector):

+		'''

+		Updates the orbital elements from a given set of position and velocity vectors.

+		@param PosVector: Position vector [m, m, m]

+		@type PosVector: Vector3D

+		@param VelVector: Velocity vector [m/s, m/s, m/s]

+		@type VelVector: Vector3D

+		'''

 		

-# 		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)

+		print("PosVector: {0}".format(PosVector))

+		print("VelVector: {0}".format(VelVector))

+		

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

+		a = -self.Primary.GravParm() / (2. * VvE)

+		HVector = PosVector * VelVector

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

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

+		print("r dot h = {0}".format(PosVector.Dot(HVector)))

+		eVector = (VelVector * HVector) / self.Primary.GravParm() - PosVector / PosVector.Magnitude()

+		e = eVector.Magnitude()

+		i = acos(HVector._k / HVector.Magnitude())

+		ANVector = kVector * HVector

+		

+		try:

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

+			if ANVector._j < 0:

+				Omega =  2. * pi - Omega

+		except ZeroDivisionError:

+			Omega =  0.0

+		

+		try:

+			n = ANVector.UnitVector()

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

+		

+			if eVector._k < 0:

+				w = 2. * pi - w

+		except ZeroDivisionError:

+			try:

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

+			except ZeroDivisionError:

+				w = 0

+			if HVector._k < 0:

+				w = 2. * pi - w

+		if e == 0 and i == 0:

+			f = acos(PosVector._i / PosVector.Magnitude())

+			if VelVector._i > 0:

+				f = 2. * pi - f

+		elif e == 0:

+			f = acos(ANVector.Dot(PosVector) / (ANVector.Magnitude() * PosVector.Magnitude()))

+			if ANVector.Dot(VelVector) > 0:

+				f = 2. * pi - f

+		else:

+			f = acos(eVector.Dot(PosVector) / (eVector.Magnitude() * PosVector.Magnitude()))

+			if PosVector.Dot(VelVector) < 0:

+				f = 2. * pi - f

+		

+		self.a = a

+		self.e = e

+		self.i = i

+		self.Omega = Omega

+		self.w = w

+		self.UpdateTrueAnomaly(f)

 	

 	def UpdateTrueAnomaly(self, f):

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

+		E = atan2((1 - self.e ** 2.) ** .5 * sin(f), self.e + cos(f))

+		if f > pi:

+			E = 2. * pi - E

+		

+		M = E - self.e * sin(E)

+		

+		self.M = M

 	

 	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  = "Orbiting '{0}'\n".format(self.Primary.Name)

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

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

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

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

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

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

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

 		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())

+		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.)

+	Orbit1 = Orbit(Kerbin, 4.35e6, 0.28, 33. * pi / 180., 0, 3. * pi / 2., pi / 2.)

+	# Orbit1.UpdateTrueAnomaly(pi / 2.)

 	print(Orbit1)

 	

-	Orbit2 = Orbit(Primary = Kerbin)

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

+	h = (Orbit1.a * (1 - Orbit1.e ** 2.) * Orbit1.Primary.GravParm()) ** .5

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

+	eps = -Orbit1.Primary.GravParm() / (2. * Orbit1.a)

+	eCheck = (1 + (2 * eps * h ** 2.) / (Orbit1.Primary.GravParm() ** 2.)) ** .5

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

+	

+	Orbit2 = Orbit(Kerbin)

+	Orbit2.UpdateFromVectors(Orbit1.PosVector(), Orbit1.VelVector())

 	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/Orbit_old.py (new)
--- /dev/null
+++ b/Orbit_old.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()))

--- a/Vector3D.py
+++ b/Vector3D.py
@@ -13,9 +13,13 @@
 	

 	@classmethod

 	def CrossProduct(cls, uVector, vVector):

+		print(uVector, vVector)

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

+		print(uVector._j * vVector._k, uVector._k * vVector._j, i)

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

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

+		print(uVector._k * vVector._i, uVector._i * vVector._k, j)

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

+		print(uVector._i * vVector._j, uVector._j * vVector._i, k)

 		return cls(i, j, k)

 	

 	def LeftCross(self, vVector):

@@ -116,6 +120,11 @@
 	def __neq__(self, vVector):

 		return not self.__eq__(vVector) 

 	

+	def __round__(self, ndigits):

+		self._i = round(self._i, ndigits)

+		self._j = round(self._j, ndigits)

+		self._k = round(self._j, ndigits)

+	

 	def __format__(self, FormatSpec):

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