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