1   /* Copyright 2023-2025 Alberto Ferrero
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    * Alberto Ferrero 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.propagation.events;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.util.Binary64;
21  import org.hipparchus.util.Binary64Field;
22  import org.hipparchus.util.FastMath;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.BeforeEach;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.Utils;
27  import org.orekit.bodies.OneAxisEllipsoid;
28  import org.orekit.frames.FramesFactory;
29  import org.orekit.orbits.FieldEquinoctialOrbit;
30  import org.orekit.orbits.FieldOrbit;
31  import org.orekit.propagation.FieldPropagator;
32  import org.orekit.propagation.FieldSpacecraftState;
33  import org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator;
34  import org.orekit.propagation.events.FieldEventsLogger.FieldLoggedEvent;
35  import org.orekit.propagation.events.handlers.FieldContinueOnEvent;
36  import org.orekit.time.FieldAbsoluteDate;
37  import org.orekit.time.TimeScale;
38  import org.orekit.time.TimeScalesFactory;
39  import org.orekit.utils.Constants;
40  import org.orekit.utils.FieldPVCoordinates;
41  import org.orekit.utils.IERSConventions;
42  import org.orekit.utils.PVCoordinates;
43  
44  /** Unit tests for {@link FieldLatitudeRangeCrossingDetector}. */
45  public class FieldLatitudeRangeCrossingDetectorTest {
46  
47      /** Arbitrary Field. */
48      private static final Binary64Field field = Binary64Field.getInstance();
49  
50      @Test
51      public void testRegularCrossing() {
52  
53          final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
54                                                              Constants.WGS84_EARTH_FLATTENING,
55                                                              FramesFactory.getITRF(IERSConventions.IERS_2010, true));
56  
57          FieldLatitudeRangeCrossingDetector<Binary64> d =
58              new FieldLatitudeRangeCrossingDetector<>(v(60.0), v(1.e-6),
59                  earth, FastMath.toRadians(50.0), FastMath.toRadians(60.0)).
60                  withHandler(new FieldContinueOnEvent<>());
61  
62          Assertions.assertEquals(60.0, d.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
63          Assertions.assertEquals(1.0e-6, d.getThreshold().getReal(), 1.0e-15);
64          Assertions.assertEquals(50.0, FastMath.toDegrees(d.getFromLatitude()), 1.0e-14);
65          Assertions.assertEquals(60.0, FastMath.toDegrees(d.getToLatitude()), 1.0e-14);
66          Assertions.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, d.getMaxIterationCount());
67  
68          final TimeScale utc = TimeScalesFactory.getUTC();
69          final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
70          final Vector3D velocity = new Vector3D(505.848, 942.781, 7435.922);
71          final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field, 2003, 9, 16, utc);
72          final FieldOrbit<Binary64> orbit = new FieldEquinoctialOrbit<>(
73                  new FieldPVCoordinates<>(v(1), new PVCoordinates(position,  velocity)),
74                  FramesFactory.getEME2000(), date,
75                  v(Constants.EIGEN5C_EARTH_MU));
76  
77          FieldPropagator<Binary64> propagator =
78              new FieldEcksteinHechlerPropagator<>(orbit,
79                                            Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
80                                            v(Constants.EIGEN5C_EARTH_MU),
81                                            Constants.EIGEN5C_EARTH_C20,
82                                            Constants.EIGEN5C_EARTH_C30,
83                                            Constants.EIGEN5C_EARTH_C40,
84                                            Constants.EIGEN5C_EARTH_C50,
85                                            Constants.EIGEN5C_EARTH_C60);
86  
87          FieldEventsLogger<Binary64> logger = new FieldEventsLogger<>();
88          propagator.addEventDetector(logger.monitorDetector(d));
89  
90          propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));
91          for (FieldLoggedEvent<Binary64> e : logger.getLoggedEvents()) {
92              FieldSpacecraftState<Binary64> state = e.getState();
93              double latitude = earth.transform(state.getPosition(earth.getBodyFrame()),
94                  earth.getBodyFrame(), date).getLatitude().getReal();
95              if (e.isIncreasing()) {
96                  if (state.getPVCoordinates().getVelocity().getZ().getReal() < 0) {
97                      // entering northward
98                      Assertions.assertEquals(60.0, FastMath.toDegrees(latitude), FastMath.toRadians(1e-4));
99                  } else {
100                     // entering southward
101                     Assertions.assertEquals(50.0, FastMath.toDegrees(latitude), FastMath.toRadians(1e-4));
102                 }
103             } else {
104                 if (state.getPVCoordinates().getVelocity().getZ().getReal() < 0) {
105                     // exiting southward
106                     Assertions.assertEquals(50.0, FastMath.toDegrees(latitude), FastMath.toRadians(1e-4));
107                 } else {
108                     // exiting northward
109                     Assertions.assertEquals(60.0, FastMath.toDegrees(latitude), FastMath.toRadians(1e-4));
110                 }
111             }
112         }
113         Assertions.assertEquals(30 * 2, logger.getLoggedEvents().size());
114 
115     }
116 
117     @Test
118     public void testNoCrossing() {
119 
120         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
121                                                             Constants.WGS84_EARTH_FLATTENING,
122                                                             FramesFactory.getITRF(IERSConventions.IERS_2010, true));
123 
124         FieldLatitudeRangeCrossingDetector<Binary64> d =
125                 new FieldLatitudeRangeCrossingDetector<>(v(10.0), v(1.e-6),
126                     earth, FastMath.toRadians(82.0), FastMath.toRadians(87.0)).
127                 withHandler(new FieldContinueOnEvent<>());
128 
129         Assertions.assertEquals(10.0, d.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
130         Assertions.assertEquals(1.0e-6, d.getThreshold().getReal(), 1.0e-15);
131         Assertions.assertEquals(82.0, FastMath.toDegrees(d.getFromLatitude()), 1.0e-14);
132         Assertions.assertEquals(87.0, FastMath.toDegrees(d.getToLatitude()), 1.0e-14);
133         Assertions.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, d.getMaxIterationCount());
134 
135         final TimeScale utc = TimeScalesFactory.getUTC();
136         final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
137         final Vector3D velocity = new Vector3D(505.848, 942.781, 7435.922);
138         final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field, 2003, 9, 16, utc);
139         final FieldOrbit<Binary64> orbit = new FieldEquinoctialOrbit<>(
140                 new FieldPVCoordinates<>(v(1), new PVCoordinates(position,  velocity)),
141                 FramesFactory.getEME2000(), date,
142                 v(Constants.EIGEN5C_EARTH_MU));
143 
144         FieldPropagator<Binary64> propagator =
145             new FieldEcksteinHechlerPropagator<>(orbit,
146                                           Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
147                                           v(Constants.EIGEN5C_EARTH_MU),
148                                           Constants.EIGEN5C_EARTH_C20,
149                                           Constants.EIGEN5C_EARTH_C30,
150                                           Constants.EIGEN5C_EARTH_C40,
151                                           Constants.EIGEN5C_EARTH_C50,
152                                           Constants.EIGEN5C_EARTH_C60);
153 
154         FieldEventsLogger<Binary64> logger = new FieldEventsLogger<Binary64>();
155         propagator.addEventDetector(logger.monitorDetector(d));
156 
157         propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));
158         Assertions.assertEquals(0, logger.getLoggedEvents().size());
159 
160     }
161 
162     /**
163      * Convert double to field value.
164      *
165      * @param value to box.
166      * @return boxed value.
167      */
168     private static Binary64 v(double value) {
169         return new Binary64(value);
170     }
171 
172     @BeforeEach
173     public void setUp() {
174         Utils.setDataRoot("regular-data");
175     }
176 
177 }
178