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 java.util.TreeMap;
20  
21  import org.apache.commons.math3.util.FastMath;
22  import org.orekit.errors.OrekitException;
23  import org.orekit.errors.OrekitMessages;
24  
25  /**
26   * This class is designed to provide coefficient from the DSST theory.
27   *
28   * @author Romain Di Costanzo
29   */
30  public class CoefficientsFactory {
31  
32      /** Internal storage of the polynomial values. Reused for further computation. */
33      private static TreeMap<NSKey, Double> VNS = new TreeMap<NSKey, Double>();
34  
35      /** Last computed order for V<sub>ns</sub> coefficients. */
36      private static int         LAST_VNS_ORDER = 2;
37  
38      /** Static initialization for the V<sub>ns</sub> coefficient. */
39      static {
40          // Initialization
41          VNS.put(new NSKey(0, 0), 1.);
42          VNS.put(new NSKey(1, 0), 0.);
43          VNS.put(new NSKey(1, 1), 0.5);
44      }
45  
46      /** Private constructor as the class is a utility class.
47       */
48      private CoefficientsFactory() {
49      }
50  
51      /** Compute the Q<sub>n,s</sub> coefficients evaluated at γ from the recurrence formula 2.8.3-(2).
52       *  <p>
53       *  Q<sub>n,s</sub> coefficients are computed for n = 0 to nMax
54       *  and s = 0 to sMax + 1 in order to also get the derivative dQ<sub>n,s</sub>/dγ = Q(n, s + 1)
55       *  </p>
56       *  @param gamma γ angle
57       *  @param nMax n max value
58       *  @param sMax s max value
59       *  @return Q<sub>n,s</sub> coefficients array
60       */
61      public static double[][] computeQns(final double gamma, final int nMax, final int sMax) {
62  
63          // Initialization
64          final int sDim = FastMath.min(sMax + 1, nMax) + 1;
65          final double[][] Qns = new double[nMax + 1][];
66          for (int i = 0; i <= nMax; i++) {
67              final int snDim = FastMath.min(i + 1, sDim);
68              Qns[i] = new double[snDim];
69          }
70  
71          // first element
72          Qns[0][0] = 1;
73  
74          for (int n = 1; n <= nMax; n++) {
75              final int snDim = FastMath.min(n + 1, sDim);
76              for (int s = 0; s < snDim; s++) {
77                  if (n == s) {
78                      Qns[n][s] = (2. * s - 1.) * Qns[s - 1][s - 1];
79                  } else if (n == (s + 1)) {
80                      Qns[n][s] = (2. * s + 1.) * gamma * Qns[s][s];
81                  } else {
82                      Qns[n][s] = (2. * n - 1.) * gamma * Qns[n - 1][s] - (n + s - 1.) * Qns[n - 2][s];
83                      Qns[n][s] /= n - s;
84                  }
85              }
86          }
87  
88          return Qns;
89      }
90  
91      /** Compute recursively G<sub>s</sub> and H<sub>s</sub> polynomials from equation 3.1-(5).
92       *  @param k x-component of the eccentricity vector
93       *  @param h y-component of the eccentricity vector
94       *  @param alpha 1st direction cosine
95       *  @param beta 2nd direction cosine
96       *  @param order development order
97       *  @return Array of G<sub>s</sub> and H<sub>s</sub> polynomials for s from 0 to order.<br>
98       *          The 1st column contains the G<sub>s</sub> values.
99       *          The 2nd column contains the H<sub>s</sub> values.
100      */
101     public static double[][] computeGsHs(final double k, final double h,
102                                          final double alpha, final double beta,
103                                          final int order) {
104         // Constant terms
105         final double hamkb = h * alpha - k * beta;
106         final double kaphb = k * alpha + h * beta;
107         // Initialization
108         final double[][] GsHs = new double[2][order + 1];
109         GsHs[0][0] = 1.;
110         GsHs[1][0] = 0.;
111 
112         for (int s = 1; s <= order; s++) {
113             // Gs coefficient
114             GsHs[0][s] = kaphb * GsHs[0][s - 1] - hamkb * GsHs[1][s - 1];
115             // Hs coefficient
116             GsHs[1][s] = hamkb * GsHs[0][s - 1] + kaphb * GsHs[1][s - 1];
117         }
118 
119         return GsHs;
120     }
121 
122     /** Compute the V<sub>n,s</sub> coefficients from 2.8.2-(1)(2).
123      * @param order Order of the computation. Computation will be done from 0 to order -1
124      * @return Map of the V<sub>n, s</sub> coefficients
125      */
126     public static TreeMap<NSKey, Double> computeVns(final int order) {
127 
128         if (order > LAST_VNS_ORDER) {
129             // Compute coefficient
130             // Need previous computation as recurrence relation is done at s + 1 and n + 2
131             final int min = (LAST_VNS_ORDER - 2 < 0) ? 0 : (LAST_VNS_ORDER - 2);
132             for (int n = min; n < order; n++) {
133                 for (int s = 0; s < n + 1; s++) {
134                     if ((n - s) % 2 != 0) {
135                         VNS.put(new NSKey(n, s), 0.);
136                     } else {
137                         // s = n
138                         if (n == s && (s + 1) < order) {
139                             VNS.put(new NSKey(s + 1, s + 1), VNS.get(new NSKey(s, s)) / (2 * s + 2.));
140                         }
141                         // otherwise
142                         if ((n + 2) < order) {
143                             VNS.put(new NSKey(n + 2, s), VNS.get(new NSKey(n, s)) * (-n + s - 1.) / (n + s + 2.));
144                         }
145                     }
146                 }
147             }
148             LAST_VNS_ORDER = order;
149         }
150         return VNS;
151     }
152 
153     /** Get the V<sub>n,s</sub><sup>m</sup> coefficient from V<sub>n,s</sub>.
154      *  <br>See § 2.8.2 in Danielson paper.
155      * @param m m
156      * @param n n
157      * @param s s
158      * @param fns (n + |s|)!
159      * @param fnm (n - m)!
160      * @return The V<sub>n, s</sub> <sup>m</sup> coefficient
161      * @throws OrekitException if m > n
162      */
163     public static double getVmns(final int m, final int n, final int s,
164                                  final double fns, final double fnm)
165         throws OrekitException {
166         if (m > n) {
167             throw new OrekitException(OrekitMessages.DSST_VMNS_COEFFICIENT_ERROR_MS, m, n);
168         }
169 
170         double result = 0;
171         // If (n - s) is odd, the Vmsn coefficient is null
172         if ((n - s) % 2 == 0) {
173             // Update the Vns coefficient
174             if ((n + 1) > LAST_VNS_ORDER) {
175                 computeVns(n + 1);
176             }
177             if (s >= 0) {
178                 result = fns  * VNS.get(new NSKey(n, s)) / fnm;
179             } else {
180                 // If s < 0 : Vmn-s = (-1)^(-s) Vmns
181                 final int mops = (s % 2 == 0) ? 1 : -1;
182                 result = mops * fns * VNS.get(new NSKey(n, -s)) / fnm;
183             }
184         }
185         return result;
186     }
187 
188     /** Key formed by two integer values. */
189     public static class NSKey implements Comparable<NSKey> {
190 
191         /** n value. */
192         private final int n;
193 
194         /** s value. */
195         private final int s;
196 
197         /** Simple constructor.
198          * @param n n
199          * @param s s
200          */
201         public NSKey(final int n, final int s) {
202             this.n = n;
203             this.s = s;
204         }
205 
206         /** Get n.
207          * @return n
208          */
209         public int getN() {
210             return n;
211         }
212 
213         /** Get s.
214          * @return s
215          */
216         public int getS() {
217             return s;
218         }
219 
220         /** {@inheritDoc} */
221         public int compareTo(final NSKey key) {
222             int result = 1;
223             if (n == key.n) {
224                 if (s < key.s) {
225                     result = -1;
226                 } else if (s == key.s) {
227                     result = 0;
228                 }
229             } else if (n < key.n) {
230                 result = -1;
231             }
232             return result;
233         }
234 
235         /** {@inheritDoc} */
236         public boolean equals(final Object key) {
237 
238             if (key == this) {
239                 // first fast check
240                 return true;
241             }
242 
243             if ((key != null) && (key instanceof NSKey)) {
244                 return (n == ((NSKey) key).n) && (s == ((NSKey) key).s);
245             }
246 
247             return false;
248 
249         }
250 
251         /** {@inheritDoc} */
252         public int hashCode() {
253             return 0x998493a6 ^ (n << 8) ^ s;
254         }
255 
256     }
257 
258     /** MNS couple's key. */
259     public static class MNSKey implements Comparable<MNSKey> {
260 
261         /** m value. */
262         private final int m;
263 
264         /** n value. */
265         private final int n;
266 
267         /** s value. */
268         private final int s;
269 
270         /** Simpleconstructor.
271          * @param m m
272          * @param n n
273          * @param s s
274          */
275         public MNSKey(final int m, final int n, final int s) {
276             this.m = m;
277             this.n = n;
278             this.s = s;
279         }
280 
281         /** {@inheritDoc} */
282         public int compareTo(final MNSKey key) {
283             int result = 1;
284             if (m == key.m) {
285                 if (n == key.n) {
286                     if (s < key.s) {
287                         result = -1;
288                     } else if (s == key.s) {
289                         result = 0;
290                     } else {
291                         result = 1;
292                     }
293                 } else if (n < key.n) {
294                     result = -1;
295                 } else {
296                     result = 1;
297                 }
298             } else if (m < key.m) {
299                 result = -1;
300             }
301             return result;
302         }
303 
304         /** {@inheritDoc} */
305         public boolean equals(final Object key) {
306 
307             if (key == this) {
308                 // first fast check
309                 return true;
310             }
311 
312             if ((key != null) && (key instanceof MNSKey)) {
313                 return (m == ((MNSKey) key).m) &&
314                        (n == ((MNSKey) key).n) &&
315                        (s == ((MNSKey) key).s);
316             }
317 
318             return false;
319 
320         }
321 
322         /** {@inheritDoc} */
323         public int hashCode() {
324             return 0x25baa451 ^ (m << 16) ^ (n << 8) ^ s;
325         }
326 
327     }
328 
329 }