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.propagation.semianalytical.dsst.utilities;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.util.FastMath;
22  
23  /** Compute the G<sub>ms</sub><sup>j</sup> and the H<sub>ms</sub><sup>j</sup>
24   *  polynomials in the equinoctial elements h, k and the direction cosines α and β
25   *  and their partial derivatives with respect to k, h, α and β.
26   *  <p>
27   *  The expressions used are equations 2.7.5-(1)(2) from the Danielson paper.
28   *  </p>
29   *  @author Romain Di Costanzo
30   *  @author Bryan Cazabonne (field translation)
31   * @param <T> type of the field elements
32   */
33  public class FieldGHmsjPolynomials <T extends CalculusFieldElement<T>> {
34      /** C<sub>j</sub>(k, h), S<sub>j</sub>(k, h) coefficient.
35       * (k, h) are the (x, y) component of the eccentricity vector in equinoctial elements
36       */
37      private final FieldCjSjCoefficient<T> cjsjKH;
38  
39      /** C<sub>j</sub>(α, β), S<sub>j</sub>(α, β) coefficient.
40       * (α, β) are the direction cosines
41       */
42      private final FieldCjSjCoefficient<T> cjsjAB;
43  
44      /** Is the orbit represented as a retrograde orbit.
45       *  I = -1 if yes, +1 otherwise.
46       */
47      private int                   I;
48  
49      /** Create a set of G<sub>ms</sub><sup>j</sup> and H<sub>ms</sub><sup>j</sup> polynomials.
50       *  @param k X component of the eccentricity vector
51       *  @param h Y component of the eccentricity vector
52       *  @param alpha direction cosine α
53       *  @param beta direction cosine β
54       *  @param retroFactor -1 if the orbit is represented as retrograde, +1 otherwise
55       *  @param field field element
56       **/
57      public FieldGHmsjPolynomials(final T k, final T h,
58                              final T alpha, final T beta,
59                              final int retroFactor,
60                              final Field<T> field) {
61          this.cjsjKH = new FieldCjSjCoefficient<>(k, h, field);
62          this.cjsjAB = new FieldCjSjCoefficient<>(alpha, beta, field);
63          this.I = retroFactor;
64      }
65  
66      /** Get the G<sub>ms</sub><sup>j</sup> coefficient.
67       * @param m m subscript
68       * @param s s subscript
69       * @param j order
70       * @return the G<sub>ms</sub><sup>j</sup>
71       */
72      public T getGmsj(final int m, final int s, final int j) {
73          final int sMj = FastMath.abs(s - j);
74          final T gms;
75          if (FastMath.abs(s) <= m) {
76              final int mMis = m - I * s;
77              gms = cjsjKH.getCj(sMj).multiply(cjsjAB.getCj(mMis)).
78                    subtract(cjsjKH.getSj(sMj).multiply(cjsjAB.getSj(mMis)).multiply(sgn(s - j)).multiply(I));
79          } else {
80              final int sMim = FastMath.abs(s - I * m);
81              gms = cjsjKH.getCj(sMj).multiply(cjsjAB.getCj(sMim)).
82                    add(cjsjKH.getSj(sMj).multiply(cjsjAB.getSj(sMim)).multiply(sgn(s - j)).multiply(sgn(s - m)));
83          }
84          return gms;
85      }
86  
87      /** Get the H<sub>ms</sub><sup>j</sup> coefficient.
88       * @param m m subscript
89       * @param s s subscript
90       * @param j order
91       * @return the H<sub>ms</sub><sup>j</sup>
92       */
93      public T getHmsj(final int m, final int s, final int j) {
94          final int sMj = FastMath.abs(s - j);
95          final T hms;
96          if (FastMath.abs(s) <= m) {
97              final int mMis = m - I * s;
98              hms = cjsjKH.getCj(sMj).multiply(cjsjAB.getSj(mMis)).multiply(I).
99                              add(cjsjKH.getSj(sMj).multiply(cjsjAB.getCj(mMis)).multiply(sgn(s - j)));
100         } else {
101             final int sMim = FastMath.abs(s - I * m);
102             hms = cjsjKH.getCj(sMj).multiply(cjsjAB.getSj(sMim)).multiply(-sgn(s - m)).
103                   add(cjsjKH.getSj(sMj).multiply(cjsjAB.getCj(sMim)).multiply(sgn(s - j)));
104         }
105         return hms;
106     }
107 
108     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>k</sub> coefficient.
109      * @param m m subscript
110      * @param s s subscript
111      * @param j order
112      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>k</sub>
113      */
114     public T getdGmsdk(final int m, final int s, final int j) {
115         final int sMj = FastMath.abs(s - j);
116         final T dGmsdk;
117         if (FastMath.abs(s) <= m) {
118             final int mMis = m - I * s;
119             dGmsdk = cjsjKH.getDcjDk(sMj).multiply(cjsjAB.getCj(mMis)).
120                      subtract(cjsjKH.getDsjDk(sMj).multiply(cjsjAB.getSj(mMis)).multiply(I).multiply(sgn(s - j)));
121         } else {
122             final int sMim = FastMath.abs(s - I * m);
123             dGmsdk = cjsjKH.getDcjDk(sMj).multiply(cjsjAB.getCj(sMim)).
124                      add(cjsjKH.getDsjDk(sMj).multiply(cjsjAB.getSj(sMim)).multiply(sgn(s - m)).multiply(sgn(s - j)));
125         }
126         return dGmsdk;
127     }
128 
129     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>h</sub> coefficient.
130      * @param m m subscript
131      * @param s s subscript
132      * @param j order
133      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>h</sub>
134      */
135     public T getdGmsdh(final int m, final int s, final int j) {
136         final int sMj = FastMath.abs(s - j);
137         final T dGmsdh;
138         if (FastMath.abs(s) <= m) {
139             final int mMis = m - I * s;
140             dGmsdh = cjsjKH.getDcjDh(sMj).multiply(cjsjAB.getCj(mMis)).
141                      subtract(cjsjKH.getDsjDh(sMj).multiply(cjsjAB.getSj(mMis)).multiply(I).multiply(sgn(s - j)));
142         } else {
143             final int sMim = FastMath.abs(s - I * m);
144             dGmsdh = cjsjKH.getDcjDh(sMj).multiply(cjsjAB.getCj(sMim)).
145                      add(cjsjKH.getDsjDh(sMj).multiply(cjsjAB.getSj(sMim)).multiply(sgn(s - m)).multiply(sgn(s - j)));
146         }
147         return dGmsdh;
148     }
149 
150     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>α</sub> coefficient.
151      * @param m m subscript
152      * @param s s subscript
153      * @param j order
154      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>α</sub>
155      */
156     public T getdGmsdAlpha(final int m, final int s, final int j) {
157         final int sMj  = FastMath.abs(s - j);
158         final T dGmsdAl;
159         if (FastMath.abs(s) <= m) {
160             final int mMis = m - I * s;
161             dGmsdAl = cjsjKH.getCj(sMj).multiply(cjsjAB.getDcjDk(mMis)).
162                       subtract(cjsjKH.getSj(sMj).multiply(cjsjAB.getDsjDk(mMis)).multiply(I).multiply(sgn(s - j)));
163         } else {
164             final int sMim = FastMath.abs(s - I * m);
165             dGmsdAl = cjsjKH.getCj(sMj).multiply(cjsjAB.getDcjDk(sMim)).
166                       add(cjsjKH.getSj(sMj).multiply(cjsjAB.getDsjDk(sMim)).multiply(sgn(s - j)).multiply(sgn(s - m)));
167         }
168         return dGmsdAl;
169     }
170 
171     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>β</sub> coefficient.
172      * @param m m subscript
173      * @param s s subscript
174      * @param j order
175      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>β</sub>
176      */
177     public T getdGmsdBeta(final int m, final int s, final int j) {
178         final int sMj = FastMath.abs(s - j);
179         final T dGmsdBe;
180         if (FastMath.abs(s) <= m) {
181             final int mMis = m - I * s;
182             dGmsdBe = cjsjKH.getCj(sMj).multiply(cjsjAB.getDcjDh(mMis)).
183                       subtract(cjsjKH.getSj(sMj).multiply(cjsjAB.getDsjDh(mMis)).multiply(I).multiply(sgn(s - j)));
184         } else {
185             final int sMim = FastMath.abs(s - I * m);
186             dGmsdBe = cjsjKH.getCj(sMj).multiply(cjsjAB.getDcjDh(sMim)).
187                       add(cjsjKH.getSj(sMj).multiply(cjsjAB.getDsjDh(sMim)).multiply(sgn(s - j)).multiply(sgn(s - m)));
188         }
189         return dGmsdBe;
190     }
191 
192     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>k</sub> coefficient.
193      * @param m m subscript
194      * @param s s subscript
195      * @param j order
196      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>k</sub>
197      */
198     public T getdHmsdk(final int m, final int s, final int j) {
199         final int sMj = FastMath.abs(s - j);
200         final T dHmsdk;
201         if (FastMath.abs(s) <= m) {
202             final int mMis = m - I * s;
203             dHmsdk = cjsjKH.getDcjDk(sMj).multiply(cjsjAB.getSj(mMis)).multiply(I).
204                      add(cjsjKH.getDsjDk(sMj).multiply(cjsjAB.getCj(mMis)).multiply(sgn(s - j)));
205         } else {
206             final int sMim = FastMath.abs(s - I * m);
207             dHmsdk = cjsjKH.getDcjDk(sMj).multiply(cjsjAB.getSj(sMim)).multiply(-sgn(s - m)).
208                      add(cjsjKH.getDsjDk(sMj).multiply(cjsjAB.getCj(sMim)).multiply(sgn(s - j)));
209         }
210         return dHmsdk;
211     }
212 
213     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>h</sub> coefficient.
214      * @param m m subscript
215      * @param s s subscript
216      * @param j order
217      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>h</sub>
218      */
219     public T getdHmsdh(final int m,  final int s, final int j) {
220         final int sMj = FastMath.abs(s - j);
221         final T dHmsdh;
222         if (FastMath.abs(s) <= m) {
223             final int mMis = m - I * s;
224             dHmsdh = cjsjKH.getDcjDh(sMj).multiply(cjsjAB.getSj(mMis)).multiply(I).
225                      add(cjsjKH.getDsjDh(sMj).multiply(cjsjAB.getCj(mMis)).multiply(sgn(s - j)));
226         } else {
227             final int sMim = FastMath.abs(s - I * m);
228             dHmsdh = cjsjKH.getDcjDh(sMj).multiply(cjsjAB.getSj(sMim)).multiply(-sgn(s - m)).
229                      add(cjsjKH.getDsjDh(sMj).multiply(cjsjAB.getCj(sMim)).multiply(sgn(s - j)));
230         }
231         return dHmsdh;
232     }
233 
234     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>α</sub> coefficient.
235      * @param m m subscript
236      * @param s s subscript
237      * @param j order
238      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>α</sub>
239      */
240     public T getdHmsdAlpha(final int m, final int s, final int j) {
241         final int sMj  = FastMath.abs(s - j);
242         final T dHmsdAl;
243         if (FastMath.abs(s) <= m) {
244             final int mMis = m - I * s;
245             dHmsdAl = cjsjKH.getCj(sMj).multiply(cjsjAB.getDsjDk(mMis)).multiply(I).
246                       add(cjsjKH.getSj(sMj).multiply(cjsjAB.getDcjDk(mMis)).multiply(sgn(s - j)));
247         } else {
248             final int sMim = FastMath.abs(s - I * m);
249             dHmsdAl = cjsjKH.getCj(sMj).multiply(cjsjAB.getDsjDk(sMim)).multiply(-sgn(s - m)).
250                       add(cjsjKH.getSj(sMj).multiply(cjsjAB.getDcjDk(sMim)).multiply(sgn(s - j)));
251         }
252         return dHmsdAl;
253     }
254 
255     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>β</sub> coefficient.
256      * @param m m subscript
257      * @param s s subscript
258      * @param j order
259      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>β</sub>
260      */
261     public T getdHmsdBeta(final int m, final int s, final int j) {
262         final int sMj = FastMath.abs(s - j);
263         final T dHmsdBe;
264         if (FastMath.abs(s) <= m) {
265             final int mMis = m - I * s;
266             dHmsdBe = cjsjKH.getCj(sMj).multiply(cjsjAB.getDsjDh(mMis)).multiply(I).
267                       add(cjsjKH.getSj(sMj).multiply(cjsjAB.getDcjDh(mMis)).multiply(sgn(s - j)));
268         } else {
269             final int sMim = FastMath.abs(s - I * m);
270             dHmsdBe = cjsjKH.getCj(sMj).multiply(cjsjAB.getDsjDh(sMim)).multiply(-sgn(s - m)).
271                       add(cjsjKH.getSj(sMj).multiply(cjsjAB.getDcjDh(sMim)).multiply(sgn(s - j)));
272         }
273         return dHmsdBe;
274     }
275 
276     /** Get the sign of an integer.
277      *  @param i number on which evaluation is done
278      *  @return -1 or +1 depending on sign of i
279      */
280     private int sgn(final int i) {
281         return (i < 0) ? -1 : 1;
282     }
283 }