1   /* Copyright 2022-2025 Luc Maisonobe
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.bodies;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.CalculusFieldElement;
23  import org.hipparchus.Field;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.geometry.euclidean.twod.FieldVector2D;
27  import org.hipparchus.util.Binary64Field;
28  import org.hipparchus.util.FastMath;
29  import org.junit.jupiter.api.Assertions;
30  import org.junit.jupiter.api.BeforeEach;
31  import org.junit.jupiter.api.Test;
32  import org.orekit.Utils;
33  import org.orekit.frames.FramesFactory;
34  import org.orekit.time.FieldAbsoluteDate;
35  import org.orekit.time.FieldTimeInterpolator;
36  import org.orekit.utils.CartesianDerivativesFilter;
37  import org.orekit.utils.Constants;
38  import org.orekit.utils.IERSConventions;
39  import org.orekit.utils.TimeStampedFieldPVCoordinates;
40  import org.orekit.utils.TimeStampedFieldPVCoordinatesHermiteInterpolator;
41  
42  public class FieldEllipseTest {
43  
44      @Test
45      void testMeridianShape() {
46          doTestMeridianShape(Binary64Field.getInstance());
47      }
48  
49      @Test
50      void testEquatorialShape() {
51          doTestEquatorialShape(Binary64Field.getInstance());
52      }
53  
54      @Test
55      void testProjectionDerivatives() {
56          doTestProjectionDerivatives(Binary64Field.getInstance());
57      }
58  
59      @Test
60      void testMinRadiusOfCurvature() {
61          doTestMinRadiusOfCurvature(Binary64Field.getInstance());
62      }
63  
64      @Test
65      void testMaxRadiusOfCurvature() {
66          doTestMaxRadiusOfCurvature(Binary64Field.getInstance());
67      }
68  
69      @Test
70      void testFlatEllipse() {
71          doTestFlatEllipse(Binary64Field.getInstance());
72      }
73  
74      private <T extends CalculusFieldElement<T>> void doTestMeridianShape(final Field<T> field) {
75          OneAxisEllipsoid model =
76                  new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
77                                       Constants.WGS84_EARTH_FLATTENING,
78                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
79          FieldEllipse<T> e = model.getPlaneSection(new FieldVector3D<>(field, new Vector3D(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0)),
80                                                    FieldVector3D.getPlusJ(field));
81          Assertions.assertEquals(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
82                                  e.getA().getReal(),
83                                  1.0e-15 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS);
84          Assertions.assertEquals(Constants.WGS84_EARTH_EQUATORIAL_RADIUS * (1 - Constants.WGS84_EARTH_FLATTENING),
85                                  e.getB().getReal(),
86                                  1.0e-15 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS);
87          Assertions.assertEquals(0.5 * FastMath.PI, FieldVector3D.angle(Vector3D.PLUS_J, e.getU()).getReal(), 1.0e-15);
88          Assertions.assertEquals(0.5 * FastMath.PI, FieldVector3D.angle(Vector3D.PLUS_K, e.getU()).getReal(), 1.0e-15);
89          Assertions.assertEquals(0.5 * FastMath.PI, FieldVector3D.angle(Vector3D.PLUS_I, e.getV()).getReal(), 1.0e-15);
90          Assertions.assertEquals(0.5 * FastMath.PI, FieldVector3D.angle(Vector3D.PLUS_J, e.getV()).getReal(), 1.0e-15);
91      }
92  
93      private <T extends CalculusFieldElement<T>> void doTestEquatorialShape(final Field<T> field) {
94          OneAxisEllipsoid model =
95                  new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
96                                       Constants.WGS84_EARTH_FLATTENING,
97                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
98          FieldEllipse<T> e = model.getPlaneSection(new FieldVector3D<>(field, new Vector3D(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0)),
99                                                    FieldVector3D.getPlusK(field));
100         Assertions.assertEquals(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
101                                 e.getA().getReal(),
102                                 1.0e-15 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS);
103         Assertions.assertEquals(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
104                                 e.getB().getReal(),
105                                 1.0e-15 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS);
106     }
107 
108     private <T extends CalculusFieldElement<T>> void doTestProjectionDerivatives(final Field<T> field) {
109         final T zero = field.getZero();
110         FieldEllipse<T> e = new FieldEllipse<>(FieldVector3D.getZero(field), FieldVector3D.getPlusI(field), FieldVector3D.getPlusJ(field),
111                                                zero.newInstance(6.4e6), zero.newInstance(6.3e6),
112                                                FramesFactory.getGCRF());
113         TimeStampedFieldPVCoordinates<T> linearMotion =
114                 new TimeStampedFieldPVCoordinates<>(FieldAbsoluteDate.getJ2000Epoch(field),
115                                                     new FieldVector3D<>(zero.newInstance(7.0e6), zero.newInstance(5.0e6), zero),
116                                                     new FieldVector3D<>(zero.newInstance(3.0e3), zero.newInstance(4.0e3), zero),
117                                                     FieldVector3D.getZero(field));
118         TimeStampedFieldPVCoordinates<T> g0 = e.projectToEllipse(linearMotion);
119         List<TimeStampedFieldPVCoordinates<T>> sample = new ArrayList<>();
120         for (double dt = -0.25; dt <= 0.25; dt += 0.125) {
121             sample.add(e.projectToEllipse(linearMotion.shiftedBy(dt)));
122         }
123 
124         // create interpolator
125         final FieldTimeInterpolator<TimeStampedFieldPVCoordinates<T>, T> interpolator =
126                 new TimeStampedFieldPVCoordinatesHermiteInterpolator<>(sample.size(), CartesianDerivativesFilter.USE_P);
127 
128         TimeStampedFieldPVCoordinates<T> ref = interpolator.interpolate(g0.getDate(), sample);
129         Assertions.assertEquals(0,
130                                 FieldVector3D.distance(g0.getPosition(), ref.getPosition()).divide(ref.getPosition().getNorm()).getReal(),
131                                 1.0e-15);
132         Assertions.assertEquals(0,
133                                 FieldVector3D.distance(g0.getVelocity(), ref.getVelocity()).divide(ref.getVelocity().getNorm()).getReal(),
134                                 6.0e-12);
135         Assertions.assertEquals(0,
136                                 FieldVector3D.distance(g0.getAcceleration(), ref.getAcceleration()).divide(ref.getAcceleration().getNorm()).getReal(),
137                                 8.0e-8);
138 
139     }
140 
141     private <T extends CalculusFieldElement<T>> void doTestMinRadiusOfCurvature(final Field<T> field) {
142         final T zero = field.getZero();
143         final double a = 100.0;
144         final double b =  50.0;
145         FieldEllipse<T> e = new FieldEllipse<>(FieldVector3D.getZero(field), FieldVector3D.getPlusI(field), FieldVector3D.getPlusJ(field),
146                                                zero.newInstance(a), zero.newInstance(b), FramesFactory.getGCRF());
147         FieldVector2D<T> point = new FieldVector2D<>(zero.newInstance(10 * a), zero);
148         Assertions.assertEquals(b * b / a,
149                                 FieldVector2D.distance(e.projectToEllipse(point), e.getCenterOfCurvature(point)).getReal(),
150                                 1.0e-15);
151     }
152 
153     private <T extends CalculusFieldElement<T>> void doTestMaxRadiusOfCurvature(final Field<T> field) {
154         final T zero = field.getZero();
155         final double a = 100.0;
156         final double b =  50.0;
157         FieldEllipse<T> e = new FieldEllipse<>(FieldVector3D.getZero(field), FieldVector3D.getPlusI(field), FieldVector3D.getPlusJ(field),
158                                                zero.newInstance(a), zero.newInstance(b), FramesFactory.getGCRF());
159         FieldVector2D<T> point = new FieldVector2D<>(zero, zero.newInstance(10 * b));
160         Assertions.assertEquals(a * a / b,
161                                 FieldVector2D.distance(e.projectToEllipse(point), e.getCenterOfCurvature(point)).getReal(),
162                                 1.0e-15);
163     }
164 
165     private <T extends CalculusFieldElement<T>> void doTestFlatEllipse(final Field<T> field) {
166         final T zero = field.getZero();
167         final double a = 0.839;
168         final double b = 0.176;
169         FieldEllipse<T> e = new FieldEllipse<>(FieldVector3D.getZero(field), FieldVector3D.getPlusI(field), FieldVector3D.getPlusJ(field),
170                         zero.newInstance(a), zero.newInstance(b), FramesFactory.getGCRF());
171         final FieldVector2D<T> close = e.projectToEllipse(new FieldVector2D<>(zero.newInstance(2.0), zero.newInstance(4.0)));
172         Assertions.assertEquals(1.0,
173                                 close.getX().multiply(close.getX()).divide(a * a).add(close.getY().multiply(close.getY()).divide(b * b)).getReal(),
174                                 1.0e-15);
175     }
176 
177     @BeforeEach
178     public void setUp() {
179         Utils.setDataRoot("regular-data");
180     }
181 
182 }
183