1   /* Copyright 2002-2013 CS Systèmes d'Information
2    * Licensed to CS Systèmes d'Information (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 java.util.ArrayList;
20  import java.util.List;
21  
22  import org.apache.commons.math3.complex.Complex;
23  import org.apache.commons.math3.util.FastMath;
24  
25  /** Compute the G<sub>ms</sub><sup>j</sup> and the H<sub>ms</sub><sup>j</sup>
26   *  polynomials in the equinoctial elements h, k and the direction cosines &alpha; and &beta;
27   *  and their partial derivatives with respect to k, h, &alpha; and &beta;.
28   *  <p>
29   *  The expressions used are equations 2.7.5-(1)(2) from the Danielson paper.
30   *  </p>
31   *  @author Romain Di Costanzo
32   */
33  public class GHmsjPolynomials {
34  
35      /** C<sub>j</sub>(k, h), S<sub>j</sub>(k, h) coefficient.
36       * (k, h) are the (x, y) component of the eccentricity vector in equinoctial elements
37       */
38      private final CjSjCoefficient cjsjKH;
39  
40      /** C<sub>j</sub>(&alpha;, &beta;), S<sub>j</sub>(&alpha;, &beta;) coefficient.
41       * (&alpha;, &beta;) are the direction cosines
42       */
43      private final CjSjCoefficient cjsjAB;
44  
45      /** Is the orbit represented as a retrograde orbit.
46       *  I = -1 if yes, +1 otherwise.
47       */
48      private int                   I;
49  
50      /** Create a set of G<sub>ms</sub><sup>j</sup> and H<sub>ms</sub><sup>j</sup> polynomials.
51       *  @param k X component of the eccentricity vector
52       *  @param h Y component of the eccentricity vector
53       *  @param alpha direction cosine &alpha;
54       *  @param beta direction cosine &beta;
55       *  @param retroFactor -1 if the orbit is represented as retrograde, +1 otherwise
56       **/
57      public GHmsjPolynomials(final double k, final double h,
58                              final double alpha, final double beta,
59                              final int retroFactor) {
60          this.cjsjKH = new CjSjCoefficient(k, h);
61          this.cjsjAB = new CjSjCoefficient(alpha, beta);
62          this.I = retroFactor;
63      }
64  
65      /** Get the G<sub>ms</sub><sup>j</sup> coefficient.
66       * @param m m subscript
67       * @param s s subscript
68       * @param j order
69       * @return the G<sub>ms</sub><sup>j</sup>
70       */
71      public double getGmsj(final int m, final int s, final int j) {
72          final int sMj = FastMath.abs(s - j);
73          double gms = 0d;
74          if (FastMath.abs(s) <= m) {
75              final int mMis = m - I * s;
76              gms = cjsjKH.getCj(sMj) * cjsjAB.getCj(mMis) -
77                    I * sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getSj(mMis);
78          } else {
79              final int sMim = FastMath.abs(s - I * m);
80              gms = cjsjKH.getCj(sMj) * cjsjAB.getCj(sMim) +
81                    sgn(s - j) * sgn(s - m) * cjsjKH.getSj(sMj) * cjsjAB.getSj(sMim);
82          }
83          return gms;
84      }
85  
86      /** Get the H<sub>ms</sub><sup>j</sup> coefficient.
87       * @param m m subscript
88       * @param s s subscript
89       * @param j order
90       * @return the H<sub>ms</sub><sup>j</sup>
91       */
92      public double getHmsj(final int m, final int s, final int j) {
93          final int sMj = FastMath.abs(s - j);
94          double hms = 0d;
95          if (FastMath.abs(s) <= m) {
96              final int mMis = m - I * s;
97              hms = I * cjsjKH.getCj(sMj) * cjsjAB.getSj(mMis) + sgn(s - j) *
98                    cjsjKH.getSj(sMj) * cjsjAB.getCj(mMis);
99          } else {
100             final int sMim = FastMath.abs(s - I * m);
101             hms = -sgn(s - m) * cjsjKH.getCj(sMj) * cjsjAB.getSj(sMim) +
102                    sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getCj(sMim);
103         }
104         return hms;
105     }
106 
107     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>k</sub> coefficient.
108      * @param m m subscript
109      * @param s s subscript
110      * @param j order
111      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>k</sub>
112      */
113     public double getdGmsdk(final int m, final int s, final int j) {
114         final int sMj = FastMath.abs(s - j);
115         double dGmsdk = 0d;
116         if (FastMath.abs(s) <= m) {
117             final int mMis = m - I * s;
118             dGmsdk = cjsjKH.getDcjDk(sMj) * cjsjAB.getCj(mMis) -
119                    I * sgn(s - j) * cjsjKH.getDsjDk(sMj) * cjsjAB.getSj(mMis);
120         } else {
121             final int sMim = FastMath.abs(s - I * m);
122             dGmsdk = cjsjKH.getDcjDk(sMj) * cjsjAB.getCj(sMim) +
123                     sgn(s - j) * sgn(s - m) * cjsjKH.getDsjDk(sMj) * cjsjAB.getSj(sMim);
124         }
125         return dGmsdk;
126     }
127 
128     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>h</sub> coefficient.
129      * @param m m subscript
130      * @param s s subscript
131      * @param j order
132      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>h</sub>
133      */
134     public double getdGmsdh(final int m, final int s, final int j) {
135         final int sMj = FastMath.abs(s - j);
136         double dGmsdh = 0d;
137         if (FastMath.abs(s) <= m) {
138             final int mMis = m - I * s;
139             dGmsdh = cjsjKH.getDcjDh(sMj) * cjsjAB.getCj(mMis) -
140                     I * sgn(s - j) * cjsjKH.getDsjDh(sMj) * cjsjAB.getSj(mMis);
141         } else {
142             final int sMim = FastMath.abs(s - I * m);
143             dGmsdh = cjsjKH.getDcjDh(sMj) * cjsjAB.getCj(sMim) +
144                     sgn(s - j) * sgn(s - m) * cjsjKH.getDsjDh(sMj) * cjsjAB.getSj(sMim);
145         }
146         return dGmsdh;
147     }
148 
149     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub> coefficient.
150      * @param m m subscript
151      * @param s s subscript
152      * @param j order
153      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub>
154      */
155     public double getdGmsdAlpha(final int m, final int s, final int j) {
156         final int sMj  = FastMath.abs(s - j);
157         double dGmsdAl = 0d;
158         if (FastMath.abs(s) <= m) {
159             final int mMis = m - I * s;
160             dGmsdAl = cjsjKH.getCj(sMj) * cjsjAB.getDcjDk(mMis) -
161                    I * sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDk(mMis);
162         } else {
163             final int sMim = FastMath.abs(s - I * m);
164             dGmsdAl = cjsjKH.getCj(sMj) * cjsjAB.getDcjDk(sMim) +
165                     sgn(s - j) * sgn(s - m) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDk(sMim);
166         }
167         return dGmsdAl;
168     }
169 
170     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub> coefficient.
171      * @param m m subscript
172      * @param s s subscript
173      * @param j order
174      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub>
175      */
176     public double getdGmsdBeta(final int m, final int s, final int j) {
177         final int sMj = FastMath.abs(s - j);
178         double dGmsdBe = 0d;
179         if (FastMath.abs(s) <= m) {
180             final int mMis = m - I * s;
181             dGmsdBe = cjsjKH.getCj(sMj) * cjsjAB.getDcjDh(mMis) -
182                     I * sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDh(mMis);
183         } else {
184             final int sMim = FastMath.abs(s - I * m);
185             dGmsdBe = cjsjKH.getCj(sMj) * cjsjAB.getDcjDh(sMim) +
186                     sgn(s - j) * sgn(s - m) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDh(sMim);
187         }
188         return dGmsdBe;
189     }
190 
191     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>k</sub> coefficient.
192      * @param m m subscript
193      * @param s s subscript
194      * @param j order
195      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>k</sub>
196      */
197     public double getdHmsdk(final int m, final int s, final int j) {
198         final int sMj = FastMath.abs(s - j);
199         double dHmsdk = 0d;
200         if (FastMath.abs(s) <= m) {
201             final int mMis = m - I * s;
202             dHmsdk = I * cjsjKH.getDcjDk(sMj) * cjsjAB.getSj(mMis) +
203                     sgn(s - j) * cjsjKH.getDsjDk(sMj) * cjsjAB.getCj(mMis);
204         } else {
205             final int sMim = FastMath.abs(s - I * m);
206             dHmsdk = -sgn(s - m) * cjsjKH.getDcjDk(sMj) * cjsjAB.getSj(sMim) +
207                     sgn(s - j) * cjsjKH.getDsjDk(sMj) * cjsjAB.getCj(sMim);
208         }
209         return dHmsdk;
210     }
211 
212     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>h</sub> coefficient.
213      * @param m m subscript
214      * @param s s subscript
215      * @param j order
216      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>h</sub>
217      */
218     public double getdHmsdh(final int m,  final int s, final int j) {
219         final int sMj = FastMath.abs(s - j);
220         double dHmsdh = 0d;
221         if (FastMath.abs(s) <= m) {
222             final int mMis = m - I * s;
223             dHmsdh = I * cjsjKH.getDcjDh(sMj) * cjsjAB.getSj(mMis) +
224                     sgn(s - j) * cjsjKH.getDsjDh(sMj) * cjsjAB.getCj(mMis);
225         } else {
226             final int sMim = FastMath.abs(s - I * m);
227             dHmsdh = -sgn(s - m) * cjsjKH.getDcjDh(sMj) * cjsjAB.getSj(sMim) +
228                     sgn(s - j) * cjsjKH.getDsjDh(sMj) * cjsjAB.getCj(sMim);
229         }
230         return dHmsdh;
231     }
232 
233     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub> coefficient.
234      * @param m m subscript
235      * @param s s subscript
236      * @param j order
237      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub>
238      */
239     public double getdHmsdAlpha(final int m, final int s, final int j) {
240         final int sMj  = FastMath.abs(s - j);
241         double dHmsdAl = 0d;
242         if (FastMath.abs(s) <= m) {
243             final int mMis = m - I * s;
244             dHmsdAl = I * cjsjKH.getCj(sMj) * cjsjAB.getDsjDk(mMis) +
245                     sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDk(mMis);
246         } else {
247             final int sMim = FastMath.abs(s - I * m);
248             dHmsdAl = -sgn(s - m) * cjsjKH.getCj(sMj) * cjsjAB.getDsjDk(sMim) +
249                     sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDk(sMim);
250         }
251         return dHmsdAl;
252     }
253 
254     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub> coefficient.
255      * @param m m subscript
256      * @param s s subscript
257      * @param j order
258      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub>
259      */
260     public double getdHmsdBeta(final int m, final int s, final int j) {
261         final int sMj = FastMath.abs(s - j);
262         double dHmsdBe = 0d;
263         if (FastMath.abs(s) <= m) {
264             final int mMis = m - I * s;
265             dHmsdBe = I * cjsjKH.getCj(sMj) * cjsjAB.getDsjDh(mMis) +
266                    sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDh(mMis);
267         } else {
268             final int sMim = FastMath.abs(s - I * m);
269             dHmsdBe = -sgn(s - m) * cjsjKH.getCj(sMj) * cjsjAB.getDsjDh(sMim) +
270                     sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDh(sMim);
271         }
272         return dHmsdBe;
273     }
274 
275     /** Get the sign of an integer.
276      *  @param i number on which evaluation is done
277      *  @return -1 or +1 depending on sign of i
278      */
279     private int sgn(final int i) {
280         return (i < 0) ? -1 : 1;
281     }
282 
283     /** Compute the S<sub>j</sub>(k, h) and the C<sub>j</sub>(k, h) series
284      *  and their partial derivatives with respect to k and h.
285      *  <p>
286      *  Those series are given in Danielson paper by expression 2.5.3-(5):
287      *  <pre>C<sub>j</sub>(k, h) + i S<sub>j</sub>(k, h) = (k+ih)<sup>j</sup> </pre>
288      *  </p>
289      *  The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) elements are store as an
290      *  {@link ArrayList} of {@link Complex} number, the C<sub>j</sub>(k, h) being
291      *  represented by the real and the S<sub>j</sub>(k, h) by the imaginary part.
292      */
293     private static class CjSjCoefficient {
294 
295         /** Last computed order j. */
296         private int jLast;
297 
298         /** Complex base (k + ih) of the C<sub>j</sub>, S<sub>j</sub> series. */
299         private final Complex kih;
300 
301         /** List of computed elements. */
302         private final List<Complex> cjsj;
303 
304         /** C<sub>j</sub>(k, h) and S<sub>j</sub>(k, h) constructor.
305          * @param k k value
306          * @param h h value
307          */
308         public CjSjCoefficient(final double k, final double h) {
309             kih  = new Complex(k, h);
310             cjsj = new ArrayList<Complex>();
311             cjsj.add(new Complex(1, 0));
312             cjsj.add(kih);
313             jLast = 1;
314         }
315 
316         /** Get the C<sub>j</sub> coefficient.
317          * @param j order
318          * @return C<sub>j</sub>
319          */
320         public double getCj(final int j) {
321             if (j > jLast) {
322                 // Update to order j
323                 updateCjSj(j);
324             }
325             return cjsj.get(j).getReal();
326         }
327 
328         /** Get the S<sub>j</sub> coefficient.
329          * @param j order
330          * @return S<sub>j</sub>
331          */
332         public double getSj(final int j) {
333             if (j > jLast) {
334                 // Update to order j
335                 updateCjSj(j);
336             }
337             return cjsj.get(j).getImaginary();
338         }
339 
340         /** Get the dC<sub>j</sub> / dk coefficient.
341          * @param j order
342          * @return dC<sub>j</sub> / d<sub>k</sub>
343          */
344         public double getDcjDk(final int j) {
345             return j == 0 ? 0 : j * getCj(j - 1);
346         }
347 
348         /** Get the dS<sub>j</sub> / dk coefficient.
349          * @param j order
350          * @return dS<sub>j</sub> / d<sub>k</sub>
351          */
352         public double getDsjDk(final int j) {
353             return j == 0 ? 0 : j * getSj(j - 1);
354         }
355 
356         /** Get the dC<sub>j</sub> / dh coefficient.
357          * @param j order
358          * @return dC<sub>i</sub> / d<sub>k</sub>
359          */
360         public double getDcjDh(final int j) {
361             return j == 0 ? 0 : -j * getSj(j - 1);
362         }
363 
364         /** Get the dS<sub>j</sub> / dh coefficient.
365          * @param j order
366          * @return dS<sub>j</sub> / d<sub>h</sub>
367          */
368         public double getDsjDh(final int j) {
369             return j == 0 ? 0 : j * getCj(j - 1);
370         }
371 
372         /** Update the cjsj up to order j.
373          * @param j order
374          */
375         private void updateCjSj(final int j) {
376             Complex last = cjsj.get(cjsj.size() - 1);
377             for (int i = jLast; i < j; i++) {
378                 final Complex next = last.multiply(kih);
379                 cjsj.add(next);
380                 last = next;
381             }
382             jLast = j;
383         }
384 
385     }
386 
387 }