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  
18  package org.orekit.estimation.iod;
19  
20  import org.hipparchus.geometry.euclidean.threed.Vector3D;
21  import org.hipparchus.util.FastMath;
22  import org.junit.jupiter.api.BeforeEach;
23  import org.orekit.Utils;
24  import org.orekit.bodies.GeodeticPoint;
25  import org.orekit.bodies.OneAxisEllipsoid;
26  import org.orekit.estimation.measurements.AngularAzEl;
27  import org.orekit.estimation.measurements.AngularRaDec;
28  import org.orekit.estimation.measurements.EstimatedMeasurementBase;
29  import org.orekit.estimation.measurements.GroundStation;
30  import org.orekit.estimation.measurements.ObservableSatellite;
31  import org.orekit.frames.Frame;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.frames.TopocentricFrame;
34  import org.orekit.orbits.Orbit;
35  import org.orekit.propagation.Propagator;
36  import org.orekit.propagation.SpacecraftState;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.utils.AbsolutePVCoordinates;
39  import org.orekit.utils.Constants;
40  import org.orekit.utils.IERSConventions;
41  import org.orekit.utils.PVCoordinatesProvider;
42  import org.orekit.utils.TimeStampedPVCoordinates;
43  
44  public abstract class AbstractIodTest {
45      /**
46       * gcrf frame.
47       */
48      protected Frame gcrf;
49  
50      /**
51       * EME2OOO frame.
52       */
53      protected Frame eme2000;
54  
55      /**
56       * ground station for the observations.
57       */
58      protected GroundStation observer;
59  
60      /**
61       * gravitational constant.
62       */
63      protected double mu;
64  
65      /**
66       * itrf frame.
67       */
68      protected Frame itrf;
69  
70      @BeforeEach
71      public void setUp() {
72          Utils.setDataRoot("regular-data");
73  
74          this.mu      = Constants.WGS84_EARTH_MU;
75          this.gcrf    = FramesFactory.getGCRF();
76          this.eme2000 = FramesFactory.getEME2000();
77  
78          this.itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
79          // The ground station is set to Austin, Texas, U.S.A
80          final OneAxisEllipsoid body = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
81                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
82          this.observer = new GroundStation(new TopocentricFrame(body,
83                                                                 new GeodeticPoint(FastMath.toRadians(40),
84                                                                                   FastMath.toRadians(-110),
85                                                                                   2000.0), ""));
86          this.observer.getPrimeMeridianOffsetDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
87          this.observer.getPolarOffsetXDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
88          this.observer.getPolarOffsetYDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
89  
90      }
91  
92      // Computation of LOS angles
93      protected Vector3D getLOSAngles(final Propagator prop, final AbsoluteDate date) {
94          final AngularRaDec raDec = new AngularRaDec(observer, gcrf, date, new double[] { 0.0, 0.0 },
95                                                      new double[] { 1.0, 1.0 },
96                                                      new double[] { 1.0, 1.0 }, new ObservableSatellite(0));
97          return (getEstimatedLineOfSight(raDec, prop, date, gcrf));
98      }
99  
100     protected AngularAzEl getAzEl(final Propagator prop, final AbsoluteDate date) {
101         ObservableSatellite satellite = new ObservableSatellite(0);
102         final AngularAzEl azEl = new AngularAzEl(observer, date, new double[] { 0.0, 0.0 },
103                                                  new double[] { 1.0, 1.0 }, new double[] { 1.0, 1.0 },
104                                                  satellite);
105         EstimatedMeasurementBase<AngularAzEl> estimated = azEl.estimateWithoutDerivatives(new SpacecraftState[] {prop.propagate(date)});
106         return new AngularAzEl(observer, date, estimated.getEstimatedValue(), azEl.getBaseWeight(),
107                                azEl.getTheoreticalStandardDeviation(), satellite);
108     }
109 
110     protected double getRelativeRangeError(final Orbit estimatedGauss, final Orbit orbitRef) {
111 
112         return FastMath.abs(estimatedGauss.getPVCoordinates().getPosition().getNorm() -
113                                     orbitRef.getPVCoordinates().getPosition().getNorm()) /
114                 FastMath.abs(orbitRef.getPVCoordinates().getPosition().getNorm());
115     }
116 
117     // Computation of the relative error in velocity
118     protected double getRelativeVelocityError(final Orbit estimatedGauss, final Orbit orbitRef) {
119 
120         return FastMath.abs(estimatedGauss.getPVCoordinates().getVelocity().getNorm() -
121                                     orbitRef.getPVCoordinates().getVelocity().getNorm()) /
122                 FastMath.abs(orbitRef.getPVCoordinates().getVelocity().getNorm());
123     }
124     
125     /** Calculate the estimated Line Of Sight of a RADEC measurement at a given date.
126      *
127      * @param pvProvider provider for satellite coordinates
128      * @param date the date for which the line of sight must be computed
129      * @param outputFrame output frame for the line of sight
130      * @return the estimate line of Sight of the measurement at the given date.
131      */
132     private Vector3D getEstimatedLineOfSight(final AngularRaDec raDec, final PVCoordinatesProvider pvProvider, final AbsoluteDate date, final Frame outputFrame) {
133         final TimeStampedPVCoordinates satPV       = pvProvider.getPVCoordinates(date, outputFrame);
134         final AbsolutePVCoordinates    satPVInGCRF = new AbsolutePVCoordinates(outputFrame, satPV);
135         final SpacecraftState[]        satState    = new SpacecraftState[] { new SpacecraftState(satPVInGCRF) };
136         final double[]                 angular     = raDec.estimateWithoutDerivatives(satState).getEstimatedValue();
137 
138         // Rotate LOS from RADEC reference frame to output frame
139         return raDec.getReferenceFrame().getStaticTransformTo(outputFrame, date)
140                         .transformVector(new Vector3D(angular[0], angular[1]));
141     }
142 }