Small correction to time period.

orbitdeterminator/propagation/cowell.py

@@ -1,7 +1,6 @@
 import numpy as np
-from orbitdeterminator.util.new_tle_kep_state import kep_to_state
 
-mu = 398600.4405
+mu = 398600.4418
 J2 = 1.08262668e-3
 Re = 6378.137
 we = 7.292115e-5
@@ -83,18 +82,42 @@ def rk4(s,t0,tf,h=30):
         s = s+(k1+2*k2+2*k3+k4)/6
         t = t+h
 
+ #       if (s[2]<0 and s[2]>-200 and s[5]>0):
+ #           dt = -s[2]/s[5]
+ #           print(t+dt)
+
     return s
 
-def propagate_state(s,t0,tf):
-    return rk4(s,t0,tf)
+def time_period(s,h=30):
+    t = 0
+
+    old_z, pass_1 = 0, None
+
+    while(True):
+        k1 = h*sdot(s)
+        k2 = h*sdot(s+k1/2)
+        k3 = h*sdot(s+k2/2)
+        k4 = h*sdot(s+k3)
+
+        s = s+(k1+2*k2+2*k3+k4)/6
+        t = t+h
 
-def propagate_kep(kep,t0,tf):
-    s = kep_to_state(kep)
+        if (s[2]>=0 and old_z<0):
+            dt = -s[2]/s[5]
+            t2 = t+dt
+
+            if pass_1 is None:
+                pass_1 = t2
+            else:
+                return t2-pass_1
+
+        old_z = s[2]
+
+def propagate_state(s,t0,tf):
     return rk4(s,t0,tf)
 
 if __name__ == "__main__":
-    #s = np.array([-1.25407719e+02,6.23615970e+03,2.67702463e+03,-5.21833653e+00,2.12226865e+00,-5.19435315e+00])
-    s = np.array([2.87327861e+03,5.22872234e+03,3.23884457e+03,-3.49536799e+00,4.87267295e+00,-4.76846910e+00])
+    s = np.array([2.87393871e+03,5.22992358e+03,3.23958865e+03,-3.49496655e+00,4.87211332e+00,-4.76792145e+00])
     t0, tf = 0, 88796.3088704
     final = propagate_state(s,t0,tf)
     print(final)

orbitdeterminator/propagation/dgsn_simulator.py

@@ -145,7 +145,7 @@ class SimParams():
     #epoch = 1529410874
     #iss_kep = np.array([6775,0.0002893,51.6382,211.1340,7.1114,148.9642])
 
-    iss_kep = np.array([6786.6420,0.0003456,51.6418,290.0933,266.6543,212.4306])
+    iss_kep = np.array([6785.6420,0.0003456,51.6418,290.0933,266.6543,212.4306])
     params = SimParams()
     params.kep = iss_kep
     params.epoch = epoch

orbitdeterminator/propagation/simulator.py

@@ -5,7 +5,7 @@
 from functools import partial
 import numpy as np
 
-from orbitdeterminator.propagation.cowell import propagate_state
+from orbitdeterminator.propagation.cowell import propagate_state,time_period
 from orbitdeterminator.util.teme_to_ecef import conv_to_ecef
 from orbitdeterminator.util.new_tle_kep_state import kep_to_state
 
@@ -105,6 +105,15 @@ def write(self,t,s):
     def close(self):
         self.f.close()
 
+class omega(OpWriter):
+    @staticmethod
+    def write(t,s):
+        #T = time_period(s)
+        r = np.linalg.norm(s[0:3])
+        v = np.linalg.norm(s[3:6])
+        T = 2*np.pi*r/v
+        print(t,T)
+
 class SimParams():
     kep = None
     epoch = None
@@ -115,12 +124,7 @@ class SimParams():
 
 if __name__ == "__main__":
     epoch = 1531152114
-    #epoch = 1530729961
-    #iss_kep = np.array([6775,0.0002893,51.6382,211.1340,7.1114,148.9642])
-    #iss_kep = np.array([6786.6787,0.0003411,51.6428,263.9950,291.0075,245.8091])
-    iss_kep = np.array([6786.6420,0.0003456,51.6418,290.0933,266.6543,212.4306])
-    #6783.1714
-    #6786.5714
+    iss_kep = np.array([6785.68682,0.0003456,51.6418,290.0933,266.6543,212.430557])
 
     params = SimParams()
     params.kep = iss_kep

orbitdeterminator/util/new_tle_kep_state.py

@@ -4,7 +4,9 @@
 import numpy as np
 from scipy.optimize import fsolve
 
-mu = 398600.4405
+from orbitdeterminator.propagation.cowell import time_period
+
+mu = 398600.4418
 
 def MtoE(M,e):
     """Calculates the eccentric anomaly from the mean anomaly.
@@ -40,19 +42,16 @@ def EtoT(E,e):
 def MtoT(M,e):
     return EtoT(MtoE(M,e),e)
 
-def ntoa(n,i):
-    T = 86400/n
-    i = math.radians(i)
-
-    Re = 6378.137
-    J2 = 1.08262668e-3
+def scale(s,T):
+    s = np.ravel(s)
+    tmp = np.empty(6)
 
-    B = 3*J2*Re**2*(8*math.sin(i)**2-6)/4
-    t = lambda a: T - 2*math.pi*(a**3/mu)**0.5*(1+B/a**2)
-    a0 = (mu*(T/2/math.pi)**2)**(1/3)
-    a = fsolve(t,a0)[0]
+    def f(x):
+        tmp[0:3] = np.multiply(x,s[0:3])
+        tmp[3:6] = np.multiply(1/(x**0.5),s[3:6])
+        return T - time_period(tmp)
 
-    return a
+    return fsolve(f,1)[0]
 
 def tle_to_state(tle):
     """ This function converts from TLE elements to the position and velocity vector
@@ -73,17 +72,26 @@ def tle_to_state(tle):
     """
 
     # unload orbital elements array
-
-    sma = ntoa(tle[5],tle[0])
+    T = 86400/tle[5]
+    sma = (mu*(T/2/math.pi)**2)**(1/3)
     ecc = tle[2]  # eccentricity
     inc = tle[0]  # inclination
     argp = tle[3]  # argument of perigee
     raan = tle[1]  # right ascension of the ascending node
     tanom = MtoT(math.radians(tle[4]), ecc)  # we use mean anomaly(kep(5)) and the function MtoT to compute true anomaly (tanom)
     tanom = math.degrees(tanom)%360
 
-    print("SMA:",sma)
-    return kep_to_state(np.array([sma,ecc,inc,argp,raan,tanom]))
+    kep = np.array([sma,ecc,inc,argp,raan,tanom])
+    s = kep_to_state(kep)
+    x = scale(s,T)   # scales sma so that time period = T
+
+    kep[0] = kep[0]*x
+    print("Keplerian Elements:")
+    print(kep)
+    s[0:3] = np.multiply(x,s[0:3])
+    s[3:6] = np.multiply(1/(x**0.5),s[3:6])
+
+    return s
 
 def kep_to_state(kep):
     """ This function converts from keplerian elements to the position and velocity vector