1   /* Copyright 2022-2025 Bryan Cazabonne
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    * Bryan Cazabonne 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.semianalytical.dsst.forces;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.util.MathArrays;
22  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
23  import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
24  
25  /**
26   * Zeis model for J2-squared second-order terms.
27   *
28   * @see "ZEIS, Eric and CEFOLA, P. Computerized algebraic utilities for the
29   *      construction of nonsingular satellite theories. Journal of Guidance and
30   *      Control, 1980, vol. 3, no 1, p. 48-54."
31   *
32   * @see "SAN-JUAN, Juan F., LÓPEZ, Rosario, et CEFOLA, Paul J. A Second-Order
33   *      Closed-Form $$ J_2 $$ Model for the Draper Semi-Analytical Satellite
34   *      Theory. The Journal of the Astronautical Sciences, 2022, p. 1-27."
35   *
36   * @author Bryan Cazabonne
37   * @since 12.0
38   */
39  public class ZeisModel implements J2SquaredModel {
40  
41      /**
42       * Retrograde factor I.
43       * <p>
44       * DSST model needs equinoctial orbit as internal representation. Classical
45       * equinoctial elements have discontinuities when inclination is close to zero.
46       * In this representation, I = +1. <br>
47       * To avoid this discontinuity, another representation exists and equinoctial
48       * elements can be expressed in a different way, called "retrograde" orbit. This
49       * implies I = -1. <br>
50       * As Orekit doesn't implement the retrograde orbit, I is always set to +1. But
51       * for the sake of consistency with the theory, the retrograde factor has been
52       * kept in the formulas.
53       * </p>
54       */
55      private static final int I = 1;
56  
57      /** Constructor. */
58      public ZeisModel() {
59          // Nothing to do...
60      }
61  
62      /** {@inheritDoc}. */
63      @Override
64      public double[] computeMeanEquinoctialSecondOrderTerms(final DSSTJ2SquaredClosedFormContext context) {
65  
66          // Auxiliary elements
67          final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
68  
69          // Zeis constant
70          final double c2z = computeC2Z(context);
71  
72          // Useful terms
73          final double s2mf    = 19.0 * context.getS2() - 15.0;
74          final double s2pIcmo = context.getS2() + I * context.getC() - 1.0;
75          final double s4mts2  = 19.0 * context.getS2() * context.getS2() - 30.0 * context.getS2() + 12.0;
76  
77          // Second-order terms (Ref [2] Eq. 37)
78          final double deltaA =  0.0;
79          final double deltaK = -c2z * auxiliaryElements.getH() * s2mf * s2pIcmo;
80          final double deltaH =  c2z * auxiliaryElements.getK() * s2mf * s2pIcmo;
81          final double deltaQ = -c2z * context.getC() * auxiliaryElements.getP() * s2mf;
82          final double deltaP =  c2z * context.getC() * auxiliaryElements.getQ() * s2mf;
83          final double deltaM =  0.5 * c2z * (2.0 * s2mf * s2pIcmo + 5.0 * s4mts2 * context.getEta());
84  
85          // Return
86          return new double[] { deltaA, deltaK, deltaH, deltaQ, deltaP, deltaM };
87  
88      }
89  
90      /** {@inheritDoc}. */
91      @Override
92      public <T extends CalculusFieldElement<T>> T[] computeMeanEquinoctialSecondOrderTerms(final FieldDSSTJ2SquaredClosedFormContext<T> context) {
93  
94          // Auxiliary elements
95          final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
96  
97          // Field
98          final Field<T> field = auxiliaryElements.getDate().getField();
99  
100         // Zeis constant
101         final T c2z = computeC2Z(context);
102 
103         // Useful terms
104         final T s2mf    = context.getS2().multiply(19.0).subtract(15.0);
105         final T s2pIcmo = context.getS2().add(context.getC().multiply(I)).subtract(1.0);
106         final T s4mts2  = context.getS2().multiply(context.getS2()).multiply(19.0).subtract(context.getS2().multiply(30.0)).add(12.0);
107 
108         // Second-order terms (Ref [2] Eq. 37)
109         final T deltaA = field.getZero();
110         final T deltaK = c2z.multiply(auxiliaryElements.getH()).multiply(s2mf).multiply(s2pIcmo).negate();
111         final T deltaH = c2z.multiply(auxiliaryElements.getK()).multiply(s2mf).multiply(s2pIcmo);
112         final T deltaQ = c2z.multiply(context.getC()).multiply(auxiliaryElements.getP()).multiply(s2mf).negate();
113         final T deltaP = c2z.multiply(context.getC()).multiply(auxiliaryElements.getQ()).multiply(s2mf);
114         final T deltaM = c2z.multiply(0.5).multiply(s2mf.multiply(s2pIcmo).multiply(2.0).add(s4mts2.multiply(context.getEta()).multiply(5.0)));
115 
116         // Return
117         final T[] terms = MathArrays.buildArray(field, 6);
118         terms[0] = deltaA;
119         terms[1] = deltaK;
120         terms[2] = deltaH;
121         terms[3] = deltaQ;
122         terms[4] = deltaP;
123         terms[5] = deltaM;
124         return terms;
125 
126     }
127 
128     /**
129      * Get the value of the Zeis constant.
130      *
131      * @param context model context
132      * @return the value of the Zeis constant
133      */
134     public double computeC2Z(final DSSTJ2SquaredClosedFormContext context) {
135         final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
136         return 0.75 * context.getAlpha4() * auxiliaryElements.getMeanMotion() / (context.getA4() * context.getEta());
137     }
138 
139     /**
140      * Get the value of the Zeis constant.
141      *
142      * @param context model context
143      * @param <T> type of the elements
144      * @return the value of the Zeis constant
145      */
146     public <T extends CalculusFieldElement<T>> T computeC2Z(final FieldDSSTJ2SquaredClosedFormContext<T> context) {
147         final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
148         return auxiliaryElements.getMeanMotion().multiply(context.getAlpha4()).multiply(0.75).divide(context.getA4().multiply(context.getEta()));
149     }
150 
151 }