Pythonized everything so there are no more PEP warnings. master
Pythonized everything so there are no more PEP warnings.

file:a/Body.py -> file:b/Body.py
--- a/Body.py
+++ b/Body.py
@@ -1,16 +1,34 @@
 #!/usr/bin/python3

+

+__all__ = ["G", "Body"]

+

+from numbers import Real

+

 

 # 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, Orbit = None):

-		self.Name = Name

-		self.Mass = Mass

-		self.Radius = Radius

-		self.Orbit = Orbit

-	

-	# Gravitational parameter µ in m³/s²

-	def GravParm(self):

-		return self.Mass * G
+    """

+    :type Name: str

+    :type Mass: Real

+    :type Radius: Real

+    :type Orbit: Orbit

+    :type mu: Real

+    """

+    # Mass in kg, Radius in m.

+    def __init__(self, name, mass=None, radius=None, orbit=None):

+        self.Name = name

+        self.Mass = mass

+        self.Radius = radius

+        self.Orbit = orbit

+

+    def gravitational_parameter(self):

+        """

+        Returns the gravitational parameter µ = M×G in m³/s²

+        :rtype: Real

+        """

+        return self.Mass * G

+

+    mu = property(gravitational_parameter)

--- a/KerbalBodies.py
+++ b/KerbalBodies.py
@@ -1,24 +1,75 @@
 #!/usr/bin/python3
+
+from math import pi
 
 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))
-Duna = Body("Duna", 4.5154812e21, 320e3, Orbit(Kerbol, 20726155264, 0.05, 0.06 * pi / 180., 135.5 * pi / 180., 0, pi))
-Ike = Body("Ike", 2.7821949e20, 130e3, Orbit(Duna, 3.2e6, 0.03, 0.2 * pi / 180., 0, 0, 1.7))
-Dres = Body("Dres", 3.2191322e20, 138e3, Orbit(Kerbol, 40839348203, 0.14, 5. * pi / 180., 280. * pi / 180., pi / 2., pi))
-Jool = Body("Jool", 4.2332635e24, 6e6, Orbit(Kerbol, 68773560320, 0.05, 1.304 * pi / 180., 52. * pi / 180., 0, 0.1))
-Laythe = Body("Laythe", 2.9397663e22, 500e3, Orbit(Jool, 27.184e6, 0, 0, 0, 0, pi))
-Vall = Body("Vall", 3.1088028e21, 300e3, Orbit(Jool, 43.152e6, 0, 0, 0, 0, 0.9))
-Tylo = Body("Tylo", 4.2332635e22, 600e3, Orbit(Jool, 68.5e6, 0, 0.025 * pi / 180., 0, 0, pi))
-Bop = Body("Bop", 3.7261536e19, 65e3, Orbit(Jool, 128.5e6, 0.24, 15. * pi / 180., 10. * pi / 180., 25. * pi / 180., 0.9))
-Pol = Body("Pol", 1.0813636e19, 44e3, Orbit(Jool, 179.89e6, 0.17, 4.25 * pi / 180., 2. * pi / 180., 15. * pi / 180., 0.9))
-Eeloo = Body("Eeloo", 1.1149358e21, 210e3, Orbit(Kerbol, 90118820000, 0.26, 6.15 * pi / 180., 50. * pi / 180., 260. * pi / 180., pi))
 
+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)
+              )
+
+Duna = Body("Duna", 4.5154812e21, 320e3,
+            Orbit(Kerbol, 20726155264, 0.05, 0.06 * pi / 180., 135.5 * pi / 180., 0, pi)
+            )
+
+Ike = Body("Ike", 2.7821949e20, 130e3,
+           Orbit(Duna, 3.2e6, 0.03, 0.2 * pi / 180., 0, 0, 1.7)
+           )
+
+Dres = Body("Dres", 3.2191322e20, 138e3,
+            Orbit(Kerbol, 40839348203, 0.14, 5. * pi / 180., 280. * pi / 180., pi / 2., pi)
+            )
+
+Jool = Body("Jool", 4.2332635e24, 6e6,
+            Orbit(Kerbol, 68773560320, 0.05, 1.304 * pi / 180., 52. * pi / 180., 0, 0.1)
+            )
+
+Laythe = Body("Laythe", 2.9397663e22, 500e3,
+              Orbit(Jool, 27.184e6, 0, 0, 0, 0, pi)
+              )
+
+Vall = Body("Vall", 3.1088028e21, 300e3,
+            Orbit(Jool, 43.152e6, 0, 0, 0, 0, 0.9)
+            )
+
+Tylo = Body("Tylo", 4.2332635e22, 600e3,
+            Orbit(Jool, 68.5e6, 0, 0.025 * pi / 180., 0, 0, pi)
+            )
+
+Bop = Body("Bop", 3.7261536e19, 65e3,
+           Orbit(Jool, 128.5e6, 0.24, 15. * pi / 180., 10. * pi / 180., 25. * pi / 180., 0.9)
+           )
+
+Pol = Body("Pol", 1.0813636e19, 44e3,
+           Orbit(Jool, 179.89e6, 0.17, 4.25 * pi / 180., 2. * pi / 180., 15. * pi / 180., 0.9)
+           )
+
+Eeloo = Body("Eeloo", 1.1149358e21, 210e3,
+             Orbit(Kerbol, 90118820000, 0.26, 6.15 * pi / 180., 50. * pi / 180., 260. * pi / 180., pi)
+             )
+
+

file:a/Orbit.py -> file:b/Orbit.py
--- a/Orbit.py
+++ b/Orbit.py
@@ -1,479 +1,539 @@
 #!/usr/bin/python3

 

+__all__ = ["Orbit"]

+

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

+from numbers import Real

+

 from Body import Body

-from Vector3D import *

-from math import acos, asin, atan2, cos, pi, sin

+from Vector3d import *

+

 

 class Orbit():

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

-		'''

-		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.Primary = Primary

-		self.a = a

-		self.e = e

-		self.i = i

-		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

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

-		'''

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

-		P, Q = self._PQVectors()

-		

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

-		

-		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.MeanMotion() / (1 - self.e * cos(E))

-		

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

-		

-		return V

-	

-	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 = -self.Primary.GravParm() / (2. * VvE)

-		

-		e = eVector.Magnitude()

-		

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

-		

-		try:

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

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

-				Omega =  2. * pi - Omega

-		except ZeroDivisionError:

-			Omega =  0.0

-		

-		try:

-			n = ANVector.UnitVector()

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

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

-				w = 2. * pi - w

-		except ZeroDivisionError:

-			try:

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

-			except ZeroDivisionError:

-				w = 0

-			if HVector._k < -1e-6:

-				w = 2. * pi - w

-		

-		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

-		

-		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.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 = VelVector.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(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(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()

-	

+    """

+    :type Primary: Body

+    :type a: Real

+    :type e: Real

+    :type i: Real

+    :type omega: Real

+    :type w: Real

+    :type mean_anomaly: Real

+    """

+    def __init__(self, primary, a=None, e=None, i=None, omega=None, w=None, mean_anomaly=None):

+        """

+        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: Real

+        @param e: Eccentricity [radians]

+        @type e: Real

+        @param i: Inclination [radians]

+        @type i: Real

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

+        @type omega: Real

+        @param w: Argument of Periapsis [radians]

+        @type w: Real

+        @param mean_anomaly: Mean Anomaly [radians]

+        """

+

+        self.Primary = primary

+        self.a = a

+        self._e = None

+        self.e = e

+        self.i = i

+        self.omega = omega

+        self.w = w

+        self.mean_anomaly = mean_anomaly

+        self._EIsDirty = True

+        self._E = None

+

+    def get_mean_anomaly_at_time(self, time):

+        return self.mean_anomaly + self.mean_motion() * time

+

+    def get_eccentric_anomaly_at_mean_anomaly(self, mean_anomaly=None):

+        """

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

+        Anomaly.

+

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

+        @type mean_anomaly: Real

+

+        TODO: Add hyperbolic code for e > 1.

+        """

+

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

+        if not self._EIsDirty and mean_anomaly is None:

+            return self._E

+

+        if mean_anomaly is None:

+            mean_anomaly = self.mean_anomaly

+

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

+        if self.e > 0.8:

+            eccentric_anomaly = pi

+        # Otherwise, start the guess at M.

+        else:

+            eccentric_anomaly = mean_anomaly

+        i = 0

+

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

+        while abs(eccentric_anomaly - self.e * sin(eccentric_anomaly) - mean_anomaly) > 1e-9:

+            eccentric_anomaly -= ((eccentric_anomaly - self.e * sin(eccentric_anomaly) - mean_anomaly) /

+                                  (1. - self.e * cos(eccentric_anomaly)))

+            i += 1

+

+        # Record E and clean the cache.

+        self._E = eccentric_anomaly

+        self._EIsDirty = False

+        return eccentric_anomaly

+

+    def get_true_anomaly_at_eccentric_anomaly(self, eccentric_anomaly=None):

+        """

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

+        Anomaly.

+

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

+        @type eccentric_anomaly: Real

+        """

+        if eccentric_anomaly is None:

+            eccentric_anomaly = self.get_eccentric_anomaly_at_mean_anomaly()

+

+        return 2. * atan2(sqrt(1. + self.e) * sin(eccentric_anomaly / 2.),

+                          sqrt(1. - self.e) * cos(eccentric_anomaly / 2.))

+

+    def mean_motion(self):

+        """

+        Returns the Mean Motion in radians / second.

+        """

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

+

+    def periapsis(self):

+        """

+        Returns the ijk vector of the periapsis point.

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

+        @rtype: Vector3d

+        """

+        return self.get_position_at_mean_anomaly(0)

+

+    def radius_at_periapsis(self):

+        """

+        Returns the radius at the apoapsis.

+        @return: Periapsis radius [m]

+        @rtype: Real

+        """

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

+

+    def apoapsis(self):

+        """

+        Returns the ijk vector of the apoapsis point.

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

+        @rtype: Vector3d

+        """

+        return self.get_position_at_mean_anomaly(pi)

+

+    def radius_at_apoapsis(self):

+        """

+        Returns the radius at the apoapsis.

+        @return: Apoapsis radius [m]

+        @rtype: Real

+        """

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

+

+    def _get_p_q_vectors(self):

+        """

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

+        """

+        p_x = cos(self.w) * cos(self.omega) - sin(self.w) * sin(self.omega) * cos(self.i)

+        p_y = cos(self.w) * sin(self.omega) + sin(self.w) * cos(self.omega) * cos(self.i)

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

+

+        q_x = -sin(self.w) * cos(self.omega) - cos(self.w) * sin(self.omega) * cos(self.i)

+        q_y = -sin(self.w) * sin(self.omega) + cos(self.w) * cos(self.omega) * cos(self.i)

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

+

+        return Vector3d(p_x, p_y, p_z), Vector3d(q_x, q_y, q_z)

+

+    def get_position_at_mean_anomaly(self, mean_anomaly=None):

+        """

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

+        Mean Anomaly.

+

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

+        @type mean_anomaly: Real

+        """

+

+        a = self.a

+        e = self.e

+        eccentric_anomaly = self.get_eccentric_anomaly_at_mean_anomaly(mean_anomaly)

+        p, q = self._get_p_q_vectors()

+

+        r = a * (cos(eccentric_anomaly) - e) * p + a * ((1 - e ** 2.) ** .5) * sin(eccentric_anomaly) * q

+

+        return r

+

+    def get_position_at_time(self, time=0):

+        """

+        Returns the position state vector in meters, optionally at the specified time since epoch.

+

+        @param time: Time since epoch [seconds] (Optional)

+        @type time: Real

+        @return Position: Vector3D

+        """

+        mean_anomaly = self.get_mean_anomaly_at_time(time)

+

+        return self.get_position_at_mean_anomaly(mean_anomaly)

+

+    def get_velocity_at_mean_anomaly(self, mean_anomaly=None):

+        """

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

+        Mean Anomaly.

+

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

+        @type mean_anomaly: Real

+        """

+        a = self.a

+        e = self.e

+        eccentric_anomaly = self.get_eccentric_anomaly_at_mean_anomaly(mean_anomaly)

+        p, q = self._get_p_q_vectors()

+        e_dot = self.mean_motion() / (1 - self.e * cos(eccentric_anomaly))

+

+        v = -a * sin(eccentric_anomaly) * e_dot * p + a * ((1 - e ** 2.) ** .5) * cos(eccentric_anomaly) * e_dot * q

+

+        return v

+

+    def get_state_at_mean_anomaly(self, mean_anomaly):

+        """

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

+

+        @param mean_anomaly: Mean anomaly [radians]

+        @type mean_anomaly: Real

+        """

+        return self.get_position_at_mean_anomaly(mean_anomaly), self.get_velocity_at_mean_anomaly(mean_anomaly)

+

+    def get_state_at_time(self, time):

+        """

+        Returns a tuple containing the R and V state vectors given time since epoch

+

+        @param time: time since epoch [s]

+        @type time: Real

+        """

+        return self.get_state_at_mean_anomaly(self.get_mean_anomaly_at_time(time))

+

+    def find_special_state_vectors(self, position, velocity):

+        """

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

+        ascending node vectors given the position and velocity vectors.

+

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

+        @type position: Vector3d

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

+        @type velocity: Vector3d

+        """

+

+        h_vector = position * velocity

+

+        e_vector = (velocity * h_vector) / self.Primary.gravitational_parameter() - position / position.magnitude()

+

+        an_vector = unit_k * h_vector

+

+        return h_vector, e_vector, an_vector

+

+    def get_elements_from_state(self, position, velocity):

+        """

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

+

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

+        @type position: Vector3d

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

+        @type velocity: Vector3d

+        """

+

+        h_vector, e_vector, an_vector = self.find_special_state_vectors(position, velocity)

+

+        vve = velocity.magnitude() ** 2. / 2. - self.Primary.gravitational_parameter() / position.magnitude()

+        a = -self.Primary.gravitational_parameter() / (2. * vve)

+

+        e = e_vector.magnitude()

+

+        i = acos(h_vector.k / h_vector.magnitude())

+

+        try:

+            omega = acos(an_vector.i / an_vector.magnitude())

+            if an_vector.j < -1e-6:

+                omega = 2. * pi - omega

+        except ZeroDivisionError:

+            omega = 0.0

+

+        try:

+            n = an_vector.unit_vector()

+            w = acos(n.dot(e_vector) / (n.magnitude() * e))

+            if e_vector.k < -1e-6:

+                w = 2. * pi - w

+        except ZeroDivisionError:

+            try:

+                w = atan2(e_vector.j / e, e_vector.i / e_vector.magnitude())

+            except ZeroDivisionError:

+                w = 0

+            if h_vector.k < -1e-6:

+                w = 2. * pi - w

+

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

+            f = acos(position.i / position.magnitude())

+            if velocity.i > 1e-6:

+                f = 2. * pi - f

+        elif e == 0:

+            f = acos(an_vector.dot(position) / (an_vector.magnitude() * position.magnitude()))

+            if an_vector.dot(velocity) > 0:

+                f = 2. * pi - f

+        else:

+            _arg = round(e_vector.dot(position) / (e_vector.magnitude() * position.magnitude()), 6)

+            f = acos(_arg)

+            if position.dot(velocity) < 0:

+                f = 2. * pi - f

+

+        return a, e, i, omega, w, f

+

+    def get_eccentric_anomaly_at_true_anomaly(self, f):

+        """

+        Returns the eccentric anomaly at a given true anomaly.

+

+        @param f: True Anomaly [radians]

+        @type f: Real

+

+        @return: Eccentric Anomaly [radians]

+        @rtype: Real

+        """

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

+        if eccentric_anomaly < 0:

+            eccentric_anomaly += 2 * pi

+

+        return eccentric_anomaly

+

+    def get_mean_anomaly_at_eccentric_anomaly(self, eccentric_anomaly):

+        """

+        Returns the mean anomaly at a given eccentric anomaly.

+

+        @param eccentric_anomaly: Eccentric anomaly [radians]

+        @type eccentric_anomaly: Real

+

+        @return: Mean anomaly [radians]

+        @rtype: Real

+        """

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

+

+    def get_mean_anomaly_at_true_anomaly(self, f):

+        """

+        Returns the mean anomaly at a given true anomaly.

+        @param f: True anomaly [radians]

+        @type f: Real

+

+        @return: Mean anomaly [radians]

+        @rtype: Real

+        """

+        return self.get_mean_anomaly_at_eccentric_anomaly(self.get_eccentric_anomaly_at_true_anomaly(f))

+

+    def get_true_anomaly_at_position(self, position, tolerance=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 position: The position vector at which to find the true anomaly

+        @type position: Vector3d

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

+        @type tolerance: Real

+

+        @return: True Anomaly

+        @type: Real

+        """

+        e = self.e

+        i = self.i

+

+        peri_pos, per_vel = self.get_state_at_mean_anomaly(0.)

+        _, e_vector, an_vector = self.find_special_state_vectors(peri_pos, per_vel)

+

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

+            f = acos(position.i / position.magnitude())

+        elif e == 0:

+            f = acos(an_vector.dot(position) / (an_vector.magnitude() * position.magnitude()))

+        else:

+            f = acos(e_vector.dot(position) / (e_vector.magnitude() * position.magnitude()))

+

+        check_pos1, _ = self.get_state_at_mean_anomaly(self.get_mean_anomaly_at_true_anomaly(f))

+        check_pos2, _ = self.get_state_at_mean_anomaly(self.get_mean_anomaly_at_true_anomaly(2. * pi - f))

+

+        if (check_pos1 - position).magnitude() < tolerance:

+            return f

+        if (check_pos2 - position).magnitude() < tolerance:

+            return 2. * pi - f

+

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

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

+        print("check_pos1: {0}".format(check_pos1), (check_pos1 - position).magnitude())

+        print("check_pos2: {0}".format(check_pos2), (check_pos2 - position).magnitude())

+

+    def update_orbit_from_true_anomaly(self, f):

+        """

+        Updates the mean anomaly from a given true anomaly.

+        @param f: True Anomaly [radians]

+        @type f: Real

+        """

+        self.mean_anomaly = self.get_mean_anomaly_at_true_anomaly(f)

+

+    def update_orbit_from_vectors(self, position, velocity):

+        """

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

+

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

+        @type position: Vector3d

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

+        @type velocity: Vector3d

+        """

+

+        a, e, i, omega, w, f = self.get_elements_from_state(position, velocity)

+

+        self.a = a

+        self.e = e

+        self.i = i

+        self.omega = omega

+        self.w = w

+        self.update_orbit_from_true_anomaly(f)

+

+    def vector_convert_ijk_to_prn_at_mean_anomaly(self, vector, mean_anomaly=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 mean_anomaly: Mean anomaly [radians] (optional)

+        @type mean_anomaly: Real

+        """

+        position = self.get_position_at_mean_anomaly(mean_anomaly)

+        velocity = self.get_velocity_at_mean_anomaly(mean_anomaly)

+

+        p = velocity.unit_vector()

+        n = (position * velocity).unit_vector()

+        r = -(p * n)

+

+        return vector.convert_reference_frame(p, r, n)

+

+    @classmethod

+    def vector_convert_prn_to_ijk(cls, vector):

+        i_part = vector.dot(unit_i)

+        j_part = vector.dot(unit_j)

+        k_part = vector.dot(unit_k)

+

+        return Vector3d(i_part, j_part, k_part)

+

+    def _get_mean_anomaly(self):

+        return self._mean_anomaly

+

+    def _set_mean_anomaly(self, value):

+        self._EIsDirty = True

+        self._E = None

+        self._mean_anomaly = value

+

+    mean_anomaly = property(_get_mean_anomaly, _set_mean_anomaly)

+

+    def _get_e(self):

+        return self._e

+

+    def _set_e(self, value):

+        if not isinstance(value, Real):

+            raise TypeError("Eccentricity must be a number type")

+        if value >= 1:

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

+        if value < 0:

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

+        self._e = value

+

+    e = property(_get_e, _set_e)

+

+    def to_string(self, mean_anomaly=None):

+        """

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

+        specified mean anomaly.

+        @param mean_anomaly: Mean anomaly [radians]

+        @type mean_anomaly: Real

+        """

+        if mean_anomaly is None:

+            mean_anomaly = self.mean_anomaly

+

+        s = []

+

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

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

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

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

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

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

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

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

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

+        s.append("\n")

+        s.append("Current Position: {0} (radius: {1})\n".format(

+            self.get_position_at_mean_anomaly(mean_anomaly),

+            self.get_position_at_mean_anomaly(mean_anomaly).magnitude()

+        ))

+        s.append("Current Velocity: {0} (magnitude: {1})\n".format(

+            self.get_velocity_at_mean_anomaly(mean_anomaly),

+            self.get_velocity_at_mean_anomaly(mean_anomaly).magnitude()

+        ))

+

+        return ''.join(s)

+

+    def __str__(self):

+        """

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

+        :rtype: str

+        """

+        return self.to_string()

+

 if __name__ == "__main__":

-	from random import random

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

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

-	

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

+    from random import random

+

+    Kerbin = Body("Kerbin", mass=5.2915793e22, radius=600e3)

+    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.update_orbit_from_vectors(Orbit1.get_position_at_mean_anomaly(), Orbit1.get_velocity_at_mean_anomaly())

+    print(Orbit2)

+

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

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

+

+    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.get_true_anomaly_at_position(Orbit3.get_position_at_mean_anomaly())

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

+    Orbit4.update_orbit_from_true_anomaly(Orbit4AtOrbit3TrueAnomaly)

+

+    print(Orbit4)

+

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

+                                                                        Orbit4.get_velocity_at_mean_anomaly() -

+                                                                        Orbit3.get_velocity_at_mean_anomaly()

+                                                                        ))

+    PRNDeltaV = Orbit4.vector_convert_ijk_to_prn_at_mean_anomaly(Orbit4.get_velocity_at_mean_anomaly() -

+                                                                 Orbit3.get_velocity_at_mean_anomaly()

+                                                                 )

+    PRNDeltaV = round(PRNDeltaV, 6)

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

 

--- a/Vector3D.py
+++ b/Vector3D.py
@@ -1,177 +1,297 @@
 #!/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):

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

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

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

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

-	

-	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 __round__(self, ndigits):

-		return Vector3D(round(self._i, ndigits), round(self._j, ndigits), 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)))

-	

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

+from numbers import Real

+from math import acos, sqrt

+

+__all__ = ['Vector3d', 'unit_i', 'unit_j', 'unit_k']

+

+

+class Vector3d:

+    """

+    :type i: Real

+    :type j: Real

+    :type k: Real

+    """

+    @classmethod

+    def cross_product(cls, left_vector, right_vector):

+        """

+        Return the cross product left_vector x right_vector.  Cross product is non-commutative.

+

+        :param left_vector: The LHS vector in the cross product

+        :type left_vector: Vector3d

+        :param right_vector: The RHS vector in the cross product

+        :type right_vector: Vector3d

+        :rtype: Vector3d

+        """

+        i = left_vector.j * right_vector.k - left_vector.k * right_vector.j

+        j = left_vector.k * right_vector.i - left_vector.i * right_vector.k

+        k = left_vector.i * right_vector.j - left_vector.j * right_vector.i

+

+        return cls(i, j, k)

+

+    @classmethod

+    def dot_product(cls, left_vector, right_vector):

+        """

+        Returns the dot product of two vectors.  Dot product is commutative.

+

+        :param left_vector: The first vector in the dot product

+        :type left_vector: Vector3d

+        :param right_vector: The second vector in the dot product

+        :type right_vector: Vector3d

+        :rtype: Real

+        """

+        return left_vector.i * right_vector.i + left_vector.j * right_vector.j + left_vector.k * right_vector.k

+

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

+        """

+        Initializes a new Vector3d object.

+

+        :param i: The "i" part of the vector

+        :type i: Real

+        :param j: The "j" part of the vector

+        :type j: Real

+        :param k: The "k" part of the vector

+        :type k: Real

+        """

+        self.i = i

+        self.j = j

+        self.k = k

+

+    def left_cross(self, right_vector):

+        """

+        Returns the cross of this vector with right_vector, self x right_vector

+

+        :param right_vector: The RHS vector in the cross product

+        :type right_vector: Vector3d

+        :rtype: Vector3d

+        """

+        return Vector3d.cross_product(self, right_vector)

+

+    def right_cross(self, left_vector):

+        """

+        Returns the cross of left_vector with this vector, left_vector x self

+

+        :param left_vector: The LHS vector in the cross product

+        :type left_vector: Vector3d

+        :rtype: Vector3d

+        """

+        return Vector3d.cross_product(left_vector, self)

+

+    def dot(self, other_vector):

+        """

+        Returns the dot product of this vector with another vector, self * other_vector

+

+        :param other_vector: The other vector in the dot product

+        :type other_vector: Vector3d

+        :rtype: Real

+        """

+        return Vector3d.dot_product(self, other_vector)

+

+    def scalar_multiply(self, scalar):

+        """

+        Returns a new vector equal to this vector multiplied by scalar.

+

+        :type scalar: Real

+        :rtype: Vector3d

+        """

+        return Vector3d(self.i * scalar, self.j * scalar, self.k * scalar)

+

+    def divide_by_scalar(self, scalar):

+        """

+        Returns a new vector equal to this vector divided by scalar

+

+        :type scalar: Real

+        :rtype: Vector3d

+        """

+        return self.scalar_multiply(1. / scalar)

+

+    def magnitude(self):

+        """

+        Returns the magnitude of this vector

+

+        :rtype: Real

+        """

+        return sqrt(self.sqr_magnitude())

+

+    def sqr_magnitude(self):

+        """

+        Returns the square of the magnitude of this vector

+

+        :rtype: Real

+        """

+        return self.i * self.i + self.j * self.j + self.k * self.k

+

+    def unit_vector(self):

+        """

+        Returns a new vector with magnitude 1 in the same direction as this vector

+

+        :rtype: Vector3d

+        """

+        return self / self.magnitude()

+

+    def angle_to(self, other_vector):

+        """

+        Returns the angle (in radians) from this vector to other_vector

+

+        :type other_vector: Vector3d

+        :rtype: Real

+        """

+        theta = acos(self.dot(other_vector) / (self.magnitude() * other_vector.magnitude()))

+        return theta

+

+    def add_vector(self, right_vector):

+        """

+        Returns a new vector equal to this vector plus right_vector

+

+        :param right_vector: The other vector being added to this vector

+        :rtype: Vector3d

+        """

+        return Vector3d(self.i + right_vector.i, self.j + right_vector.j, self.k + right_vector.k)

+

+    def negate(self):

+        """

+        Returns a new vector with equal but opposite components to this vector

+

+        :rtype: Vector3d

+        """

+        return Vector3d(-self.i, -self.j, -self.k)

+

+    def subtract_vector(self, right_vector):

+        """

+        Returns a new vector equal to this vector minus right_vector

+

+        :param right_vector: The RHS vector in the subtraction

+        :type right_vector: Vector3d

+        :rtype: Vector3d

+        """

+        return self.add_vector(-right_vector)

+

+    def convert_reference_frame(self, frame_i, frame_j, frame_k):

+        """

+        Returns a new vector in the arbitrary reference frame given by frame_i, frame_j, frame_k

+

+        :type frame_i: Vector3d

+        :type frame_j: Vector3d

+        :type frame_k: Vector3d

+        :rtype: Vector3d

+        """

+        if abs(frame_i.magnitude() - 1.) > .0001:

+            frame_i = frame_i.unit_vector()

+        if abs(frame_j.magnitude() - 1.) > .0001:

+            frame_j = frame_j.unit_vector()

+        if abs(frame_k.magnitude() - 1.) > .0001:

+            frame_k = frame_k.unit_vector()

+

+        i_part = self.dot(frame_i)

+        j_part = self.dot(frame_j)

+        k_part = self.dot(frame_k)

+

+        return Vector3d(i_part, j_part, k_part)

+

+    def get_as_tuple(self):

+        """

+        Converts this vector to a tuple (i, j, k)

+        :rtype: tuple[Real]

+        """

+        return self.i, self.j, self.k

+

+    def __neg__(self):

+        return self.negate()

+

+    def __mul__(self, right_vector):

+        if isinstance(right_vector, Vector3d):

+            return self.left_cross(right_vector)

+        if isinstance(right_vector, Real):

+            return self.scalar_multiply(right_vector)

+

+        return NotImplemented

+

+    def __rmul__(self, left_vector):

+        if isinstance(left_vector, Vector3d):

+            return self.right_cross(left_vector)

+        if isinstance(left_vector, Real):

+            return self.scalar_multiply(left_vector)

+

+        return NotImplemented

+

+    def __add__(self, right_vector):

+        if isinstance(right_vector, Vector3d):

+            return self.add_vector(right_vector)

+

+        return NotImplemented

+

+    def __sub__(self, right_vector):

+        if isinstance(right_vector, Vector3d):

+            return self.subtract_vector(right_vector)

+

+        return NotImplemented

+

+    def __truediv__(self, scalar):

+        if isinstance(scalar, Real):

+            return self.divide_by_scalar(scalar)

+

+        return NotImplemented

+

+    def __eq__(self, other_vector):

+        if isinstance(other_vector, Vector3d):

+            return (self.i == other_vector.i) and (self.j == other_vector.j) and (self.k == other_vector.k)

+

+        return NotImplemented

+

+    def __neq__(self, other_vector):

+        return not self.__eq__(other_vector)

+

+    def __round__(self, ndigits):

+        return Vector3d(round(self.i, ndigits), round(self.j, ndigits), round(self.j, ndigits))

+

+    def __format__(self, format_spec):

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

+                                            self.k.__format__(format_spec)))

+

+    def __str__(self):

+        return self.__format__("")

+

+    def __repr__(self):

+        return "vector3d({0}, {1}, {2})".format(self.i, self.j, self.k)

+

+

+unit_i = Vector3d(1, 0, 0)

+unit_j = Vector3d(0, 1, 0)

+unit_k = 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))
+    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.unit_vector()))

+    print("|û| = {0}".format(uVector.unit_vector().magnitude()))

+

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

+    print("Angle between v and u: {0}".format(vVector.angle_to(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))

file:b/Vector3d.py (new)
--- /dev/null
+++ b/Vector3d.py
@@ -1,1 +1,297 @@
-
+#!/usr/bin/python3
+
+from numbers import Real
+from math import acos, sqrt
+
+__all__ = ['Vector3d', 'unit_i', 'unit_j', 'unit_k']
+
+
+class Vector3d:
+    """
+    :type i: Real
+    :type j: Real
+    :type k: Real
+    """
+    @classmethod
+    def cross_product(cls, left_vector, right_vector):
+        """
+        Return the cross product left_vector x right_vector.  Cross product is non-commutative.
+
+        :param left_vector: The LHS vector in the cross product
+        :type left_vector: Vector3d
+        :param right_vector: The RHS vector in the cross product
+        :type right_vector: Vector3d
+        :rtype: Vector3d
+        """
+        i = left_vector.j * right_vector.k - left_vector.k * right_vector.j
+        j = left_vector.k * right_vector.i - left_vector.i * right_vector.k
+        k = left_vector.i * right_vector.j - left_vector.j * right_vector.i
+
+        return cls(i, j, k)
+
+    @classmethod
+    def dot_product(cls, left_vector, right_vector):
+        """
+        Returns the dot product of two vectors.  Dot product is commutative.
+
+        :param left_vector: The first vector in the dot product
+        :type left_vector: Vector3d
+        :param right_vector: The second vector in the dot product
+        :type right_vector: Vector3d
+        :rtype: Real
+        """
+        return left_vector.i * right_vector.i + left_vector.j * right_vector.j + left_vector.k * right_vector.k
+
+    def __init__(self, i=0, j=0, k=0):
+        """
+        Initializes a new Vector3d object.
+
+        :param i: The "i" part of the vector
+        :type i: Real
+        :param j: The "j" part of the vector
+        :type j: Real
+        :param k: The "k" part of the vector
+        :type k: Real
+        """
+        self.i = i
+        self.j = j
+        self.k = k
+
+    def left_cross(self, right_vector):
+        """
+        Returns the cross of this vector with right_vector, self x right_vector
+
+        :param right_vector: The RHS vector in the cross product
+        :type right_vector: Vector3d
+        :rtype: Vector3d
+        """
+        return Vector3d.cross_product(self, right_vector)
+
+    def right_cross(self, left_vector):
+        """
+        Returns the cross of left_vector with this vector, left_vector x self
+
+        :param left_vector: The LHS vector in the cross product
+        :type left_vector: Vector3d
+        :rtype: Vector3d
+        """
+        return Vector3d.cross_product(left_vector, self)
+
+    def dot(self, other_vector):
+        """
+        Returns the dot product of this vector with another vector, self * other_vector
+
+        :param other_vector: The other vector in the dot product
+        :type other_vector: Vector3d
+        :rtype: Real
+        """
+        return Vector3d.dot_product(self, other_vector)
+
+    def scalar_multiply(self, scalar):
+        """
+        Returns a new vector equal to this vector multiplied by scalar.
+
+        :type scalar: Real
+        :rtype: Vector3d
+        """
+        return Vector3d(self.i * scalar, self.j * scalar, self.k * scalar)
+
+    def divide_by_scalar(self, scalar):
+        """
+        Returns a new vector equal to this vector divided by scalar
+
+        :type scalar: Real
+        :rtype: Vector3d
+        """
+        return self.scalar_multiply(1. / scalar)
+
+    def magnitude(self):
+        """
+        Returns the magnitude of this vector
+
+        :rtype: Real
+        """
+        return sqrt(self.sqr_magnitude())
+
+    def sqr_magnitude(self):
+        """
+        Returns the square of the magnitude of this vector
+
+        :rtype: Real
+        """
+        return self.i * self.i + self.j * self.j + self.k * self.k
+
+    def unit_vector(self):
+        """
+        Returns a new vector with magnitude 1 in the same direction as this vector
+
+        :rtype: Vector3d
+        """
+        return self / self.magnitude()
+
+    def angle_to(self, other_vector):
+        """
+        Returns the angle (in radians) from this vector to other_vector
+
+        :type other_vector: Vector3d
+        :rtype: Real
+        """
+        theta = acos(self.dot(other_vector) / (self.magnitude() * other_vector.magnitude()))
+        return theta
+
+    def add_vector(self, right_vector):
+        """
+        Returns a new vector equal to this vector plus right_vector
+
+        :param right_vector: The other vector being added to this vector
+        :rtype: Vector3d
+        """
+        return Vector3d(self.i + right_vector.i, self.j + right_vector.j, self.k + right_vector.k)
+
+    def negate(self):
+        """
+        Returns a new vector with equal but opposite components to this vector
+
+        :rtype: Vector3d
+        """
+        return Vector3d(-self.i, -self.j, -self.k)
+
+    def subtract_vector(self, right_vector):
+        """
+        Returns a new vector equal to this vector minus right_vector
+
+        :param right_vector: The RHS vector in the subtraction
+        :type right_vector: Vector3d
+        :rtype: Vector3d
+        """
+        return self.add_vector(-right_vector)
+
+    def convert_reference_frame(self, frame_i, frame_j, frame_k):
+        """
+        Returns a new vector in the arbitrary reference frame given by frame_i, frame_j, frame_k
+
+        :type frame_i: Vector3d
+        :type frame_j: Vector3d
+        :type frame_k: Vector3d
+        :rtype: Vector3d
+        """
+        if abs(frame_i.magnitude() - 1.) > .0001:
+            frame_i = frame_i.unit_vector()
+        if abs(frame_j.magnitude() - 1.) > .0001:
+            frame_j = frame_j.unit_vector()
+        if abs(frame_k.magnitude() - 1.) > .0001:
+            frame_k = frame_k.unit_vector()
+
+        i_part = self.dot(frame_i)
+        j_part = self.dot(frame_j)
+        k_part = self.dot(frame_k)
+
+        return Vector3d(i_part, j_part, k_part)
+
+    def get_as_tuple(self):
+        """
+        Converts this vector to a tuple (i, j, k)
+        :rtype: tuple[Real]
+        """
+        return self.i, self.j, self.k
+
+    def __neg__(self):
+        return self.negate()
+
+    def __mul__(self, right_vector):
+        if isinstance(right_vector, Vector3d):
+            return self.left_cross(right_vector)
+        if isinstance(right_vector, Real):
+            return self.scalar_multiply(right_vector)
+
+        return NotImplemented
+
+    def __rmul__(self, left_vector):
+        if isinstance(left_vector, Vector3d):
+            return self.right_cross(left_vector)
+        if isinstance(left_vector, Real):
+            return self.scalar_multiply(left_vector)
+
+        return NotImplemented
+
+    def __add__(self, right_vector):
+        if isinstance(right_vector, Vector3d):
+            return self.add_vector(right_vector)
+
+        return NotImplemented
+
+    def __sub__(self, right_vector):
+        if isinstance(right_vector, Vector3d):
+            return self.subtract_vector(right_vector)
+
+        return NotImplemented
+
+    def __truediv__(self, scalar):
+        if isinstance(scalar, Real):
+            return self.divide_by_scalar(scalar)
+
+        return NotImplemented
+
+    def __eq__(self, other_vector):
+        if isinstance(other_vector, Vector3d):
+            return (self.i == other_vector.i) and (self.j == other_vector.j) and (self.k == other_vector.k)
+
+        return NotImplemented
+
+    def __neq__(self, other_vector):
+        return not self.__eq__(other_vector)
+
+    def __round__(self, ndigits):
+        return Vector3d(round(self.i, ndigits), round(self.j, ndigits), round(self.j, ndigits))
+
+    def __format__(self, format_spec):
+        return ("[{0}i, {1}j, {2}k]".format(self.i.__format__(format_spec), self.j.__format__(format_spec),
+                                            self.k.__format__(format_spec)))
+
+    def __str__(self):
+        return self.__format__("")
+
+    def __repr__(self):
+        return "vector3d({0}, {1}, {2})".format(self.i, self.j, self.k)
+
+
+unit_i = Vector3d(1, 0, 0)
+unit_j = Vector3d(0, 1, 0)
+unit_k = 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.unit_vector()))
+    print("|û| = {0}".format(uVector.unit_vector().magnitude()))
+
+    print("Angle between u and v: {0}".format(uVector.angle_to(vVector)))
+    print("Angle between v and u: {0}".format(vVector.angle_to(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))