1   /* Copyright 2002-2024 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.forces.radiation;
18  
19  import java.util.ArrayList;
20  import java.util.Collections;
21  import java.util.List;
22  
23  import org.hipparchus.CalculusFieldElement;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.util.FastMath;
27  import org.hipparchus.util.FieldSinCos;
28  import org.hipparchus.util.SinCos;
29  import org.orekit.annotation.DefaultDataContext;
30  import org.orekit.bodies.OneAxisEllipsoid;
31  import org.orekit.frames.FramesFactory;
32  import org.orekit.propagation.FieldSpacecraftState;
33  import org.orekit.propagation.SpacecraftState;
34  import org.orekit.utils.ExtendedPVCoordinatesProvider;
35  import org.orekit.utils.ParameterDriver;
36  
37  /**
38   * The Empirical CODE Orbit Model 2 (ECOM2) of the Center for Orbit Determination in Europe (CODE).
39   * <p>
40   * The drag acceleration is computed as follows :
41   * γ = γ<sub>0</sub> + D(u)e<sub>D</sub> + Y(u)e<sub>Y</sub> + B(u)e<sub>B</sub>
42   * </p> <p>
43   * In the above equation, γ<sub>0</sub> is a selectable a priori model. Since 2013, no
44   * a priori model is used for CODE IGS contribution (i.e. γ<sub>0</sub> = 0). Moreover,
45   * u denotes the satellite's argument of latitude.
46   * </p> <p>
47   * D(u), Y(u) and B(u) are three functions of the ECOM2 model that can be represented
48   * as Fourier series. The coefficients of the Fourier series are estimated during the
49   * estimation process. he ECOM2 model has user-defines upper limits <i>nD</i> and
50   * <i>nB</i> for the Fourier series (i.e. <i>nD</i> for D(u) and <i>nB</i> for
51   * B(u). Y(u) is defined as a constant value).
52   * </p> <p>
53   * It exists several configurations to initialize <i>nD</i> and <i>nB</i> values. However,
54   * Arnold et al recommend to use <b>D2B1</b> (i.e. <i>nD</i> = 1 and <i>nB</i> = 1) and
55   * <b>D4B1</b> (i.e. <i>nD</i> = 2 an <i>nB</i> = 1) configurations. At the opposite, in Arnold paper, it
56   * is recommend to not use <b>D2B0</b> (i.e. <i>nD</i> = 1 and <i>nB</i> = 0) configuration.
57   * </p> <p>
58   * Since Orekit 11.0, it is possible to take into account
59   * the eclipses generated by Moon in the solar radiation
60   * pressure force model using the
61   * {@link #addOccultingBody(ExtendedPVCoordinatesProvider, double)}
62   * method.<br>
63   * <code> ECOM2 srp =</code>
64   * <code>       new ECOM2(1, 1, 0.0, CelestialBodyFactory.getSun(), Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS);</code><br>
65   * <code> srp.addOccultingBody(CelestialBodyFactory.getMoon(), Constants.MOON_EQUATORIAL_RADIUS);</code><br>
66   *
67   * @see "Arnold, Daniel, et al, CODE’s new solar radiation pressure model for GNSS orbit determination,
68   *       Journal of geodesy 89.8 (2015): 775-791."
69   *
70   * @see "Tzu-Pang tseng and Michael Moore, Impact of solar radiation pressure mis-modeling on
71   *       GNSS satellite orbit determination, IGS Worshop, Wuhan, China, 2018."
72   *
73   * @author David Soulard
74   * @since 10.2
75   */
76  public class ECOM2 extends AbstractRadiationForceModel {
77  
78      /** Parameter name for ECOM model coefficients enabling Jacobian processing. */
79      public static final String ECOM_COEFFICIENT = "ECOM coefficient";
80  
81      /** Minimum value for ECOM2 estimated parameters. */
82      private static final double MIN_VALUE = Double.NEGATIVE_INFINITY;
83  
84      /** Maximum value for ECOM2 estimated parameters. */
85      private static final double MAX_VALUE = Double.POSITIVE_INFINITY;
86  
87      /** Parameters scaling factor.
88       * <p>
89       * We use a power of 2 to avoid numeric noise introduction
90       * in the multiplications/divisions sequences.
91       * </p>
92       */
93      private final double SCALE = FastMath.scalb(1.0, -22);
94  
95      /** Highest order for parameter along eD axis (satellite --> sun direction). */
96      private final int nD;
97  
98      /** Highest order for parameter along eB axis. */
99      private final int nB;
100 
101     /** Estimated acceleration coefficients.
102      * <p>
103      * The 2 * nD first driver are Fourier driver along eD, axis,
104      * then along eY, then 2*nB following are along eB axis.
105      * </p>
106      */
107     private final List<ParameterDriver> coefficients;
108 
109     /** Sun model. */
110     private final ExtendedPVCoordinatesProvider sun;
111 
112     /**
113      * Constructor.
114      * @param nD truncation rank of Fourier series in D term.
115      * @param nB truncation rank of Fourier series in B term.
116      * @param value parameters initial value
117      * @param sun provide for Sun parameter
118      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
119      */
120     @DefaultDataContext
121     public ECOM2(final int nD, final int nB, final double value,
122                  final ExtendedPVCoordinatesProvider sun, final double equatorialRadius) {
123         super(sun, new OneAxisEllipsoid(equatorialRadius, 0.0, FramesFactory.getGCRF()));
124         this.nB = nB;
125         this.nD = nD;
126         this.coefficients = new ArrayList<>(2 * (nD + nB) + 3);
127 
128         // Add parameter along eB axis in alphabetical order
129         coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " B0", value, SCALE, MIN_VALUE, MAX_VALUE));
130         for (int i = 1; i < nB + 1; i++) {
131             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Bcos" + Integer.toString(i - 1), value, SCALE, MIN_VALUE, MAX_VALUE));
132         }
133         for (int i = nB + 1; i < 2 * nB + 1; i++) {
134             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Bsin" + Integer.toString(i - (nB + 1)), value, SCALE, MIN_VALUE, MAX_VALUE));
135         }
136         // Add driver along eD axis in alphabetical order
137         coefficients.add(2 * nB + 1, new ParameterDriver(ECOM_COEFFICIENT + " D0", value, SCALE, MIN_VALUE, MAX_VALUE));
138         for (int i = 2 * nB + 2; i < 2 * nB + 2 + nD; i++) {
139             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Dcos" + Integer.toString(i - (2 * nB + 2)), value, SCALE, MIN_VALUE, MAX_VALUE));
140         }
141         for (int i = 2 * nB + 2 + nD; i < 2 * (nB + nD) + 2; i++) {
142             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Dsin" + Integer.toString(i - (2 * nB + nD + 2)), value, SCALE, MIN_VALUE, MAX_VALUE));
143         }
144         // Add  Y0
145         coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Y0", value, SCALE, MIN_VALUE, MAX_VALUE));
146 
147         // For ECOM2 model, all parameters are estimated
148         coefficients.forEach(parameter -> parameter.setSelected(true));
149         this.sun = sun;
150     }
151 
152     /** {@inheritDoc} */
153     @Override
154     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {
155 
156         // Spacecraft and Sun position vectors (expressed in the spacecraft's frame)
157         final Vector3D satPos = s.getPosition();
158         final Vector3D sunPos = sun.getPosition(s.getDate(), s.getFrame());
159 
160         // Build the coordinate system
161         final Vector3D Z = s.getPVCoordinates().getMomentum();
162         final Vector3D Y = Z.crossProduct(sunPos).normalize();
163         final Vector3D X = Y.crossProduct(Z).normalize();
164 
165         // Build eD, eY, eB vectors
166         final Vector3D eD = sunPos.subtract(satPos).normalize();
167         final Vector3D eY = eD.crossProduct(satPos).normalize();
168         final Vector3D eB = eD.crossProduct(eY);
169 
170         // Angular argument difference u_s - u
171         final double delta_u = FastMath.atan2(satPos.dotProduct(Y), satPos.dotProduct(X));
172 
173         // Compute B(u)
174         double b_u = parameters[0];
175         for (int i = 1; i < nB + 1; i++) {
176             final SinCos sc = FastMath.sinCos((2 * i - 1) * delta_u);
177             b_u += parameters[i] * sc.cos() + parameters[i + nB] * sc.sin();
178         }
179         // Compute D(u)
180         double d_u = parameters[2 * nB + 1];
181         for (int i = 1; i < nD + 1; i++) {
182             final SinCos sc = FastMath.sinCos(2 * i * delta_u);
183             d_u += parameters[2 * nB + 1 + i] * sc.cos() + parameters[2 * nB + 1 + i + nD] * sc.sin();
184         }
185         // Return acceleration
186         return new Vector3D(d_u, eD, parameters[2 * (nD + nB) + 2], eY, b_u, eB);
187     }
188 
189     /** {@inheritDoc} */
190     @Override
191     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s, final T[] parameters) {
192 
193         // Spacecraft and Sun position vectors (expressed in the spacecraft's frame)
194         final FieldVector3D<T> satPos = s.getPosition();
195         final FieldVector3D<T> sunPos = sun.getPosition(s.getDate(), s.getFrame());
196 
197         // Build the coordinate system
198         final FieldVector3D<T> Z = s.getPVCoordinates().getMomentum();
199         final FieldVector3D<T> Y = Z.crossProduct(sunPos).normalize();
200         final FieldVector3D<T> X = Y.crossProduct(Z).normalize();
201 
202         // Build eD, eY, eB vectors
203         final FieldVector3D<T> eD = sunPos.subtract(satPos).normalize();
204         final FieldVector3D<T> eY = eD.crossProduct(satPos).normalize();
205         final FieldVector3D<T> eB = eD.crossProduct(eY);
206 
207         // Angular argument difference u_s - u
208         final T  delta_u = FastMath.atan2(satPos.dotProduct(Y), satPos.dotProduct(X));
209 
210         // Compute B(u)
211         T b_u =  parameters[0];
212         for (int i = 1; i < nB + 1; i++) {
213             final FieldSinCos<T> sc = FastMath.sinCos(delta_u.multiply(2 * i - 1));
214             b_u = b_u.add(sc.cos().multiply(parameters[i])).add(sc.sin().multiply(parameters[i + nB]));
215         }
216         // Compute D(u)
217         T d_u = parameters[2 * nB + 1];
218 
219         for (int i = 1; i < nD + 1; i++) {
220             final FieldSinCos<T> sc = FastMath.sinCos(delta_u.multiply(2 * i));
221             d_u =  d_u.add(sc.cos().multiply(parameters[2 * nB + 1 + i])).add(sc.sin().multiply(parameters[2 * nB + 1 + i + nD]));
222         }
223         // Return the acceleration
224         return new FieldVector3D<>(d_u, eD, parameters[2 * (nD + nB) + 2], eY, b_u, eB);
225     }
226 
227     /** {@inheritDoc} */
228     @Override
229     public List<ParameterDriver> getParametersDrivers() {
230         return Collections.unmodifiableList(coefficients);
231     }
232 
233 }
234