Daily repeat groundtrack Mars orbit

[Meuge]

• Added solar rotation period constants for some bodies .
• Integrate the files, and the constants.
• Documentation
• Merge tests.
• Replace the constant "Mean rotation rate from the concerning body around the sun, n_sol"
• Add some tests for other attractors.

[codecov]

Codecov Report

Merging #1002 into master will decrease coverage by 0.21%.
The diff coverage is 84.14%.

@@            Coverage Diff             @@
##           master    #1002      +/-   ##
==========================================
- Coverage   89.52%   89.30%   -0.22%     
==========================================
  Files          74       74              
  Lines        3941     4104     +163     
  Branches      343      369      +26     
==========================================
+ Hits         3528     3665     +137     
- Misses        322      335      +13     
- Partials       91      104      +13     

Impacted Files	Coverage Δ
src/poliastro/util.py	84.96% <81.48%> (-15.04%)	⬇️
src/poliastro/twobody/orbit.py	87.27% <89.28%> (+0.22%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c78984e...3a3a7a7. Read the comment docs.

[astrojuanlu]

Gave a first pass and left some comments, will give a second pass tomorrow!

[astrojuanlu]

I wonder if there is a way to derive the solar rotational period given the rotational period and the orbital period around the Sun. Am I missing something?

[Meuge]

Maybe it's possible, but I don't have the right knowledge to do so. The paper sets that variable as a constant.

[Meuge]

Doing some research I found here the relation between Psol, Psid and the Period of the attractor around the sun

And Period of the attractor around the sun you can obtain it from

detailed in the book An introduction to modern astrophysics. 2nd edition. Pearson
So there is no need to add this as a constant

[astrojuanlu]

Any reference of where this algorithm came from?

[Meuge]

This is explained in the following comment

[astrojuanlu]

Is this algorithm better than minimizing |y - f(x)| using SciPy minimize_scalar?

[Meuge]

Juan Luis, I've been doing some research and also tried using scipy.optimize.minimize_scalar.
The algorithm applied, as I have said before, it's linked to the LMS algorithm, particularly I used a VSS-LMS (Variable Step Linear Mean Square) variant linked to the Harris et al version.
The VSSLMS is an enhanced version of LMS that improves the convergence rate and mean-square error.
It is a very well known algorithm, one of the principal competitors to the Kalman Filter performance and accuracy, but with the need for less computation power and easier maintainability.
Here you can see that they are trying to link both algorithms to get a better hybrid solution.
You can see the algorithm is explained here
The algorithm implemented is one of the most simple cases, because L=1, α = 2, m0 = 1, m1 = 2, u = 1, d = Pq.

After consulting with some experts in signal processing academics of my university, they explained to me that you can think of this problem as a constant signal of u = [1....1] and with a fixed target/desired signal d = [Pq, Pq, Pq...]. With those two concepts in mind, you will find that W will result in the solution of the problem. Another interesting point that they have advised me, is that due to d and u are fixed, you can avoid w = w + µ * u * e multiplication, and just use the sign version w = w + µ * sign(u * e) => w = w +/- µ. In this particular case, this will result in less computation and faster convergence.

The problem when I was solving the issue, it's that when you use scipy.optimize.minimize_scalar, it uses one of three methods Brent, Bounded, or Golden, being Brent the fastest of the three. Nonetheless, when you don't define a bracket for minimize_scalar, it needs to find one because all of the three algorithms need a bracket with two different signs. So in order to do it, it uses the bracket method to find a viable bracket while losing time and computing power. On the other hand, the VSS-LMS algorithm achieves the solution without the necessity of finding a bracket.

The next point that it's necessary to explain, is that the Brent Algorithm it's a hybrid algorithm that links the secant and the bisect method. Brent method attempts to exploit some advantages of the secant but there are some known difficulties related to the convergence rate. It tries to minimize this convergence defect by switching the secant method to the bisect method and continues from there.
However, like newton raphson, in secant and brent algorithm, you need to know a little bit about the function you want to optimize because there are cases that may not converge
In the last slide, you can see that the problem doesn't converge for f(x) = 1/x. And that's the reason why I have some issues with EQ 9 from the Noreen paper, because nt(a) ~ 1/a^(7/2).

When you analyze the VSSLMS, you can realize that the algorithm performs quite similar to the bisect algorithm. First when trying to approach the d value (Pq) by increasing the step 2^n and then when w > d, VSSLMS starts pivoting around d, and reducing its step every time w is greater or smaller than d, so the step is being reduced by 2^-n. In this case, you do not need to find a bracket, and the problem has a complexity of log2(n).
This is a great advantage, first of all, for the simplicity of the algorithm, the less usage of computational cost, and then for reducing the time needed to converge.

[Meuge]

[Meuge]

I continued doing some research and tried to use other root-finding algorithms. I've tested Scipy algorithms and I reached the following conclusions:
Scipy.optimize.minimize_scalar: it's not viable to use it since it finds the local minimums and not the roots, in the final attached picture you'll see that the function gets stuck in a local minimum, so we don't end up with the right solution.

What we can use it's, Scipy.optimize.root_scalar, specifically should be Bisect, Brentq, Brenth, Ridder, Toms748, and Secant. Newton Raphson and Halley use the d/da (Psol/Q - 2π/nt(a)), but I didn't implement them because it won't get much difference from the Secant method and the complexity is greater. If you want to implement the methods as well I can get into it.

The problem with Scipy.optimize.root_scalar, is that you need to give them a bracket where the function has a different sign. So by analyzing the total mean motion and the nodal period you can find an interval for a:

As for the elapsed time, LMS takes a little bit more of time since it has under consideration Astropy Quantities, and that the algorithm knows nothing about the bracket, which it's a big win when you don't know anything about the function.

Notice that minimize_scalar doesn't solve the problem and stays in a local minimum.

[astrojuanlu]

Thanks for the thorough explanations about *LMS in the PR @Meuge ! Being an algorithm this powerful I expect that we would try to contribute it to SciPy at some point, like we did with the DOP853 integration scheme couple of years ago.

As a last chance to simplify the approach, I was still not satisfied with having two non-linear equations, so by doing some algebraic manipulations on equations 15 and 18 from Noreen I arrived to this 7th order polynomial for a:

In [1]: from sympy import symbols, Eq, solve, sin, cos, sqrt, Rational 
   ...: from sympy.abc import i, a 
   ...: from sympy import init_printing 
   ...: init_printing() 
   ...:  
   ...: R_M, μ, J_2, n_sol = symbols("R_M, μ, J_2, n_sol") 
   ...: P_Omega, P_K = symbols("P_Ω, P_K") 
   ...:  
   ...: eq1 = Eq(cos(i), -(2 * a ** Rational(7,2) * n_sol) / (3 * J_2 * sqrt(μ) * R_M**2)) 
   ...: eq2 = Eq(P_Omega, P_K / (1 + (3 * R_M**2 * J_2 * (3 - 4 * sin(i) ** 2) / (2 * a**2)))) 
   ...:  
   ...: eq2_subs = eq2.subs(i, solve(eq1, i)[1]).simplify() 
   ...:  
   ...: _den = eq2_subs.rhs.as_numer_denom()[1] 
   ...: eq2_subs_noden = Eq(eq2_subs.lhs * _den, eq2_subs.rhs * _den) 
   ...:  
   ...: eq2_subs_noden.expand().collect(a) 
   ...:                                                                                                                                                                                                                                       
Out[1]: 
      2        4                 2  2             7     2               2  2  
- 9⋅J₂ ⋅P_Ω⋅R_M ⋅μ + 6⋅J₂⋅P_Ω⋅R_M ⋅a ⋅μ + 16⋅P_Ω⋅a ⋅nₛₒₗ  = 6⋅J₂⋅P_K⋅R_M ⋅a ⋅μ

Even if this must be solved numerically, do you think we can use this to our advantage somehow?

[astrojuanlu]

Also, the fact that we now have tests will make everything much easier 👏

[Meuge]

Thanks for the thorough explanations about *LMS in the PR @Meuge ! Being an algorithm this powerful I expect that we would try to contribute it to SciPy at some point, like we did with the DOP853 integration scheme couple of years ago.

As a last chance to simplify the approach, I was still not satisfied with having two non-linear equations, so by doing some algebraic manipulations on equations 15 and 18 from Noreen I arrived to this 7th order polynomial for a:

In [1]: from sympy import symbols, Eq, solve, sin, cos, sqrt, Rational 
   ...: from sympy.abc import i, a 
   ...: from sympy import init_printing 
   ...: init_printing() 
   ...:  
   ...: R_M, μ, J_2, n_sol = symbols("R_M, μ, J_2, n_sol") 
   ...: P_Omega, P_K = symbols("P_Ω, P_K") 
   ...:  
   ...: eq1 = Eq(cos(i), -(2 * a ** Rational(7,2) * n_sol) / (3 * J_2 * sqrt(μ) * R_M**2)) 
   ...: eq2 = Eq(P_Omega, P_K / (1 + (3 * R_M**2 * J_2 * (3 - 4 * sin(i) ** 2) / (2 * a**2)))) 
   ...:  
   ...: eq2_subs = eq2.subs(i, solve(eq1, i)[1]).simplify() 
   ...:  
   ...: _den = eq2_subs.rhs.as_numer_denom()[1] 
   ...: eq2_subs_noden = Eq(eq2_subs.lhs * _den, eq2_subs.rhs * _den) 
   ...:  
   ...: eq2_subs_noden.expand().collect(a) 
   ...:                                                                                                                                                                                                                                       
Out[1]: 
      2        4                 2  2             7     2               2  2  
- 9⋅J₂ ⋅P_Ω⋅R_M ⋅μ + 6⋅J₂⋅P_Ω⋅R_M ⋅a ⋅μ + 16⋅P_Ω⋅a ⋅nₛₒₗ  = 6⋅J₂⋅P_K⋅R_M ⋅a ⋅μ

Even if this must be solved numerically, do you think we can use this to our advantage somehow?

I understand what you are trying to achieve, here are the steps to substitute inc(a). However, you can see that the algebraic complexity increases exponentially.

You can try to solve PΩ equation, but you can see that the terms involved in the final equation are not simple to determine "a".
We can obtain a polynomial approximation but we cannot assure that we reach a proper solution of "a" or that "a" has enough terms in order to assume that the error is too small to be considered relevant.
This problem is similar in a way to:

[astrojuanlu]

I think part of the confusion comes because we are discussing different problems. I was using equations 15 and 18 which solve the "circular Sun synchronous orbits", while I think you're solving a more general case. It's made more confusing by the fact that equations 5 to 8 are under the "circular equatorial orbits" of the paper (therefore not Sun-sync inclination, but inc = 0), however this assumption is only introduced in equations 9 and 12-14.

My superficial review of the code was not enough and I still need to dive deeper into this, sorry for slowing you down 🙏