User Manual 4.7 Celestial bodies
Introduction
Scope
This section presents the Trigonometric Polynomials implemented in PATRIUS as well as Fourier Series.
Javadoc
The trigonometric polynomial (and related) objects are available in the package fr.cnes.sirius.patrius.math.analysis.polynomials
and the FFT algorithms in fr.cnes.sirius.patrius.math.transform
Library | Javadoc |
---|---|
Patrius | Package fr.cnes.sirius.patrius.math.geometry.analysis.polynomials |
Patrius | Package fr.cnes.sirius.patrius.math.transform |
Links
Useful resources for this theme can be found here :
Useful Documents
None.
Package Overview
Features Description
Trigonometric Polynomials
The trigonometric polynomials are defined as :
[math]\forall t \in \mathbb{R}, P(t) = a_0 + \sum_{k=0}^n \left( a_k \cos(kt) + b_k \sin(kt) \right)[/math]
And the primitive of such a function is :
[math]\forall t \in \mathbb{R}, P(t) = a_0 t + c^{te} + \sum_{k=0}^n \left( {-b_k \over k} \cos(kt) + {a_k \over k} \sin(kt) \right)[/math]
Classes that represent these functions are described in this section.
Fourier Series
Let [math]f[/math] be a periodic, real variable, real function with period [math]T[/math]. The partial sums of the Fourier Series of [math]f[/math] are given by
[math] \forall N \in \mathbb{R}^{+*}, \, \left(S_Nf\right)(x) = a_0 + \sum_{n=1}^N\left(a_n\cos\left(n t {2\pi\over T}\right) + b_n\sin\left(n t {2\pi\over T}\right)\right) [/math]
The Fourier coefficients (n > 0) of [math]f[/math] are given by :
[math] a_0 = {1\over T}\int_{-T/2}^{T/2} \! f(t) \, \mathrm{d}t [/math]
[math] a_n = {2\over T}\int_{-T/2}^{T/2} \! f(t) \cos\left(n t {2\pi\over T}\right) \, \mathrm{d}t [/math]
[math] b_n = {2\over T}\int_{-T/2}^{T/2} \! f(t) \sin\left(n t {2\pi\over T}\right) \, \mathrm{d}t [/math]
Classes that allow decomposing functions into finite Fourier series are described.
Fast Fourier Transforms
A fast Fourier transform (FFT) is an algorithm to compute the discrete Fourier transform (DFT) and its inverse. Fourier analysis converts time (or space) to frequency and vice versa. An FFT rapidly computes such transformations by factorizing the DFT matrix into a product of sparse (mostly zero) factors.
Let [math]( X_0, ..., X_{N-1})[/math] be complex numbers. The DFT is defined by the formula
[math]X_k = \sum^{N-1}_{n=0} x_n e^{-2i\pi k n/N}, \quad k=\left\{0, ..., N-1\right\}.[/math]
The class FastFourierTransformer.java
owns two methods
transform(final double[] f, final TransformType type)
transform(final Complex[] f, final TransformType type)
where TransformType
is an enumeration that defines the type of transform which is to be computed. It can be FORWARD or INVERSE.
Depending on the size, this method will use a radix-2 algorithm if the size is a power-of-two, and the Bluestein algorithm otherwise. Both algorithms are explained above.
The Cooley–Tukey & radix-2 the algorithm
By far the most commonly used FFT algorithm. This is a divide and conquer algorithm that recursively breaks down a DFT of any composite size [math]N = N_1N_2[/math] into many smaller DFTs of sizes [math]N_1[/math] and [math]N_2[/math], along with [math]O(N)[/math] multiplications by complex roots of unity.
The radix-2 algorithme is the most common use of FFT algorithm and is based on the Cooley–Tukey algorithm. It divides the transform into two pieces of size [math]N/2[/math] at each step, and is therefore limited to power-of-two sizes.
The Bluestein algorithm
It is commonly called the chirp z-transform algorithm. It computes the discrete Fourier transform (DFT) of arbitrary sizes (including prime sizes) by re-expressing the DFT as a convolution. For more information, see http://en.wikipedia.org/wiki/Bluestein%27s_FFT_algorithm
Getting Started
Usage of a trigonometric polynomial
- To create the polynomial [math]P(x) =-2 + \cos(x) + 2\sin(3x)[/math], use the following code snippet :
double a0 =-2; double[] a = new double[] {1, 0, 0}; double[] b = new double[] {0, 0, 2}; TrigonometricPolynomialFunction myPol2 = new TrigonometricPolynomialFunction(a0, a, b);
- To get the value of the second derivative at [math]x = 5[/math], use the following snippet :
double result = myPol2.value(2,5); // equivalent to double result = myPol2.polynomialDerivative(2).value(5);
- To get the primitive :
// With an integration constant c = 1 : TrigonometricPolynomialPrimitive result = myPol2.polynomialPrimitive(1.);
- To multiply two trigonometric polynomials [math]p1[/math] and [math]p2[/math] :
TrigonometricPolynomialFunction result = p1.multiply(p2);
Decomposing a UnivariateFunction into a Fourier Series
Given a user function with period T :
UnivariateFunction f = new UnivariateFunction () { public double value(double x) { // ... return result; } };
The user must build an instance of the FourierDecompositionEngine (using an integrator) :
FourierDecompositionEngine engine = new FourierDecompositionEngine( integrator );
And pass the function, period and approximation order to it.NOTE : It is left up to the user to make sure that the specified period is coherent with the function as no internal check is conducted.
engine.setOrder(10); engine.setFunction(f, T);
Then, the user can call the decompose( )
method :
FourierSeriesApproximation approximation = engine.decompose(); FourierSeries fourier = approximation.getFourier();
Example
The following code snippet shows how to decompose a square function of period 2 to the 10th order :
// function definition UnivariateFunction function = new UnivariateFunction() { public double value(final double x) { final double local = x- FastMath.floor(x / 2) * 2; if (local >= 1) { return-1; } else { return 1; } } }; // Engine parameters UnivariateIntegrator integrator = new LegendreGaussIntegrator(5, 1e-14, 1e-14); FourierDecompositionEngine engine = new FourierDecompositionEngine(integrator); engine.setMaxEvals(Integer.MAX_VALUE); // decompose double period = 2; engine.setOrder(10); engine.setFunction(function, period); FourierSeriesApproximation result = engine.decompose(); FourierSeries fourier = result.getFourier();
The function and its approximation are shown hereunder. The parasite undulations are known as the Gibbs Phenomenon :
Contents
Interfaces
Interface | Summary | Javadoc |
---|---|---|
UnivariateDifferentiableFunction | Extension of UnivariateFunction representing a differentiable univariate function. | ... |
IntegrableUnivariateFunction | Extension of UnivariateRealFunction representing an integrable univariate function. | ... |
UnivariateFunction | Interface representing an univariate function. | ... |
IFastFourierTransformer | Interface representing a FFT. | ... |
Classes
Class | Summary | Javadoc | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
FourierDecompositionEngine | Decompose a UnivariateFunction as a Fourier Series using TrigonometricPolynomialFunction representation. | ... | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
FourierSeries | This class represents a finite Fourier Series. | ... | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
FourierSeriesApproximation | Holder for a UnivariateFunction and its FourierSeries approximation. | ... | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
TrigonometricPolynomialFunction | This class is the Trigonometric Polynomial Function class. | ... | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
TrigonometricPolynomialPrimitive | This class represents a Trigonometric Polynomial Primitive. | ... | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
AbstractFastFourierTransformer | Abstract class representing a FFT. | [http://patrius.cnes.fr/uploads/JavaDocs/V4.7/fr/cnes/sa
IntroductionScopeThe celestial bodies are described by their main features : position and geometry. The positions are ephemeris that must be loaded from models, the geometries are created as one axis ellipsoids. The package provides a factory able to create any celestial body of the solar system. JavadocThe classes for bodies description are available in the package
LinksOrekit bodies : Orekit Bodies architecture description IAU report : Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009 Useful DocumentsNone as of now. Package OverviewThis package overview is from the OREKIT website, under apache license. Features DescriptionCelestial bodiesThe Moon, the Sun and planets of the solar system are all represented by the CelestialBody interface. This class associates a name (eg Sun) to :
Inertially-oriented and body-oriented frames are defined in the following way:
To build a celestial body, the user can:
These methods are detailed in the following sections. Ephemeris LoaderFor any celestial body of the Solar System, the actual computation of its position and velocity relies on the JPL planetary ephemerides files. These files are binary files and loaded thanks to the JPLEphemeridesLoader object. For the moment, this object is only able to load the DE 405 or the DE 406 files which are the most commonly used data files. The former covers the years 1600 to 2200 at maximum precision, the latter covers the years-3000 to +3000 at only slightly reduced precision. The DE 421 file which is the latest JPL ephemeris with fully consistent treatment of planetary and lunar laser ranging data (Folkner et al 2009) is not yet readable by JPLEphemeridesLoader. However the DE 405 file is the basis for the Astronomical Almanac and leads to sufficiently accurate results and, for most purposes, even the accuracy of DE 406 is sufficient. JPLEphemeridesLoaderIn order to generate the ephemerides of one celestial body of the Solar System, one has to use the JPLEphemeridesLoader as follows : final JPLEphemeridesLoader loaderEMB = new JPLEphemeridesLoader(fileName, JPLEphemeridesLoader.EphemerisType.EARTH_MOON); final JPLEphemeridesLoader loaderSSB = new JPLEphemeridesLoader(fileName, JPLEphemeridesLoader.EphemerisType.SOLAR_SYSTEM_BARYCENTER); CelestialBodyFactory.addCelestialBodyLoader(CelestialBodyFactory.EARTH_MOON, loaderEMB); CelestialBodyFactory.addCelestialBodyLoader(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER, loaderSSB); // Reference frame Frame icrf = FramesFactory.getICRF(); // Ephemeris of the Sun JPLEphemeridesLoader loader = new JPLEphemeridesLoader("unxp2000.405", JPLEphemeridesLoader.EphemerisType.SUN); // Creation of the Sun CelestialBody sun = loader.loadCelestialBody(CelestialBodyFactory.SUN); // Coordinates of the Sun given a date and a reference frame PVCoordinates pvSun = sun.getPVCoordinates(date, icrf); When the user wants to create a JPLEphemeridesLoader, first of all, he must supply the folder where the DE 405 and DE 406 files are stored. Then he has to give the name of the file which contains the data (DE 405 or DE 406), the type of celestial body and a date (desired central date) as entries of the JPLEphemeridesLoader. The first argument (name of the data file) may be Once the data is loaded, thanks to the JPLEphemeridesLoader, the user can create the celestial body with the method loadCelestialBody() of the JPLEphemeridesLoader class. To do so, the user has to be sure that all required data is loaded. Indeed, apart from the Moon, the Earth and the Earth-Moon barycenter, the creation of a celestial body requires some data on the Earth Moon barycenter and the Solar System barycenter in order to instanciate the frame in which the ephemerides of the celestial body will be defined. Therefore, the user has to create a JPLEphemeridesLoader for the Earth-Moon barycenter and the Solar System barycenter prior to creating the celestial body, otherwise Patrius will arbitrarily load a DE file to generate the corresponding ephemerides that are used afterwards for the generation of the frame. Then the user has to complete the list of the loaders of the CelestialBodyFactory with these two loaders before calling the method loadCelestialBody(). NB : The user has to clean the CelestialBodyFactory memory if he does not want to work with the previously defined celestial bodies. JPL ephemeridesThe JPL ephemerides data is available on the JPL FTP server [R3]. Available ephemerides :
Simplified analytical modelsMeeus ModelThe Meeus Model is a simplified model which gives the position of the Sun and the Moon with respect to the time T expressed in centuries (TT time scale). This model is an analytical model less precise than the DE ephemerides given by JPL. It is adapted for onboard applications. The class implementing the Meeus Model allows three different models computing the position of the Sun with appropriate equations : the standard model (provided by Jean Meeus), the Stela model and an onboard model. The main differences between these model is the computation of the obliquity of the ecliptic : indeed, its value is fixed to 0 for the standard model, given by an expression involving [math]T, T^{2}[/math]for Stela model and [math]T, T^{2}, T^{3}[/math] for the onboard model. Moreover, the onboard model computed the position in J2000 frame, whereas standard and Stela models use respectively EOD and MOD frame. Regarding the precision, one has to expect a maximum difference of 25593km in position for the Sun and a maximum angular difference of 34.7 (wrt DE405 ephemerides). For the Moon, one has to expect a maximum difference of 26 km in position and a maximum angular difference of 15.2 (wrt DE405 ephemerides). As for the performances (in terms of execution time), the Meeus model for the Sun is faster than the DE405 ephemerides. However, we did not come to the same conclusion for the Moon even by decreasing the degree of precision (number of terms taken into account to compute the latitude, longitude and distance which are needed to compute in fine the position of the Moon).
References for the tables :
In order to build such a Sun or Moon, one has to use the object MeeusSun or MeeusMoon which both extend AbstractCelestialBody. Note that it is not possible to build those celestial bodies from the CelestialBodyFactory, for the moment. Basic board Sun modelThe basic board Sun model is a simplified analytical model which gives the direction of the Sun (the normalized position) with respect to time. This model is similar to the Meeus model (the constants of the model are different) and is adapted for onboard applications. The reference inertial frame is the CIRF. User-defined celestial bodiesUser can defined its own celestial body by using
Example: building Ceres. According to [FDY_BODY_Links IAU report], Ceres has the following parameters:
With d being the interval in days from the standard epoch (the standard epoch is JD 2451545.0, i.e. 2000 January 1 12 hours TDB) Ceres gravitational constant is 6.263E10 m^^3^^kg^^-1^^s^^-2^^. If one knows Ceres motion given for instance by a // Gravitational parameter final double gm = 6.263E10; // Pole motion final IAUPole pole = new IAUPole() { @Override public double getPrimeMeridianAngle(final AbsoluteDate date) { // W final double d = date.durationFrom(AbsoluteDate.J2000_EPOCH) / Constants.JULIAN_DAY; return FastMath.toRadians(170.90) + FastMath.toRadians(952.1532)* d; } @Override public Vector3D getPole(final AbsoluteDate date) { // Pole bias: alpha and delta return new Vector3D(FastMath.toRadians(291), FastMath.toRadians(59)); } }; // Build Ceres body final CelestialBody ceres = new UserDefinedCelestialBody("Ceres", pv, gm, pole); Then
Body shapesThe one-axis ellipsoid is a good approximate model for most planet-size and larger natural bodies. It is the equilibrium shape reached by a fluid body under its own gravity field when it rotates. The symmetry axis is the rotation or polar axis. OneAxisEllipsoidThis type is used to represent the shape of a planet. One very useful implementation represents an ellipsoid (OneAxisEllipsoid). It is constructed from an equatorial radius, a flattening coefficient, and a reference frame that will be used to localize Geodetic points on the shape. It could be interesting to obtained the geodetic coordinates of the nadir point of the satellite. To that purpose, one can use the method getIntersectionPoint(Line, Vector3D, Frame, AbsoluteDate) of the object OneAxisEllipsoid. AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 03, 21), TimeComponents.H12, TimeScalesFactory.getUTC()); Frame frame = FramesFactory.getITRF(); // Body shape model OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, frame); // Satellite on any position circ = new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(50.), FastMath.toRadians(0.), FastMath.toRadians(90.), PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu); // Transform satellite position to position/velocity parameters in EME2000 and ITRF200B PVCoordinates pvSatEME2000 = circ.getPVCoordinates(); PVCoordinates pvSatItrf = frame.getTransformTo(FramesFactory.getEME2000(), date).transformPVCoordinates(pvSatEME2000); Vector3D pSatItrf = pvSatItrf.getPosition(); // Nadir point of the satellite Vector3D pointItrf = new Vector3D.ZERO; Vector3D direction = Vector3D(1., pSatItrf,-1., pointItrf); Line line = new Line(pSatItrf, direction); // intersection point between the ellipsoid and the line that joins the satellite and the center of the body GeodeticPoint nadir = earth.getIntersectionPoint(line, pSatItrf, frame, date); GeometricBodyShape and ExtendedOneAxisEllipsoidThe GeometricBodyShape interface extends BodyShape for more complex computations. ExtendedOneAxisEllipsoid implements it. See the [MIS_SENSORS_PatriusBodySpheroid specific page] for more details. Facet celestial body
This class:
Mesh providerA mesh is described by a list of facets
Available methodsNoe that most methods return exact results. For example,
This method is recursive and is hence in O(log(n)).
If main direction of the field of view is provided, this method is recursive and is in O(log(n)). Otherwise, it is in O(n) and is much slower.
GeodeticPointThe geodetic point is defined by a latitude, a longitude and an altitude in the frame associated to the body. It could be interesting to know the position of a satellite in terms of geodetic coordinates rather than Cartesian ones and vice versa (the corresponding methods in OneAxisEllipsoid). // equatorial radius of the celestial body double ae = 6378137.0; // flatness of the celestial body double f = 1.0 / 298.257222101; // date AbsoluteDate date = AbsoluteDate.J2000_EPOCH; // reference frame attached to the body Frame frame = FramesFactory.getITRF(); // body shape model (ellipsoid) OneAxisEllipsoid model = new OneAxisEllipsoid(ae, f, frame); // transformation with jacobian matrix : cartesian to geodetic // initial cartesian point that will be transformed Vector3D cp = new Vector3D(4637885.347, 121344.708, 4362452.869); // coresponding geodetic point GeodeticPoint gp = model.transform(cp, frame, date); // transformation with jacobian matrix : geodetic to cartesian // initial geodetic point that will be transformed GeodeticPoint gp2 = new GeodeticPoint(0.852479154923577, 0.0423149994747243, 111.6); // coresponding Cartesian point Vector3D cp2 = model.transform(gp2); It could be also interesting to obtain the jacobian matrix of the transformation. // transformation : cartesian to geodetic double[][] jacobianGeodesicWrtCartesian = new double[3][3]; gp = model.transformAndComputeJacobian(cp, frame, date, jacobianGeodesicWrtCartesian); // transformation : geodetic to cartesian double[][] jacobianCartesianWrtGeodesic = new double[3][3]; Vector3D cp2 = model.transformAndComputeJacobian(gp2, jacobianCartesianWrtGeodesic ); Getting StartedTBD ContentsInterfaces
Classes
irius/patrius/math/transform/AbstractFastFourierTransformer.html ...] | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
FastFourierTransformer | Computes the FFT for real and complex arrays. | ... |