Orbit: So many changes. Fixed the update and output code in both
Orbit: So many changes. Fixed the update and output code in both
directions. Broke out some of the update functionality into useful
"get" functions. Added some other utility functions.
Body: Added an "Orbit" parameter.
KerbalBodies: Began seeding KSP data.

file:a/Body.py -> file:b/Body.py
--- a/Body.py
+++ b/Body.py
@@ -5,10 +5,11 @@
 

 class Body():

 	# Mass in kg, Radius in m.

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

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

 		self.Name = Name

 		self.Mass = Mass

 		self.Radius = Radius

+		self.Orbit = Orbit

 	

 	# Gravitational parameter µ in m³/s²

 	def GravParm(self):


file:b/KerbalBodies.py (new)
--- /dev/null
+++ b/KerbalBodies.py
@@ -1,1 +1,13 @@
+#!/usr/bin/python3
 
+from Body import Body
+from Orbit import Orbit
+from math import pi
+
+Kerbol = Body("Kerbol", 1.7565670e28, 261.6e6)
+Moho = Body("Moho", 2.5263617e21, 250e3, Orbit(Kerbol, 5263138304, 0.2, 7. * pi / 180., 70. * pi / 180., 15. * pi / 180., pi))
+Eve = Body("Eve", 1.2244127e23, 700e3, Orbit(Kerbol, 9832684544, 0.01, 2.1 * pi / 180., 15. * pi / 180., 0, pi))
+Gilly = Body("Gilly", 1.2420512e17, 13e3, Orbit(Eve, 31.5e6, 0.55, 12. * pi / 180., 80. * pi / 180., 10. * pi / 180., 0.9))
+Kerbin = Body("Kerbin", 5.2915793e22, 600e3, Orbit(Kerbol, 13599840256, 0, 0, 0, 0, pi))
+Mun = Body("Mün", 9.7600236e20, 200e3, Orbit(Kerbin, 12e6, 0, 0, 0, 0, 0.9))
+Minmus = Body("Minmus", 2.6457897e19, 60e3, Orbit(Kerbin, 47e6, 0, 6. * pi / 180., 78. * pi / 180., 38. * pi / 180., 1.7))

file:a/Orbit.py -> file:b/Orbit.py
--- a/Orbit.py
+++ b/Orbit.py
@@ -24,6 +24,7 @@
 		@type w: Number

 		@param M: Mean Anomaly [radians]

 		'''

+		

 		self.Primary = Primary

 		self.a = a

 		self.e = e

@@ -31,42 +32,120 @@
 		self.Omega = Omega

 		self.w = w

 		self.M = M

+		self._EIsDirty = True

 	

 	def EccentricAnomaly(self, M = None):

+		'''

+		Returns the Eccentric Anomaly in radians, optionally at a specified Mean

+		Anomaly.

+		

+		@param M: Mean Anomaly [radians] (Optional)

+		@type M: Number

+		

+		TODO: Add hyperbolic code for e > 1.

+		'''

+		

+		# Check if we have a clean cache to avoid the loop.

+		if not self._EIsDirty and M is None:

+			return self._E

+		

 		if M is None:

 			M = self.M

 		

+		# If the orbit is highly eccentric, start the guess with E=pi.

 		if self.e > 0.8:

 			E = pi

+		# Otherwise, start the guess at M.

 		else:

 			E = M

-		

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

+		i = 0

+		

+		# Loop through a Newton solver to find E from Kepler's equation.

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

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

-		

+			i += 1

+		

+		# Record E and clean the cache.

+		self._E = E

+		self._EIsDirty = False

 		return E

 	

 	def TrueAnomaly(self, E = None):

+		'''

+		Returns the True Anomaly in radians, optionally at a specified Eccentric

+		Anomaly.

+		

+		@param E: Eccentric Anomaly [radians] (Optional)

+		@type E: Number

+		'''

 		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 SemiParameter(self):

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

+	def MeanMotion(self):

+		'''

+		Returns the Mean Motion in radians / second.

+		'''

+		return (self.Primary.GravParm() / self.a ** 3.) ** .5

+	

+	def PeriapsisVector(self):

+		'''

+		Returns the ijk vector of the periapsis point.

+		@return: Periapsis vector [m, m, m]

+		@rtype: Vector3D

+		'''

+		return self.PosVector(0)

+	

+	def Periapsis(self):

+		'''

+		Returns the radius at the apoapsis.

+		@return: Periapsis radius [m]

+		@rtype: Number

+		'''

+		return self.a * (1. - self.e)

+	

+	def ApoapsisVector(self):

+		'''

+		Returns the ijk vector of the apoapsis point.

+		@return: Apoapsis vector [m, m, m]

+		@rtype: Vector3D

+		'''

+		return self.PosVector(pi)

+	

+	def Apoapsis(self):

+		'''

+		Returns the radius at the apoapsis.

+		@return: Apoapsis radius [m]

+		@rtype: Number

+		'''

+		return self.a * (1. + self.e)

 	

 	def _PQVectors(self):

-		Px = cos(self.Omega) * cos(self.w) - sin(self.Omega) * sin(self.w) * cos(self.i)

-		Py = cos(self.Omega) * sin(self.w) + sin(self.Omega) * cos(self.w) * cos(self.i)

-		Pz = sin(self.Omega) * sin(self.i)

-		

-		Qx = -sin(self.Omega) * cos(self.w) - cos(self.Omega) * sin(self.w) * cos(self.i)

-		Qy = -sin(self.Omega) * sin(self.w) + cos(self.Omega) * cos(self.w) * cos(self.i)

-		Qz = sin(self.i) * cos(self.Omega)

+		'''

+		Returns "P" and "Q" vectors for rotating the orbital frame about the

+		primary when determining the cartesian position of the orbiting body.

+		

+		@rtype: (Vector3D, Vector3D)

+		'''

+		Px = cos(self.w) * cos(self.Omega) - sin(self.w) * sin(self.Omega) * cos(self.i)

+		Py = cos(self.w) * sin(self.Omega) + sin(self.w) * cos(self.Omega) * cos(self.i)

+		Pz = sin(self.w) * sin(self.i)

+		

+		Qx = -sin(self.w) * cos(self.Omega) - cos(self.w) * sin(self.Omega) * cos(self.i)

+		Qy = -sin(self.w) * sin(self.Omega) + cos(self.w) * cos(self.Omega) * cos(self.i)

+		Qz = sin(self.i) * cos(self.w)

 		

 		return (Vector3D(Px, Py, Pz), Vector3D(Qx, Qy, Qz))

 	

 	def PosVector(self, M = None):

+		'''

+		Returns the position state vector in meters, optionally at specified

+		Mean Anomaly.

+		

+		@param M: Mean Anomaly [radians] (Optional)

+		@type M: Number 

+		'''

 		a = self.a

 		e = self.e

 		E = self.EccentricAnomaly(M)

@@ -77,35 +156,69 @@
 		return R

 	

 	def VelVector(self, M = None):

+		'''

+		Returns the velocity state vector in meters, optionally at specified

+		Mean Anomaly.

+		

+		@param M: Mean Anomaly [radians] (Optional)

+		@type M: Number

+		'''

 		a = self.a

 		e = self.e

 		E = self.EccentricAnomaly(M)

 		P, Q = self._PQVectors()

-		EDot = ((self.Primary.GravParm() / self.a ** 3.) ** .5) / (1 - self.e * cos(E))

+		EDot = self.MeanMotion() / (1 - self.e * cos(E))

 		

 		V = -a * sin(E) * EDot * P + a * ((1 - e ** 2.) ** .5) * cos(E) * EDot * Q

 		

 		return V

 	

-	def UpdateFromVectors(self, PosVector, VelVector):

-		'''

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

+	def GetStateFromMeanAnomaly(self, M):

+		'''

+		Returns a tuple containing the R and V state vectors given mean anomaly "M". 

+		

+		@param M: Mean anomaly [radians]

+		@type M: Number

+		'''

+		return (self.PosVector(M), self.VelVector(M))

+	

+	def FindAdditionalStateVectors(self, PosVector, VelVector):

+		'''

+		Calculates and returns the specific angular momentum, eccentricity, and

+		ascending node vectors given the 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

 		'''

 		

+		HVector = PosVector * VelVector

+		

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

+		

+		ANVector = kVector * HVector

+		

+		return (HVector, eVector, ANVector)

+	

+	def GetElementsFromStateVectors(self, PosVector, VelVector):

+		'''

+		Returns 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

+		'''

+		

+		HVector, eVector, ANVector = self.FindAdditionalStateVectors(PosVector, VelVector)

+		

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

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

-		HVector = PosVector * VelVector

-		

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

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

+		

 		e = eVector.Magnitude()

 		

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

-		

-		ANVector = kVector * HVector

 		

 		try:

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

@@ -114,16 +227,10 @@
 		except ZeroDivisionError:

 			Omega =  0.0

 		

-		if i != 0:

-			Omega = atan2(HVector._j, -HVector._i)

-		else:

-			Omega = 0

-		

 		try:

 			n = ANVector.UnitVector()

-			w = asin(n.Dot(eVector) / (n.Magnitude() * e))

-			

-			if eVector._j < -1e-6:

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

+			if eVector._k < -1e-6:

 				w = 2. * pi - w

 		except ZeroDivisionError:

 			try:

@@ -146,49 +253,227 @@
 			if PosVector.Dot(VelVector) < 0:

 				f = 2. * pi - f

 		

+		return (a, e, i, Omega, w, f)

+	

+	def GetEccentricAnomalyFromTrueAnomaly(self, f):

+		'''

+		Returns the eccentric anomaly at a given true anomaly.

+		

+		@param f: True Anomaly [radians]

+		@type f: Number

+		

+		@return: Eccentric Anomaly [radians]

+		@rtype: Number

+		'''

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

+		if E < 0:

+			E = 2 * pi + E

+		

+		return E

+	

+	def GetMeanAnomalyFromEccentricAnomaly(self, E):

+		'''

+		Returns the mean anomaly at a given eccentric anomaly.

+		

+		@param E: Eccentric anomaly [radians]

+		@type E: Number

+		

+		@return: Mean anomaly [radians]

+		@rtype: Number

+		'''

+		return (E - self.e * sin(E))

+	

+	def GetMeanAnomalyFromTrueAnomaly(self, f):

+		'''

+		Returns the mean anomaly at a given true anomaly.

+		@param f: True anomaly [radians]

+		@type f: Number

+		

+		@return: Mean anomaly [radians]

+		@rtype: Number

+		'''

+		return self.GetMeanAnomalyFromEccentricAnomaly(self.GetEccentricAnomalyFromTrueAnomaly(f))

+	

+	def GetTrueAnomalyFromPositionVector(self, PosVector, Tol = 1e-3):

+		'''

+		Attempts to find and return the true anomaly at a given position vector.

+		If it fails, it prints an error and returns None.

+		

+		@param PosVector: The position vector at which to find the true anomaly

+		@type PosVector: Number

+		@param Tol: The tolerance to which the positions should match, in meters.  Defaults to .001.

+		@type Tol: Number

+		

+		@return: True Anomaly

+		@type: Number

+		'''

+		e = self.e

+		i = self.i

+		

+		PeriapsisPosVector, PeriapsisVelVector = self.GetStateFromMeanAnomaly(0.)

+		_, eVector, ANVector = self.FindAdditionalStateVectors(PeriapsisPosVector, PeriapsisVelVector)

+		

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

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

+		#	if VelVector._i > 1e-6:

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

+		

+		CheckPos1, _ = self.GetStateFromMeanAnomaly(self.GetMeanAnomalyFromTrueAnomaly(f))

+		CheckPos2, _ = self.GetStateFromMeanAnomaly(self.GetMeanAnomalyFromTrueAnomaly(2. * pi - f))

+		

+		if (CheckPos1 - PosVector).Magnitude() < Tol:

+			return f

+		if (CheckPos2 - PosVector).Magnitude() < Tol:

+			return (2. * pi - f)

+		

+		print("GetTrueAnomalyFromPositionVector: Position not on orbit")

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

+		print("CheckPos1: {0}".format(CheckPos1), (CheckPos1 - PosVector).Magnitude())

+		print("CheckPos2: {0}".format(CheckPos2), (CheckPos2 - PosVector).Magnitude())

+	

+	def UpdateFromTrueAnomaly(self, f):

+		'''

+		Updates the mean anomaly from a given true anomaly.

+		@param f: True Anomaly [radians]

+		@type f: Number

+		'''

+		self.M = self.GetMeanAnomalyFromTrueAnomaly(f)

+	

+	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

+		'''

+		

+		a, e, i, Omega, w, f = self.GetElementsFromStateVectors(PosVector, VelVector)

+		

 		self.a = a

 		self.e = e

 		self.i = i

 		self.Omega = Omega

 		self.w = w

-		self.UpdateTrueAnomaly(f)

-		

-	

-	def UpdateTrueAnomaly(self, 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):

+		self.UpdateFromTrueAnomaly(f)

+	

+	def ConvertVectortoPRNFrame(self, Vector, M = None):

+		'''

+		Returns a Vector3D converted from the cartesian i, j, k frame to the

+		prograde, radial, normal frame described by the orbital state,

+		optionally at the specified mean anomaly.

+		 

+		@param Vector: The x, y, z vector to conver to P, R, N.

+		@type Vector: Vector3D

+		@param M: Mean anomaly [radians] (optional)

+		@type M: Number

+		'''

+		PosVector = self.PosVector(M)

+		VelVector = self.VelVector(M)

+		

+		PVector = self.VelVector(M).UnitVector()

+		NVector = (PosVector * VelVector).UnitVector()

+		RVector = PVector * NVector

+		

+		return Vector.ConvertToFrame(PVector, RVector, NVector)

+	

+	def __setattr__(self, Name, Value):

+		'''

+		Overrides the default attribute setter.

+		@param Name: The name of the class attribute being set.

+		@param Value: The value to which the named class attribute will be set.

+		'''

+		

+		# When M is set, dirty and wipe the cached eccentric anomaly.

+		if Name == 'M':

+			self.__dict__['_EIsDirty'] = True

+			self.__dict__['_E'] = None

+			return super().__setattr__(Name, Value)

+		# When e is set, check if the orbit is circular or elliptical.  If not,

+		# report NYI.  Also, make sure the eccentricity isn't negative.

+		if Name == 'e':

+			if Value is None:

+				return super().__setattr__(Name, Value)

+			if Value >= 1:

+				raise NotImplementedError("Parabolic and hyperbolic trajectories are not yet implemented.")

+			if Value < 0:

+				raise ValueError("The eccentricity must be in the domain [0, 1).")

+			return super().__setattr__(Name, Value)

+		

+		

+		return super().__setattr__(Name, Value)

+	

+	def GetString(self, M = None):

+		if M is None:

+			M = self.M

+		'''

+		Returns a long, multi-line string representing the orbit, optionally at

+		specified mean anomaly.

+		@param M: Mean anomaly [radians]

+		@type M: Number 

+		'''

 		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 += "M = {0} [radians]\n".format(M)

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

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

 		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(M), self.PosVector(M).Magnitude())

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

 		

 		return s

 	

+	def __str__(self):

+		'''

+		Returns a long, multi-line string representing the orbit.

+		'''

+		

+		return self.GetString()

+	

 if __name__ == "__main__":

+	from random import random

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

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

-	Orbit1.UpdateTrueAnomaly(pi / 2.)

+	Orbit1 = Orbit(Kerbin, random() * 3e6 + 7e5, random(), random() * pi / 2, random() * 2 * pi, 3 * pi / 2 + random() * pi / 2, random() * 2 * pi)

 	print(Orbit1)

 	

 	Orbit2 = Orbit(Kerbin)

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

-	print()

 	print(Orbit2)

 	

-	print(Orbit1.PosVector() - Orbit2.PosVector())

-		
+	PosError = Orbit1.PosVector() - Orbit2.PosVector()

+	VelError = Orbit1.VelVector() - Orbit2.VelVector()

+	

+	if PosError.Magnitude() > 1e-3:

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

+	if VelError.Magnitude() > 1e-3:

+		print("PosError: {0}".format(VelError))

+	

+	Orbit3 = Orbit(Kerbin, 1803823.21032378, 0.21141372855476, 55. * pi / 180., 0, 0, pi)

+	print(Orbit3)

+	

+	Orbit4 = Orbit(Kerbin, 2185176.20087195, 0, 55. * pi / 180., 0, 0, 0)

+	

+	Orbit4AtOrbit3TrueAnomaly = Orbit4.GetTrueAnomalyFromPositionVector(Orbit3.PosVector())

+	print("Current position on Orbit3 lies on Orbit4 at true anomaly {0}.".format(Orbit4AtOrbit3TrueAnomaly))

+	Orbit4.UpdateFromTrueAnomaly(Orbit4AtOrbit3TrueAnomaly)

+	

+	print(Orbit4)

+	

+	print("Delta-V from Orbit3 to Orbit4 at {0} radians: {1:.6}".format(Orbit4AtOrbit3TrueAnomaly, Orbit4.VelVector() - Orbit3.VelVector()))

+	PRNDeltaV = Orbit4.ConvertVectortoPRNFrame(Orbit4.VelVector() - Orbit3.VelVector())

+	PRNDeltaV = round(PRNDeltaV, 6)

+	print("Delta-V in PRN: {0:.6}".format(PRNDeltaV))

+

--- a/Vector3D.py
+++ b/Vector3D.py
@@ -73,6 +73,9 @@
 		kPart = self.Dot(kVector1)

 		

 		return Vector3D(iPart, jPart, kPart)

+	

+	def GetTuple(self):

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

 	

 	def __neg__(self):

 		return self.Negate()