1 /* Copyright 2002-2025 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.forces.gravity.potential;
18
19 import org.hipparchus.util.SinCos;
20
21 /** Time-dependent part of a single component of spherical harmonics.
22 * @author Luc Maisonobe
23 * @since 11.1
24 */
25 class TimeDependentHarmonic {
26
27 /** Index of the trend reference in the gravity field. */
28 private final int trendReferenceIndex;
29
30 /** Base of the cosine coefficient. */
31 private final double cBase;
32
33 /** Base of the sine coefficient. */
34 private final double sBase;
35
36 /** Secular trend of the cosine coefficient. */
37 private double cTrend;
38
39 /** Secular trend of the sine coefficient. */
40 private double sTrend;
41
42 /** Indices of the reference dates in the gravity field. */
43 private int[] cosReferenceIndices;
44
45 /** Indices of the harmonic pulsations in the gravity field. */
46 private int[] cosPulsationIndices;
47
48 /** Cosine component of the cosine coefficient. */
49 private double[] cosC;
50
51 /** Cosine component of the sine coefficient. */
52 private double[] cosS;
53
54 /** Indices of the reference dates in the gravity field. */
55 private int[] sinReferenceIndices;
56
57 /** Indices of the harmonic pulsations in the gravity field. */
58 private int[] sinPulsationIndices;
59
60 /** Sine component of the cosine coefficient. */
61 private double[] sinC;
62
63 /** Sine component of the sine coefficient. */
64 private double[] sinS;
65
66 /** Build a part with only base.
67 * @param trendReferenceIndex index of the trend reference in the gravity field
68 * @param cBase base of the cosine coefficient
69 * @param sBase base of the sine coefficient
70 */
71 TimeDependentHarmonic(final int trendReferenceIndex, final double cBase, final double sBase) {
72 this(trendReferenceIndex, cBase, sBase, 0, 0);
73 }
74
75 /** Build a rescaled component.
76 * @param scale scaling factor to apply to all coefficients elements
77 * @param original original component
78 */
79 TimeDependentHarmonic(final double scale, final TimeDependentHarmonic original) {
80
81 // rescale base
82 this(original.trendReferenceIndex, scale * original.cBase, scale * original.sBase,
83 original.cosReferenceIndices.length, original.sinReferenceIndices.length);
84
85 // rescale trend
86 cTrend = scale * original.cTrend;
87 sTrend = scale * original.sTrend;
88
89 // rescale cosine
90 for (int i = 0; i < cosReferenceIndices.length; ++i) {
91 cosReferenceIndices[i] = original.cosReferenceIndices[i];
92 cosPulsationIndices[i] = original.cosPulsationIndices[i];
93 cosC[i] = scale * original.cosC[i];
94 cosS[i] = scale * original.cosS[i];
95 }
96
97 // rescale sine
98 for (int i = 0; i < sinReferenceIndices.length; ++i) {
99 sinReferenceIndices[i] = original.sinReferenceIndices[i];
100 sinPulsationIndices[i] = original.sinPulsationIndices[i];
101 sinC[i] = scale * original.sinC[i];
102 sinS[i] = scale * original.sinS[i];
103 }
104
105 }
106
107 /** Build a part with only base.
108 * @param trendReferenceIndex index of the trend reference in the gravity field
109 * @param cBase base of the cosine coefficient
110 * @param sBase base of the sine coefficient
111 * @param cSize initial size of the cosine arrays
112 * @param sSize initial size of the sine arrays
113 */
114 private TimeDependentHarmonic(final int trendReferenceIndex, final double cBase, final double sBase,
115 final int cSize, final int sSize) {
116
117 // linear part
118 this.trendReferenceIndex = trendReferenceIndex;
119 this.cBase = cBase;
120 this.sBase = sBase;
121 this.cTrend = 0.0;
122 this.sTrend = 0.0;
123
124 // cosine component
125 this.cosReferenceIndices = new int[cSize];
126 this.cosPulsationIndices = new int[cSize];
127 this.cosC = new double[cSize];
128 this.cosS = new double[cSize];
129
130 // sine component
131 this.sinReferenceIndices = new int[sSize];
132 this.sinPulsationIndices = new int[sSize];
133 this.sinC = new double[sSize];
134 this.sinS = new double[sSize];
135
136 }
137
138 /** Set the trend part.
139 * @param cDot secular trend of the cosine coefficient (s⁻¹)
140 * @param sDot secular trend of the sine coefficient (s⁻¹)
141 */
142 public void setTrend(final double cDot, final double sDot) {
143 this.cTrend = cDot;
144 this.sTrend = sDot;
145 }
146
147 /** Add a cosine component.
148 * @param cosReferenceIndex index of the reference date in the gravity field
149 * (if negative, use the trend reference index)
150 * @param cosPulsationIndex index of the harmonic pulsation in the gravity field
151 * @param cosineC cosine component of the cosine coefficient
152 * @param cosineS cosine component of the sine coefficient
153 */
154 public void addCosine(final int cosReferenceIndex, final int cosPulsationIndex,
155 final double cosineC, final double cosineS) {
156 final int refIndex = cosReferenceIndex < 0 ? trendReferenceIndex : cosReferenceIndex;
157 this.cosReferenceIndices = addInt(refIndex, this.cosReferenceIndices);
158 this.cosPulsationIndices = addInt(cosPulsationIndex, this.cosPulsationIndices);
159 this.cosC = addDouble(cosineC, this.cosC);
160 this.cosS = addDouble(cosineS, this.cosS);
161 }
162
163 /** Add a sine component.
164 * @param sinReferenceIndex index of the reference date in the gravity field
165 * (if negative, use the trend reference index)
166 * @param sinPulsationIndex index of the harmonic pulsation in the gravity field
167 * @param sineC sine component of the cosine coefficient
168 * @param sineS sine component of the sine coefficient
169 */
170 public void addSine(final int sinReferenceIndex, final int sinPulsationIndex,
171 final double sineC, final double sineS) {
172 final int refIndex = sinReferenceIndex < 0 ? trendReferenceIndex : sinReferenceIndex;
173 this.sinReferenceIndices = addInt(refIndex, this.sinReferenceIndices);
174 this.sinPulsationIndices = addInt(sinPulsationIndex, this.sinPulsationIndices);
175 this.sinC = addDouble(sineC, this.sinC);
176 this.sinS = addDouble(sineS, this.sinS);
177 }
178
179 /** Add an integer to an array.
180 * <p>
181 * Expanding the array one element at a time may seem a waste of time,
182 * but we expect the array to be 0, 1 or 2 elements long only, and this
183 * if performed only when reading gravity field, so its is worth doing
184 * it this way.
185 * </p>
186 * @param n integer to add
187 * @param array array where to add the integer
188 * @return new array
189 */
190 private static int[] addInt(final int n, final int[] array) {
191 final int[] newArray = new int[array.length + 1];
192 System.arraycopy(array, 0, newArray, 0, array.length);
193 newArray[array.length] = n;
194 return newArray;
195 }
196
197 /** Add a double number to an array.
198 * <p>
199 * Expanding the array one element at a time may seem a waste of time,
200 * but we expect the array to be 0, 1 or 2 elements long only, and this
201 * if performed only when reading gravity field, so its is worth doing
202 * it this way.
203 * </p>
204 * @param d double number to add
205 * @param array array where to add the double number
206 * @return new array
207 */
208 private static double[] addDouble(final double d, final double[] array) {
209 final double[] newArray = new double[array.length + 1];
210 System.arraycopy(array, 0, newArray, 0, array.length);
211 newArray[array.length] = d;
212 return newArray;
213 }
214
215 /** Compute the time-dependent part of a spherical harmonic cosine coefficient.
216 * @param offsets offsets to reference dates in the gravity field
217 * @param pulsations angular pulsations in the gravity field
218 * @return raw coefficient Cnm
219 */
220 public double computeCnm(final double[] offsets, final SinCos[][] pulsations) {
221
222 // trend effect
223 double cnm = cBase + offsets[trendReferenceIndex] * cTrend;
224
225 for (int i = 0; i < cosPulsationIndices.length; ++i) {
226 // cosine effect
227 cnm += cosC[i] * pulsations[cosReferenceIndices[i]][cosPulsationIndices[i]].cos();
228 }
229
230 for (int i = 0; i < sinPulsationIndices.length; ++i) {
231 // sine effect
232 cnm += sinC[i] * pulsations[sinReferenceIndices[i]][sinPulsationIndices[i]].sin();
233 }
234
235 return cnm;
236
237 }
238
239 /** Compute the time-dependent part of a spherical harmonic sine coefficient.
240 * @param offsets offsets to reference dates in the gravity field
241 * @param pulsations angular pulsations in the gravity field
242 * @return raw coefficient Snm
243 */
244 public double computeSnm(final double[] offsets, final SinCos[][] pulsations) {
245
246 // trend effect
247 double snm = sBase + offsets[trendReferenceIndex] * sTrend;
248
249 for (int i = 0; i < cosPulsationIndices.length; ++i) {
250 // cosine effect
251 snm += cosS[i] * pulsations[cosReferenceIndices[i]][cosPulsationIndices[i]].cos();
252 }
253
254 for (int i = 0; i < sinPulsationIndices.length; ++i) {
255 // sine effect
256 snm += sinS[i] * pulsations[sinReferenceIndices[i]][sinPulsationIndices[i]].sin();
257 }
258
259 return snm;
260
261 }
262
263 }