A whole bunch of possibly-wrong changes to Orbit.
A whole bunch of possibly-wrong changes to Orbit.

file:a/Orbit.py -> file:b/Orbit.py
--- a/Orbit.py
+++ b/Orbit.py
@@ -2,7 +2,7 @@
 

 from Body import Body

 from Vector3D import *

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

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

 

 class Orbit():

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

@@ -41,7 +41,7 @@
 		else:

 			E = M

 		

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

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

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

 		

 		return E

@@ -55,32 +55,37 @@
 	def SemiParameter(self):

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

 	

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

+		

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

+	

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

+		a = self.a

+		e = self.e

+		E = self.EccentricAnomaly(M)

+		P, Q = self._PQVectors()

 		

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

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

 		

-		return Vector3D(Rx, Ry, Rz)

+		return R

 	

 	def VelVector(self, M = None):

-		p = self.SemiParameter()

+		a = self.a

 		e = self.e

-		i = self.i

-		Omega = self.Omega

-		w = self.w

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

+		E = self.EccentricAnomaly(M)

+		P, Q = self._PQVectors()

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

 		

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

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

 		

-		return Vector3D(Vx, Vy, Vz)

+		return V

 	

 	def UpdateFromVectors(self, PosVector, VelVector):

 		'''

@@ -91,43 +96,46 @@
 		@type VelVector: Vector3D

 		'''

 		

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

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

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

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

+		HVector = PosVector * 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:

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

 				Omega =  2. * pi - Omega

 		except ZeroDivisionError:

 			Omega =  0.0

 		

+		if i != 0:

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

+		else:

+			Omega = 0

+		

 		try:

 			n = ANVector.UnitVector()

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

-		

-			if eVector._k < 0:

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

+			

+			if eVector._j < -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 < 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 > 0:

+			if VelVector._i > 1e-6:

 				f = 2. * pi - f

 		elif e == 0:

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

@@ -144,6 +152,7 @@
 		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))

@@ -172,17 +181,14 @@
 	

 if __name__ == "__main__":

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

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

-	# Orbit1.UpdateTrueAnomaly(pi / 2.)

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

+	Orbit1.UpdateTrueAnomaly(pi / 2.)

 	print(Orbit1)

-	

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

 	print(Orbit2)

+	

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

 		

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

 	@classmethod

 	def CrossProduct(cls, uVector, vVector):

-		print(uVector, vVector)

+# 		print(uVector, vVector)

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

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

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

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

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

 		return cls(i, j, k)

 	

 	def LeftCross(self, vVector):

@@ -121,9 +121,7 @@
 		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)

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