1   /* Copyright 2002-2021 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 java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.Field;
23  import org.hipparchus.CalculusFieldElement;
24  import org.hipparchus.complex.Complex;
25  import org.hipparchus.exception.NullArgumentException;
26  
27  /** Compute the S<sub>j</sub>(k, h) and the C<sub>j</sub>(k, h) series
28   *  and their partial derivatives with respect to k and h.
29   *  <p>
30   *  Those series are given in Danielson paper by expression 2.5.3-(5):
31   *
32   *  <p> C<sub>j</sub>(k, h) + i S<sub>j</sub>(k, h) = (k+ih)<sup>j</sup>
33   *
34   *  <p>
35   *  The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) elements are store as an
36   *  {@link ArrayList} of {@link Complex} number, the C<sub>j</sub>(k, h) being
37   *  represented by the real and the S<sub>j</sub>(k, h) by the imaginary part.
38   */
39  public class FieldCjSjCoefficient <T extends CalculusFieldElement<T>> {
40  
41      /** Zero for initialization. /*/
42      private final T zero;
43  
44      /** Last computed order j. */
45      private int jLast;
46  
47      /** Complex base (k + ih) of the C<sub>j</sub>, S<sub>j</sub> series. */
48      private final FieldComplex<T> kih;
49  
50      /** List of computed elements. */
51      private final List<FieldComplex<T>> cjsj;
52  
53      /** C<sub>j</sub>(k, h) and S<sub>j</sub>(k, h) constructor.
54       * @param k k value
55       * @param h h value
56       * @param field field for fieldElements
57       */
58      public FieldCjSjCoefficient(final T k, final T h, final Field<T> field) {
59          zero = field.getZero();
60          kih  = new FieldComplex<>(k, h);
61          cjsj = new ArrayList<FieldComplex<T>>();
62          cjsj.add(new FieldComplex<>(zero.add(1.), zero));
63          cjsj.add(kih);
64          jLast = 1;
65      }
66  
67      /** Get the C<sub>j</sub> coefficient.
68       * @param j order
69       * @return C<sub>j</sub>
70       */
71      public T getCj(final int j) {
72          if (j > jLast) {
73              // Update to order j
74              updateCjSj(j);
75          }
76          return cjsj.get(j).getReal();
77      }
78  
79      /** Get the S<sub>j</sub> coefficient.
80       * @param j order
81       * @return S<sub>j</sub>
82       */
83      public T getSj(final int j) {
84          if (j > jLast) {
85              // Update to order j
86              updateCjSj(j);
87          }
88          return cjsj.get(j).getImaginary();
89      }
90  
91      /** Get the dC<sub>j</sub> / dk coefficient.
92       * @param j order
93       * @return dC<sub>j</sub> / d<sub>k</sub>
94       */
95      public T getDcjDk(final int j) {
96          return j == 0 ? zero : getCj(j - 1).multiply(j);
97      }
98  
99      /** Get the dS<sub>j</sub> / dk coefficient.
100      * @param j order
101      * @return dS<sub>j</sub> / d<sub>k</sub>
102      */
103     public T getDsjDk(final int j) {
104         return j == 0 ? zero : getSj(j - 1).multiply(j);
105     }
106 
107     /** Get the dC<sub>j</sub> / dh coefficient.
108      * @param j order
109      * @return dC<sub>i</sub> / d<sub>k</sub>
110      */
111     public T getDcjDh(final int j) {
112         return j == 0 ? zero : getSj(j - 1).multiply(-j);
113     }
114 
115     /** Get the dS<sub>j</sub> / dh coefficient.
116      * @param j order
117      * @return dS<sub>j</sub> / d<sub>h</sub>
118      */
119     public T getDsjDh(final int j) {
120         return j == 0 ? zero : getCj(j - 1).multiply(j);
121     }
122 
123     /** Update the cjsj up to order j.
124      * @param j order
125      */
126     private void updateCjSj(final int j) {
127         FieldComplex<T> last = cjsj.get(cjsj.size() - 1);
128         for (int i = jLast; i < j; i++) {
129             final FieldComplex<T> next = last.multiply(kih);
130             cjsj.add(next);
131             last = next;
132         }
133         jLast = j;
134     }
135 
136     private static class FieldComplex <T extends CalculusFieldElement<T>> {
137 
138         /** The imaginary part. */
139         private final T imaginary;
140 
141         /** The real part. */
142         private final T real;
143 
144         /**
145          * Create a complex number given the real and imaginary parts.
146          *
147          * @param real Real part.
148          * @param imaginary Imaginary part.
149          */
150         FieldComplex(final T real, final T imaginary) {
151             this.real = real;
152             this.imaginary = imaginary;
153         }
154 
155         /**
156          * Access the real part.
157          *
158          * @return the real part.
159          */
160         public T getReal() {
161             return real;
162         }
163 
164         /**
165          * Access the imaginary part.
166          *
167          * @return the imaginary part.
168          */
169         public T getImaginary() {
170             return imaginary;
171         }
172 
173         /**
174          * Create a complex number given the real and imaginary parts.
175          *
176          * @param realPart Real part.
177          * @param imaginaryPart Imaginary part.
178          * @return a new complex number instance.
179          *
180          * @see #valueOf(double, double)
181          */
182         protected FieldComplex<T> createComplex(final T realPart, final T imaginaryPart) {
183             return new FieldComplex<>(realPart, imaginaryPart);
184         }
185 
186         /**
187          * Returns a {@code Complex} whose value is {@code this * factor}.
188          * Implements preliminary checks for {@code NaN} and infinity followed by
189          * the definitional formula:
190          * <p>
191          *   {@code (a + bi)(c + di) = (ac - bd) + (ad + bc)i}
192          * </p>
193          * <p>
194          * Returns finite values in components of the result per the definitional
195          * formula in all remaining cases.</p>
196          *
197          * @param  factor value to be multiplied by this {@code Complex}.
198          * @return {@code this * factor}.
199          * @throws NullArgumentException if {@code factor} is {@code null}.
200          */
201         public FieldComplex<T> multiply(final FieldComplex<T> factor) throws NullArgumentException {
202             return createComplex(real.multiply(factor.real).subtract(imaginary.multiply(factor.imaginary)),
203                                  real.multiply(factor.imaginary).add(imaginary.multiply(factor.real)));
204         }
205     }
206 }