import scipy as simport numpy as n q1 = a1 * (1.0 + e1) q2 = a2 * (1.0 + e2) I21 = s.arccos( s.cos(i1) * s.cos(i2) + s.sin(i1) * s.sin(i2) * s.cos(O2 - O1)) if n.absolute(O2 - O1) > 180: II21 = w2 - w1 - 2*s.arcsin( s.cos((i2 + i1)/2.0) * s.sin((O2 - O1)/2.0) * 1.0/s.cos(I21/2.0) ) else: II21 = w2 - w1 + 2*s.arcsin( s.cos((i2 + i1)/2.0) * s.sin((O2 - O1)/2.0) * 1.0/s.cos(I21/2.0) ) D2 = ((e2-e1) ** 2) + (((q2 - q1) / (q2 + q1)) ** 2) + ((2*s.sin(I21/2.0)) ** 2) + (((e2+e1)/2.0) ** 2) * ((2*s.sin(II21/2.0)) ** 2)