1   /* Copyright 2002-2025 CS GROUP
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.estimation.measurements.gnss;
18  
19  import java.util.SortedSet;
20  
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.util.FastMath;
23  import org.junit.jupiter.api.AfterEach;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.BeforeEach;
26  import org.junit.jupiter.api.Test;
27  import org.orekit.Utils;
28  import org.orekit.attitudes.AttitudeProvider;
29  import org.orekit.attitudes.FrameAlignedProvider;
30  import org.orekit.bodies.CelestialBodyFactory;
31  import org.orekit.bodies.GeodeticPoint;
32  import org.orekit.bodies.OneAxisEllipsoid;
33  import org.orekit.estimation.measurements.EstimatedMeasurementBase;
34  import org.orekit.estimation.measurements.GroundStation;
35  import org.orekit.estimation.measurements.ObservableSatellite;
36  import org.orekit.estimation.measurements.generation.EventBasedScheduler;
37  import org.orekit.estimation.measurements.generation.GatheringSubscriber;
38  import org.orekit.estimation.measurements.generation.Generator;
39  import org.orekit.estimation.measurements.generation.InterSatellitesPhaseBuilder;
40  import org.orekit.estimation.measurements.generation.SignSemantic;
41  import org.orekit.frames.FramesFactory;
42  import org.orekit.frames.TopocentricFrame;
43  import org.orekit.gnss.PredefinedGnssSignal;
44  import org.orekit.gnss.SatelliteSystem;
45  import org.orekit.gnss.attitude.GPSBlockIIA;
46  import org.orekit.gnss.attitude.GPSBlockIIR;
47  import org.orekit.orbits.CartesianOrbit;
48  import org.orekit.orbits.Orbit;
49  import org.orekit.propagation.Propagator;
50  import org.orekit.propagation.SpacecraftState;
51  import org.orekit.propagation.analytical.KeplerianPropagator;
52  import org.orekit.propagation.events.InterSatDirectViewDetector;
53  import org.orekit.propagation.events.handlers.ContinueOnEvent;
54  import org.orekit.time.AbsoluteDate;
55  import org.orekit.time.FixedStepSelector;
56  import org.orekit.time.GNSSDate;
57  import org.orekit.time.TimeScalesFactory;
58  import org.orekit.utils.Constants;
59  import org.orekit.utils.IERSConventions;
60  import org.orekit.utils.TimeStampedPVCoordinates;
61  
62  public class InterSatellitesWindUpTest {
63  
64      @Test
65      public void testYawSteering() {
66          // this test corresponds to a classical yaw steering attitude far from turns
67          // where Sun remains largely below orbital plane during the turn (β is about -18.8°)
68          // in this case, yaw does not evolve a lot, so wind-up changes only about 0.14 cycle
69          doTest(new CartesianOrbit(new TimeStampedPVCoordinates(new GNSSDate(1206, 307052.670, SatelliteSystem.GPS).getDate(),
70                                                                 new Vector3D( 8759594.455119, 12170903.262908, 21973798.932235),
71                                                                 new Vector3D(-2957.570165356,  2478.252315039,  -263.042027935)),
72                                    FramesFactory.getGCRF(),
73                                    Constants.EIGEN5C_EARTH_MU),
74                 new GPSBlockIIA(GPSBlockIIA.getDefaultYawRate(17), GPSBlockIIA.DEFAULT_YAW_BIAS,
75                                 AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
76                                 CelestialBodyFactory.getSun(), FramesFactory.getGCRF()),
77                 SatelliteSystem.GPS, 17,
78                 new GroundStation(new TopocentricFrame(earth,
79                                                        new GeodeticPoint(FastMath.toRadians(55.0 + ( 1.0 + 10.0 / 60.0) / 60.0),
80                                                                          FastMath.toRadians(82.0 + (55.0 + 22.0 / 60.0) / 60.0),
81                                                                          160.0),
82                                                        "Новосибирск")),
83                 -0.082134, 0.060814);
84      }
85  
86      @Test
87      public void testMidnightTurn() {
88          // this test corresponds to a Block II-A midnight turn (prn = 07, satellite G37)
89          // where Sun crosses the orbital plane during the turn (β increases from negative to positive values)
90          // this very special case seems to trigger a possible bug in Block IIA attitude model.
91          // The spacecraft starts its turn at about 0.1272°/s and the Sun changes
92          // side, so the satellites keeps turning for about 70 minutes, completing one turn and an half
93          // instead of only one half of a turn. One turn and an half seems unrealistic.
94          // The wind-up effect changes therefore almost linearly by about 1.3 cycle
95          doTest(new CartesianOrbit(new TimeStampedPVCoordinates(new GNSSDate(1218, 287890.543, SatelliteSystem.GPS).getDate(),
96                                                                 new Vector3D(-17920092.444521, -11889104.443797, -15318905.173501),
97                                                                 new Vector3D(   231.983556337,  -3232.849996931,   2163.378049467)),
98                                    FramesFactory.getGCRF(),
99                                    Constants.EIGEN5C_EARTH_MU),
100                new GPSBlockIIA(GPSBlockIIA.getDefaultYawRate(7), GPSBlockIIA.DEFAULT_YAW_BIAS,
101                                AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
102                                CelestialBodyFactory.getSun(), FramesFactory.getGCRF()),
103                SatelliteSystem.GPS, 7,
104                new GroundStation(new TopocentricFrame(earth,
105                                                       new GeodeticPoint(FastMath.toRadians( -(25.0 + 4.0 / 60.0)),
106                                                                         FastMath.toRadians(-(130.0 + 6.0 / 60.0)),
107                                                                         0.0),
108                                                       "Adamstown")),
109                -0.961925, 0.360695);
110     }
111 
112     @Test
113     public void testNoonTurn() {
114         // this test corresponds to a Block II-R noon turn (prn = 11, satellite G46)
115         // where Sun remains slightly above orbital plane during the turn (β is about +1.5°)
116         // this is a regular turn, corresponding to a half turn, so wind-up effect changes by about 0.6 cycle
117         doTest(new CartesianOrbit(new TimeStampedPVCoordinates(new GNSSDate(1225, 509000.063, SatelliteSystem.GPS).getDate(),
118                                                                new Vector3D( 2297608.196826, 20928500.842189, 16246321.092008),
119                                                                new Vector3D(-2810.598090399,  1819.511241767, -1939.009527296)),
120                                   FramesFactory.getGCRF(),
121                                   Constants.EIGEN5C_EARTH_MU),
122                new GPSBlockIIR(GPSBlockIIR.DEFAULT_YAW_RATE,
123                                AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
124                                CelestialBodyFactory.getSun(), FramesFactory.getGCRF()),
125                SatelliteSystem.GPS, 11,
126                new GroundStation(new TopocentricFrame(earth,
127                                                       new GeodeticPoint(FastMath.toRadians(   19.0 + (49.0 + 20.0 / 60.0) / 60.0),
128                                                                         FastMath.toRadians(-(155.0 + (28.0 + 30.0 / 60.0) / 60.0)),
129                                                                         4205.0),
130                                                       "Mauna Kea")),
131                0.349123, 0.972542);
132     }
133 
134     private void doTest(final Orbit emitterOrbit, final AttitudeProvider emitterAttitudeProvider,
135                         final SatelliteSystem emitterSystem, final int emitterPrn,
136                         final GroundStation receiverStation,
137                         final double expectedMin, final double expectedMax) {
138 
139         // create a fake orbit 10000km above the ground station
140         final GeodeticPoint gpOnGround = receiverStation.getBaseFrame().getPoint();
141         final GeodeticPoint gpAltitude = new GeodeticPoint(gpOnGround.getLatitude(),
142                                                            gpOnGround.getLongitude(),
143                                                            10000000.0);
144         final Vector3D position = earth.transform(gpAltitude);
145         final Vector3D velocity = new Vector3D(FastMath.sqrt(emitterOrbit.getMu() / position.getNorm()),
146                                                position.orthogonal());
147         final CartesianOrbit receiverOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(emitterOrbit.getDate(),
148                                                                                              position, velocity),
149                                                                 emitterOrbit.getFrame(),
150                                                                 emitterOrbit.getMu());
151         final AttitudeProvider receiverAttitudeProvider = new FrameAlignedProvider(emitterOrbit.getFrame());
152 
153         // generate phase measurements
154         Generator                   generator   = new Generator();
155         ObservableSatellite         emitterSat  = generator.addPropagator(new KeplerianPropagator(emitterOrbit,
156                                                                                                   emitterAttitudeProvider));
157         ObservableSatellite         receiverSat = generator.addPropagator(new KeplerianPropagator(receiverOrbit,
158                                                                                                   receiverAttitudeProvider));
159         InterSatellitesPhaseBuilder builder     = new InterSatellitesPhaseBuilder(null,
160                                                                                   receiverSat, emitterSat,
161                                                                                   PredefinedGnssSignal.G01.getWavelength(),
162                                                                                   0.01 * PredefinedGnssSignal.G01.getWavelength(),
163                                                                                   1.0,
164                                                                                   new AmbiguityCache());
165         generator.addScheduler(new EventBasedScheduler<>(builder,
166                                                          new FixedStepSelector(60.0, TimeScalesFactory.getUTC()),
167                                                          generator.getPropagator(emitterSat),
168                                                          new InterSatDirectViewDetector(earth, receiverOrbit).
169                                                          withHandler(new ContinueOnEvent()),
170                                                          SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE));
171         final GatheringSubscriber gatherer = new GatheringSubscriber();
172         generator.addSubscriber(gatherer);
173         generator.generate(emitterOrbit.getDate(), emitterOrbit.getDate().shiftedBy(7200));
174         SortedSet<EstimatedMeasurementBase<?>> measurements = gatherer.getGeneratedMeasurements();
175         Assertions.assertEquals(120, measurements.size());
176 
177         InterSatellitesWindUp windUp  = new InterSatellitesWindUpFactory().getWindUp(emitterSystem,  emitterPrn,
178                                                                                      Dipole.CANONICAL_I_J,
179                                                                                      SatelliteSystem.USER_DEFINED_V, 89,
180                                                                                      Dipole.CANONICAL_I_J);
181         Propagator emitterPropagator  = new KeplerianPropagator(emitterOrbit,  emitterAttitudeProvider);
182         Propagator receiverPropagator = new KeplerianPropagator(receiverOrbit, receiverAttitudeProvider);
183         double min = Double.POSITIVE_INFINITY;
184         double max = Double.NEGATIVE_INFINITY;
185         for (EstimatedMeasurementBase<?> m : measurements) {
186             InterSatellitesPhase phase = (InterSatellitesPhase) m.getObservedMeasurement();
187             @SuppressWarnings("unchecked")
188             EstimatedMeasurementBase<InterSatellitesPhase> estimated =
189             (EstimatedMeasurementBase<InterSatellitesPhase>) m.
190                 getObservedMeasurement().
191                 estimateWithoutDerivatives(new SpacecraftState[] {
192                                                receiverPropagator.propagate(phase.getDate()),
193                                                emitterPropagator.propagate(phase.getDate())
194                                            });
195             final double original = estimated.getEstimatedValue()[0];
196             windUp.modifyWithoutDerivatives(estimated);
197             final double modified = estimated.getEstimatedValue()[0];
198             final double correction = modified - original;
199             min = FastMath.min(min, correction);
200             max = FastMath.max(max, correction);
201             Assertions.assertEquals(1,
202                                     estimated.getAppliedEffects().entrySet().stream().
203                                     filter(e -> e.getKey().getEffectName().equals("wind-up")).count());
204         }
205         Assertions.assertEquals(expectedMin, min, 1.0e-5);
206         Assertions.assertEquals(expectedMax, max, 1.0e-5);
207 
208     }
209 
210     @BeforeEach
211     public void setUp() {
212         Utils.setDataRoot("regular-data");
213         earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
214                                      Constants.WGS84_EARTH_FLATTENING,
215                                      FramesFactory.getITRF(IERSConventions.IERS_2010, true));
216     }
217 
218     @AfterEach
219     public void tearDown() {
220         earth = null;
221     }
222 
223     private OneAxisEllipsoid earth;
224 
225 }