update midpoint function

pages/inference.py

@@ -130,7 +130,6 @@ def high_level(option_time, slider_sample_orbit, progress_bar):
 def midpoint(option_time, slider_sample_orbit, progress_bar):
     time.sleep(0.1)
     mid_point_total_start_time = time.time()
-    mid_point_30000_start_time = time.time()
 
     midpointhelper = MidPoint(j=slider_sample_orbit)
     j, h, b, n, x, y, z, xa, ya, za, px, py, pz, pxa, pya, pza = midpointhelper.initial()
@@ -147,25 +146,26 @@ def midpoint(option_time, slider_sample_orbit, progress_bar):
 
     original_mid_point_data = []
 
+    mid_point_30000_start_time = time.time()
+
     while t < float(option_time):
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a1 * h, x, y, z, px, py, pz, b)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a2 * h, x, y, z, px, py, pz, b)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a1 * h, x, y, z, px, py, pz, b)
 
         t = t + h
 
         if jn % 10 == 0:
             original_mid_point_data.append([b, x, y, z, px, py, pz])
+            # Update progress bar
+            progress_percentage = int((current_iteration / total_iterations) * 100)
+            progress_bar.progress(progress_percentage)
         if jn == 300000:
             mid_point_30000_end_time = time.time()
             mid_point_30000_execute_time = mid_point_30000_end_time - mid_point_30000_start_time
             mid_point_2000_start_time = time.time()
         jn = jn + 1
 
-        # Update progress bar
-        progress_percentage = int((current_iteration / total_iterations) * 100)
-        progress_bar.progress(progress_percentage)
         current_iteration += 1
 
     #mid_point_df.to_excel('mid_point_df_output.xlsx', index=False)
@@ -277,10 +277,9 @@ def mid_point_lstm(slider_sample_orbit, lstm_progress_bar):
     mid_point_start_time = time.time()
 
     while t < float(30000):
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa,pya, pza)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a1 * h, x, y, z, px, py, pz, b)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a2 * h, x, y, z, px, py, pz, b)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a1 * h, x, y, z, px, py, pz, b)
 
         t = t + h
 
@@ -429,10 +428,9 @@ def mid_point_tcn(slider_sample_orbit, tcn_progress_bar):
 
     # Perform classical method prediction for the initial segment
     while t < float(30000):
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a2, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.symplectic(h * a1, x, y, z, px, py, pz, xa, ya,za, pxa, pya, pza, b)
-        x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza = midpointhelper.rejust(x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a1 * h, x, y, z, px, py, pz, b)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a2 * h, x, y, z, px, py, pz, b)
+        x, y, z, px, py, pz = midpointhelper.implecitsymplectic(a1 * h, x, y, z, px, py, pz, b)
 
         t = t + h
 
@@ -490,18 +488,18 @@ with st.sidebar:
     st.write(f'Total Time Step: {option_time}')
     options_method = st.multiselect(
         'Compared Methods',
-        ['Low-Level', 'High-Level', 'Midpoint', 'Low-Level with LSTM', 'Low-Level with TCN', 'Midpoint with LSTM', 'Midpoint with TCN'],
-        ['Low-Level'])
+        ['EPS', 'Midpoint', 'EPS with LSTM', 'EPS with TCN', 'Midpoint with LSTM', 'Midpoint with TCN'],
+        ['EPS'])
     btn_go = st.button("Go", type="primary", use_container_width=True)
 
 if btn_go:
-    if 'Low-Level' in options_method:
+    if 'EPS' in options_method:
         with container1:
-            st.write('Low Level Progress Bar')
+            st.write('EPS Progress Bar')
             low_level_progress_bar = st.progress(0)
         low_level_30000_time, low_level_2000_time, low_level_total_time, low_level_result = low_level(option_time, slider_sample_orbit, low_level_progress_bar)
         with container2:
-            st.table(pd.DataFrame({'Model':"Low Level", '30000 Time Steps (s)': [low_level_30000_time], '2000 Time Steps (s)': [low_level_2000_time], 'Total Time (s)': [low_level_total_time]}))
+            st.table(pd.DataFrame({'Model':"EPS", '30000 Time Steps (s)': [low_level_30000_time], '2000 Time Steps (s)': [low_level_2000_time], 'Total Time (s)': [low_level_total_time]}))
     if 'High-Level' in options_method:
         with container1:
             st.write('High Level Progress Bar')
@@ -516,20 +514,20 @@ if btn_go:
         mid_point_30000_time, mid_point_2000_time, mid_point_total_time, mid_point_result = midpoint(option_time, slider_sample_orbit, mid_point_progress_bar)
         with container2:
             st.table(pd.DataFrame({'Model':"Midpoint", '30000 Time Steps (s)': [mid_point_30000_time], '2000 Time Steps (s)': [mid_point_2000_time], 'Total Time (s)': [mid_point_total_time]}))
-    if 'Low-Level with LSTM' in options_method:
+    if 'EPS with LSTM' in options_method:
         with container1:
-            st.write('Low Level LSTM Progress Bar')
+            st.write('EPS LSTM Progress Bar')
             low_level_lstm_progress_bar = st.progress(0)
         lstm_30000_time, lstm_2000_time, lstm_total_time, lstm_result = low_level_lstm(slider_sample_orbit, low_level_lstm_progress_bar)
         with container2:
-            st.table(pd.DataFrame({'Model':"Low Level + LSTM", '30000 Time Steps (s)': [lstm_30000_time], '2000 Time Steps (s)': [lstm_2000_time], 'Total Time (s)': [lstm_total_time]}))
-    if 'Low-Level with TCN' in options_method:
+            st.table(pd.DataFrame({'Model':"EPS + LSTM", '30000 Time Steps (s)': [lstm_30000_time], '2000 Time Steps (s)': [lstm_2000_time], 'Total Time (s)': [lstm_total_time]}))
+    if 'EPS with TCN' in options_method:
         with container1:
-            st.write('Low Level TCN Progress Bar')
+            st.write('EPS TCN Progress Bar')
             low_level_tcn_progress_bar = st.progress(0)
         tcn_30000_time, tcn_2000_time, tcn_total_time, tcn_result = low_level_tcn(slider_sample_orbit, low_level_tcn_progress_bar)
         with container2:
-            st.table(pd.DataFrame({'Model':"Low Level + TCN", '30000 Time Steps (s)': [tcn_30000_time], '2000 Time Steps (s)': [tcn_2000_time], 'Total Time (s)': [tcn_total_time]}))
+            st.table(pd.DataFrame({'Model':"EPS + TCN", '30000 Time Steps (s)': [tcn_30000_time], '2000 Time Steps (s)': [tcn_2000_time], 'Total Time (s)': [tcn_total_time]}))
     if 'Midpoint with LSTM' in options_method:
         with container1:
             st.write('Midpoint LSTM Progress Bar')

utils/midpoint.py

@@ -1,3 +1,4 @@
+import numpy as np
 
 class MidPoint():
 
@@ -91,74 +92,114 @@ class MidPoint():
 
         e = ht + hv + h1pn
 
-        vnpx=px
-        v1pnpx=4*px*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (px*(n + 3) + (n*x*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
-        vpx=vnpx+v1pnpx
-
-        vnpy=py
-        v1pnpy=4*py*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (py*(n + 3) + (n*y*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
-        vpy=vnpy+v1pnpy
-
-        vnpz=pz
-        v1pnpz=4*pz*((3*n)/8 - 1/8)*(px**2 + py**2 + pz**2) - (pz*(n + 3) + (n*z*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2)))/(x**2 + y**2 + z**2)**(1/2))/(x**2 + y**2 + z**2)**(1/2)
-        vpz=vnpz+v1pnpz
-
-        vnx=x/(x**2 + y**2 + z**2)**(3/2)
-        v1pnx=(x*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*x)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((px*x**2)/(x**2 + y**2 + z**2)**(3/2) -px/(x**2 + y**2 + z**2)**(1/2) + (py*x*y)/(x**2 + y**2 + z**2)**(3/2) + (pz*x*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
-        vx=vnx+v1pnx
-
-        vny=y/(x**2 + y**2 + z**2)**(3/2)
-        v1pny=(y*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*y)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((py*y**2)/(x**2 + y**2 + z**2)**(3/2) -py/(x**2 + y**2 + z**2)**(1/2) + (px*x*y)/(x**2 + y**2 + z**2)**(3/2) + (pz*y*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
-        vy=vny+v1pny
-
-        vnz=z/(x**2 + y**2 + z**2)**(3/2)
-        v1pnz=(z*((n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 +y**2 + z**2)**(1/2))**2)/2 + ((n + 3)*(px**2 + py**2 + pz**2))/2))/(x**2 + y**2 + z**2)**(3/2) -(4*z)/(2*x**2 + 2*y**2 + 2*z**2)**2 + (n*((px*x)/(x**2 + y**2 + z**2)**(1/2) + (py*y)/(x**2 + y**2 + z**2)**(1/2) + (pz*z)/(x**2 + y**2 + z**2)**(1/2))*((pz*z**2)/(x**2 + y**2 + z**2)**(3/2) -pz/(x**2 + y**2 + z**2)**(1/2) + (px*x*z)/(x**2 + y**2 + z**2)**(3/2) + (py*y*z)/(x**2 + y**2 +z**2)**(3/2)))/(x**2 + y**2 + z**2)**(1/2)
-        vz=vnz+v1pnz
-
-        return vx,vy,vz,vpx,vpy,vpz,e
-
-    def rejust(self, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza):
-
-        x = (x + xa) / 2
-        y = (y + ya) / 2
-        z = (z + za) / 2
-
-        px = (px + pxa) / 2
-        py = (py + pya) / 2
-        pz = (pz + pza) / 2
-        xa = x
-        ya = y
-        za = z
-        pxa = px
-        pya = py
-        pza = pz
-
-        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza
-
-    def symplectic(self, h, x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza, b):
-
-        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
-        x = x + h / 2 * vpx
-        y = y + h / 2 * vpy
-        z = z + h / 2 * vpz
-        pxa = pxa - h / 2 * vxa
-        pya = pya - h / 2 * vya
-        pza = pza - h / 2 * vza
-
-        vx, vy, vz, vpxa, vpya, vpza, e = self.f(x, y, z, pxa, pya, pza, b)
-        xa = xa + h * vpxa
-        ya = ya + h * vpya
-        za = za + h * vpza
-        px = px - h * vx
-        py = py - h * vy
-        pz = pz - h * vz
-
-        vxa, vya, vza, vpx, vpy, vpz, e = self.f(xa, ya, za, px, py, pz, b)
-        x = x + h / 2 * vpx
-        y = y + h / 2 * vpy
-        z = z + h / 2 * vpz
-        pxa = pxa - h / 2 * vxa
-        pya = pya - h / 2 * vya
-        pza = pza - h / 2 * vza
-
-        return x, y, z, px, py, pz, xa, ya, za, pxa, pya, pza
+        vnpx = px
+        v1pnpx = 4 * px * ((3 * n) / 8 - 1 / 8) * (px ** 2 + py ** 2 + pz ** 2) - (px * (n + 3) + (n * x * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2))) / (x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                               1 / 2)) / (
+                             x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        vpx = vnpx + v1pnpx
+
+        vnpy = py
+        v1pnpy = 4 * py * ((3 * n) / 8 - 1 / 8) * (px ** 2 + py ** 2 + pz ** 2) - (py * (n + 3) + (n * y * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2))) / (x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                               1 / 2)) / (
+                             x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        vpy = vnpy + v1pnpy
+
+        vnpz = pz
+        v1pnpz = 4 * pz * ((3 * n) / 8 - 1 / 8) * (px ** 2 + py ** 2 + pz ** 2) - (pz * (n + 3) + (n * z * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2))) / (x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                               1 / 2)) / (
+                             x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        vpz = vnpz + v1pnpz
+
+        vnx = x / (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)
+        v1pnx = (x * ((n * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) ** 2) / 2 + (
+                                  (n + 3) * (px ** 2 + py ** 2 + pz ** 2)) / 2)) / (x ** 2 + y ** 2 + z ** 2) ** (
+                            3 / 2) - (4 * x) / (2 * x ** 2 + 2 * y ** 2 + 2 * z ** 2) ** 2 + (n * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) * ((px * x ** 2) / (
+                    x ** 2 + y ** 2 + z ** 2) ** (3 / 2) - px / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * x * y) / (
+                                                                                                  x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                                  3 / 2) + (
+                                                                                                  pz * x * z) / (
+                                                                                                  x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                                  3 / 2))) / (
+                            x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        vx = vnx + v1pnx
+
+        vny = y / (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)
+        v1pny = (y * ((n * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) ** 2) / 2 + (
+                                  (n + 3) * (px ** 2 + py ** 2 + pz ** 2)) / 2)) / (x ** 2 + y ** 2 + z ** 2) ** (
+                            3 / 2) - (4 * y) / (2 * x ** 2 + 2 * y ** 2 + 2 * z ** 2) ** 2 + (n * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) * ((py * y ** 2) / (
+                    x ** 2 + y ** 2 + z ** 2) ** (3 / 2) - py / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (px * x * y) / (
+                                                                                                  x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                                  3 / 2) + (
+                                                                                                  pz * y * z) / (
+                                                                                                  x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                                  3 / 2))) / (
+                            x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        vy = vny + v1pny
+
+        vnz = z / (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)
+        v1pnz = (z * ((n * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) ** 2) / 2 + (
+                                  (n + 3) * (px ** 2 + py ** 2 + pz ** 2)) / 2)) / (x ** 2 + y ** 2 + z ** 2) ** (
+                            3 / 2) - (4 * z) / (2 * x ** 2 + 2 * y ** 2 + 2 * z ** 2) ** 2 + (n * (
+                    (px * x) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (py * y) / (x ** 2 + y ** 2 + z ** 2) ** (
+                        1 / 2) + (pz * z) / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2)) * ((pz * z ** 2) / (
+                    x ** 2 + y ** 2 + z ** 2) ** (3 / 2) - pz / (x ** 2 + y ** 2 + z ** 2) ** (1 / 2) + (px * x * z) / (
+                                                                                                  x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                                  3 / 2) + (
+                                                                                                  py * y * z) / (
+                                                                                                  x ** 2 + y ** 2 + z ** 2) ** (
+                                                                                                  3 / 2))) / (
+                            x ** 2 + y ** 2 + z ** 2) ** (1 / 2)
+        vz = vnz + v1pnz
+        return vx, vy, vz, vpx, vpy, vpz, e
+
+    def implecitsymplectic(self, h, x, y, z, px, py, pz, b):
+
+        x1 = x
+        y1 = y
+        z1 = z
+        px1 = px
+        py1 = py
+        pz1 = pz
+        num = 0
+        d = 1.1
+        while (d > 1e-16 and num < 1000):
+            px2 = px
+            py2 = py
+            pz2 = pz
+            x2 = x
+            y2 = y
+            z2 = z
+            x = (x1 + x) / 2
+            y = (y1 + y) / 2
+            z = (z1 + z) / 2
+            px = (px1 + px) / 2
+            py = (py1 + py) / 2
+            pz = (pz1 + pz) / 2
+            vx, vy, vz, vpx, vpy, vpz, e = self.f(x, y, z, px, py, pz, b)
+            x = x1 + h * vpx
+            y = y1 + h * vpy
+            z = z1 + h * vpz
+            px = px1 - h * vx
+            py = py1 - h * vy
+            pz = pz1 - h * vz
+            d = np.sqrt(np.abs(x - x2) + np.abs(y - y2) + np.abs(z - z2) + np.abs(px - px2) + np.abs(py - py2) + np.abs(
+                pz - pz2))
+            num += 1
+            #print(num, d)
+        return x, y, z, px, py, pz