ECOM2.java

  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.forces.radiation;

  19. import java.util.ArrayList;
  20. import java.util.Collections;
  21. import java.util.List;

  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.ExtendedPositionProvider;
  35. import org.orekit.utils.ParameterDriver;

  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(ExtendedPositionProvider, 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 {

  78.     /** Parameter name for ECOM model coefficients enabling Jacobian processing. */
  79.     public static final String ECOM_COEFFICIENT = "ECOM coefficient";

  81.     /** Minimum value for ECOM2 estimated parameters. */
  82.     private static final double MIN_VALUE = Double.NEGATIVE_INFINITY;

  84.     /** Maximum value for ECOM2 estimated parameters. */
  85.     private static final double MAX_VALUE = Double.POSITIVE_INFINITY;

  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);

  95.     /** Highest order for parameter along eD axis (satellite --> sun direction). */
  96.     private final int nD;

  98.     /** Highest order for parameter along eB axis. */
  99.     private final int nB;

  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;

  109.     /** Sun model. */
  110.     private final ExtendedPositionProvider sun;

  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 ExtendedPositionProvider sun, final double equatorialRadius) {
  123.         super(sun, new OneAxisEllipsoid(equatorialRadius, 0.0, FramesFactory.getGCRF()),
  124.                 getDefaultEclipseDetectionSettings());
  125.         this.nB = nB;
  126.         this.nD = nD;
  127.         this.coefficients = new ArrayList<>(2 * (nD + nB) + 3);

  129.         // Add parameter along eB axis in alphabetical order
  130.         coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " B0", value, SCALE, MIN_VALUE, MAX_VALUE));
  131.         for (int i = 1; i < nB + 1; i++) {
  132.             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Bcos" + Integer.toString(i - 1), value, SCALE, MIN_VALUE, MAX_VALUE));
  133.         }
  134.         for (int i = nB + 1; i < 2 * nB + 1; i++) {
  135.             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Bsin" + Integer.toString(i - (nB + 1)), value, SCALE, MIN_VALUE, MAX_VALUE));
  136.         }
  137.         // Add driver along eD axis in alphabetical order
  138.         coefficients.add(2 * nB + 1, new ParameterDriver(ECOM_COEFFICIENT + " D0", value, SCALE, MIN_VALUE, MAX_VALUE));
  139.         for (int i = 2 * nB + 2; i < 2 * nB + 2 + nD; i++) {
  140.             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Dcos" + Integer.toString(i - (2 * nB + 2)), value, SCALE, MIN_VALUE, MAX_VALUE));
  141.         }
  142.         for (int i = 2 * nB + 2 + nD; i < 2 * (nB + nD) + 2; i++) {
  143.             coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Dsin" + Integer.toString(i - (2 * nB + nD + 2)), value, SCALE, MIN_VALUE, MAX_VALUE));
  144.         }
  145.         // Add  Y0
  146.         coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Y0", value, SCALE, MIN_VALUE, MAX_VALUE));

  148.         // For ECOM2 model, all parameters are estimated
  149.         coefficients.forEach(parameter -> parameter.setSelected(true));
  150.         this.sun = sun;
  151.     }

  153.     /** {@inheritDoc} */
  154.     @Override
  155.     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {

  157.         // Spacecraft and Sun position vectors (expressed in the spacecraft's frame)
  158.         final Vector3D satPos = s.getPosition();
  159.         final Vector3D sunPos = sun.getPosition(s.getDate(), s.getFrame());

  161.         // Build the coordinate system
  162.         final Vector3D Z = s.getPVCoordinates().getMomentum();
  163.         final Vector3D Y = Z.crossProduct(sunPos).normalize();
  164.         final Vector3D X = Y.crossProduct(Z).normalize();

  166.         // Build eD, eY, eB vectors
  167.         final Vector3D eD = sunPos.subtract(satPos).normalize();
  168.         final Vector3D eY = eD.crossProduct(satPos).normalize();
  169.         final Vector3D eB = eD.crossProduct(eY);

  171.         // Angular argument difference u_s - u
  172.         final double delta_u = FastMath.atan2(satPos.dotProduct(Y), satPos.dotProduct(X));

  174.         // Compute B(u)
  175.         double b_u = parameters[0];
  176.         for (int i = 1; i < nB + 1; i++) {
  177.             final SinCos sc = FastMath.sinCos((2 * i - 1) * delta_u);
  178.             b_u += parameters[i] * sc.cos() + parameters[i + nB] * sc.sin();
  179.         }
  180.         // Compute D(u)
  181.         double d_u = parameters[2 * nB + 1];
  182.         for (int i = 1; i < nD + 1; i++) {
  183.             final SinCos sc = FastMath.sinCos(2 * i * delta_u);
  184.             d_u += parameters[2 * nB + 1 + i] * sc.cos() + parameters[2 * nB + 1 + i + nD] * sc.sin();
  185.         }
  186.         // Return acceleration
  187.         return new Vector3D(d_u, eD, parameters[2 * (nD + nB) + 2], eY, b_u, eB);
  188.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s, final T[] parameters) {

  194.         // Spacecraft and Sun position vectors (expressed in the spacecraft's frame)
  195.         final FieldVector3D<T> satPos = s.getPosition();
  196.         final FieldVector3D<T> sunPos = sun.getPosition(s.getDate(), s.getFrame());

  198.         // Build the coordinate system
  199.         final FieldVector3D<T> Z = s.getPVCoordinates().getMomentum();
  200.         final FieldVector3D<T> Y = Z.crossProduct(sunPos).normalize();
  201.         final FieldVector3D<T> X = Y.crossProduct(Z).normalize();

  203.         // Build eD, eY, eB vectors
  204.         final FieldVector3D<T> eD = sunPos.subtract(satPos).normalize();
  205.         final FieldVector3D<T> eY = eD.crossProduct(satPos).normalize();
  206.         final FieldVector3D<T> eB = eD.crossProduct(eY);

  208.         // Angular argument difference u_s - u
  209.         final T  delta_u = FastMath.atan2(satPos.dotProduct(Y), satPos.dotProduct(X));

  211.         // Compute B(u)
  212.         T b_u =  parameters[0];
  213.         for (int i = 1; i < nB + 1; i++) {
  214.             final FieldSinCos<T> sc = FastMath.sinCos(delta_u.multiply(2 * i - 1));
  215.             b_u = b_u.add(sc.cos().multiply(parameters[i])).add(sc.sin().multiply(parameters[i + nB]));
  216.         }
  217.         // Compute D(u)
  218.         T d_u = parameters[2 * nB + 1];

  220.         for (int i = 1; i < nD + 1; i++) {
  221.             final FieldSinCos<T> sc = FastMath.sinCos(delta_u.multiply(2 * i));
  222.             d_u =  d_u.add(sc.cos().multiply(parameters[2 * nB + 1 + i])).add(sc.sin().multiply(parameters[2 * nB + 1 + i + nD]));
  223.         }
  224.         // Return the acceleration
  225.         return new FieldVector3D<>(d_u, eD, parameters[2 * (nD + nB) + 2], eY, b_u, eB);
  226.     }

  228.     /** {@inheritDoc} */
  229.     @Override
  230.     public List<ParameterDriver> getParametersDrivers() {
  231.         return Collections.unmodifiableList(coefficients);
  232.     }

  234. }
Created with JaCoCo 0.8.12.202403310830