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 }