In the orekit-python-wrapper tutorial Track_Corridor, the subsatellite point center is calculated as below

date    = currentState.getDate()
pvInert = currentState.getPVCoordinates()
t       = currentState.getFrame().getTransformTo(self.earth.getBodyFrame(), date)
p       = t.transformPosition(pvInert.getPosition())
v       = t.transformVector(pvInert.getVelocity())
center  = self.earth.transform(p, self.earth.getBodyFrame(), date)

which firstly transfers the satellite position pvInert to the Earth body frame (ITRF, remembering this by focusing on T for terrestrial) , and secondly transform the resulting Cartesian vector to (latitude, longitude, range) in WGS system used by self.earth ellipsoid.

For a validation as well as a simpler approach, I used the following function projectToGround() provided by Orekit natively:

subsatellite = self.earth.transform(
  self.earth.projectToGround(pvInert.getPosition(), date, currentState.frame),
  currentState.frame, date)

To verify they produce the same latitude and longitude, simply using

{subsatellite .__getattribute__(i) - center.__getattribute__(i) 
  for i in {'latitude', 'longitude', 'altitude'}}

Notice that altitude are different in two approaches. But since we are calculating subsatellite points only, the results are deemed as identical.

The two method has slightly numerical error on latitude, need to check the detailed implementation of projectToGround().

https://www.orekit.org/site-orekit-9.2/architecture/propagation.html

Taylor algebra
A very important implementation of the RealFieldElement interface is the DerivativeStructure class, which in addition to compute the result of the canonical operation (add, multiply, sin, atanh…) also computes its derivatives, with respect to any number of variables and to any derivation order. If for example a user starts a computation with 6 canonical variables px, py, pz, vx, vy, vz to represent an initial state and then performs a propagation. At the end for all produced results (final position, final velocity but also geodetic altitude with respect to an ellipsoid body or anything that Orekit computes), then for these results one can retrieve its partial derivatives up to the computed order with respect to the 6 canonical variables. So if for example in a step handler you compute a geodetic altitude h, you also have ∂³h/∂px²∂vz or any of the 84 components computed at order 3 for each value (1 value, 6 first order derivatives, 14 second order derivatives, 56 third order derivatives). The DerivativeStructure class also provides Taylor expansion, which allow to extrapolate the result accurately to close values. This is an implementation of Taylor Algebra. Its two main uses in space flight dynamics are

https://www.orekit.org/site-orekit-9.3/architecture/estimation.html

Kalman filter
A simplified model expanding an error ellipsoid with principal axes in the along-track, across-track and normal directions and using polynomials lengths for the ellipsoid axes is expected to be included in the library. Finding the right coefficients for the polynomials will remain the responsibility of user. One way to tune it properly for a family of orbits is to use a first run of Taylor algebra during the mission analysis phase, to monitor how the high order uncertainties evolve, and to fit an ellipsoid on this evolution. Then a Kalman filter using a covariance matrix provider with such a polynomial ellipsoid model will be more realistic than a basic constant matrix.

似乎没有找到相应的代码，猜测是没有实现。

下一步工作是否需要依赖这部分功能？（不确定是否有更好的办法实现uncertainty propagation，暂时收尾，到时候再考虑。）

[Orekit Developers] Taylor differential algebra, https://www.orekit.org/mailing-list-archives/orekit-developers/msg00236.html

The approach we followed was to create generic Field-based propagators, so users can do things like:

FieldPropagator<SomeField> propagator = new FieldNumericalPropagator<SomeField>(...);

(or some other supported analytical propagator if desired) and then to use it just like a regular double-based propagator.

In order to use Taylor differential algebra, users just have to use DerivativeStructure as the field for propagated elements.

This mean you propagate uncertainties to a higher order than classical state transition matrices (or covariance matrices) which are inherently first order only.

Of course, the initial propagation is much slower, especially if you have high order, but the overhead is done only once, so the complete cost is easily paid for for any realistic Monte-Carlo size.

The field output is mainly a FieldSpacecraftState, but of course user can also use derived parameters (position/velocity/acceleration, in any frame, inertial or rotating, geodetic coordinates …).

Andrea already implemented Keplerian, Eckstein-Hechler and numerical propagators with a bunch of force models. We may add SGP4/SDP4 also, it is quite straightforward to do. （有添加SGP4支持吗？暂时没有找到……）

A side-effect of our approach is that we can use it to implement other things apart from Taylor differential algebra (high accuracy, interval arithmetics, …).

Re: [Orekit Developers] Taylor differential algebra, https://www.orekit.org/mailing-list-archives/orekit-developers/msg00259.html

[Orekit Users] Taylor algebra and force model parameters, https://www.orekit.org/mailing-list-archives/orekit-users/msg00124.html

to analyze the effects on the orbit of the dispersion of force model parameters. Propagating uncertainties on the orbit knowledge from OD uncertainties is quite straightforward with “FieldPropagation” but I don't see the way to analyze the effect of the uncertainty of a force model parameter with taylor algebra in Orekit.

Re: [Orekit Users] Taylor algebra and force model parameters, https://www.orekit.org/mailing-list-archives/orekit-users/msg00119.html

Orekit-tutorials mentioning “Taylor”: https://gitlab.orekit.org/search?search=taylor&project_id=35&group_id=3
https://gitlab.orekit.org/orekit/orekit-tutorials/-/blob/master/src/main/java/org/orekit/tutorials/propagation/FieldPropagation.java

COMMERCIAL SPACE TRACKING SERVICE FOR SMALL SATELLITES, https://digitalcommons.usu.edu/cgi/viewcontent.cgi?filename=0&article=4412&context=smallsat&type=additional

Used Orekit open source orbital dynamics

In Apache Common Math, http://commons.apache.org/proper/commons-math/userguide/analysis.html#a4.7_Differentiation

The second method is to write a classical UnivariateFunction and to pass it to an existing implementation of the UnivariateFunctionDifferentiator interface to retrieve a differentiated version of the same function.

Apache Commons Math provides one implementation of the UnivariateFunctionDifferentiator interface: FiniteDifferencesDifferentiator. This class creates a wrapper that will call the user-provided function on a grid sample and will use finite differences to compute the derivatives. It takes care of boundaries if the variable is not defined on the whole real line.

Another implementation of the UnivariateFunctionDifferentiator interface is under development in the related project Apache Commons Nabla. This implementation uses automatic code analysis and generation at binary level. However, at time of writing (end 2012), this project is not yet suitable for production use.

The one in Hipparchus is more advanced, with many old constructor APIs in Common Math deprecated and replaced by a DSFactory paradigm. See https://hipparchus.org/hipparchus-core/analysis.html#Differentiation.

Hipparchus also provides different integrators for solving Ordinary Differential Equations. All of them are already able to integrate systems of differentials equations based on DSs. [1]

The differentiation framework is based on the doubly recursive multivariate method [2] for the computation of DS, with an additional compilation step that folds the computation rules found through recursion into extremely fast single pass iterative rules using indirection tables. [1]

The main advantages of the automatic differentiation is in the between of the two precedent ways. It has the precision of an analytical method, being of course slower in the computations, and the adaptability of the numerical, without problems on the boundaries, and the different scale of its values are irrelevant. [1]

The translation from the primitive double number-based Orekit to the DS-based Orekit has covered all the different types of orbits representation, some of the attitudes, some of the force models, all the analytical propagators and the numerical propagator.
The force models for the numerical propagator already translated include the main body gravity field, with all the tesseral and zonal terms, the isotropic drag, the third body attraction and the solar radiation pressure. If needed also all the other force models on Orekit may be translated, but in order to ease the validation of the global project it has been chosen to translate only in case of need. The remaining force models will be translated in the future (solid tides, ocean tides, relativity, nonisotropic non-conservative forces, ). [1]

Once the DS had been implemented and tested in Orekit, it had its firsts applications. It has been used in an internal project of C-S to validate the orbit restitution of Orekit. The goal of this application was to see how the error on some of the orbital elements would evolve in 48h.

To show the performances of the DS propagation we chose to compare it with the propagation time of a Real Based propagation to show how the DS propagation slows the procedure. The table 1 is showing that in case of a MC simulation its worth using the DS propagation if the number of samples needed for theMontecarlo simulation, is above the number wrote in the table.

Table-1 from [1].

According to the developer Luc's reply, to my understanding, Hipparchus is newer and better maintained, at least for the purpose of Orekit. Not much needed to be concerned if I only need it for Orekit.

Field<T> is the mathematical concept field. Field_(mathematics)

Field<T>定义了一个有加法和乘法封闭的域，其中的元素是FieldElement<T>。

CalculusFieldElement<T> – extended by → RealFieldElement<T> – extended by → T which is the type tinput of FieldAbsoluteDate<>

It is an analytical propagator considering J2, J3, J4, J5, J6 potential zonal coefficients, where higher order Jn parameters can be set to zero. See https://www.orekit.org/site-orekit-9.1/architecture/propagation.html#Eckstein-Hechler_propagation