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.models.earth.atmosphere;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21 import org.hipparchus.geometry.euclidean.threed.Vector3D;
22 import org.hipparchus.util.FastMath;
23 import org.hipparchus.util.Precision;
24 import org.hipparchus.util.SinCos;
25 import org.orekit.bodies.OneAxisEllipsoid;
26 import org.orekit.errors.OrekitException;
27 import org.orekit.errors.OrekitMessages;
28 import org.orekit.frames.Frame;
29 import org.orekit.time.AbsoluteDate;
30 import org.orekit.time.FieldAbsoluteDate;
31 import org.orekit.utils.ExtendedPositionProvider;
32
33
34 /** This atmosphere model is the realization of the Modified Harris-Priester model.
35 * <p>
36 * This model is a static one that takes into account the diurnal density bulge.
37 * It doesn't need any space weather data but a density vs. altitude table, which
38 * depends on solar activity.
39 * </p>
40 * <p>
41 * The implementation relies on the book:<br>
42 * <b>Satellite Orbits</b><br>
43 * <i>Oliver Montenbruck, Eberhard Gill</i><br>
44 * Springer 2005
45 * </p>
46 * @author Pascal Parraud
47 */
48 public class HarrisPriester extends AbstractSunInfluencedAtmosphere {
49
50 // Constants :
51
52 /** Default cosine exponent value. */
53 private static final int N_DEFAULT = 4;
54
55 /** Minimal value for calculating poxer of cosine. */
56 private static final double MIN_COS = 1.e-12;
57
58 /** Lag angle for diurnal bulge. */
59 private static final double LAG = FastMath.toRadians(30.0);
60
61 /** Lag angle sine and cosine. */
62 private static final SinCos SCLAG = FastMath.sinCos(LAG);
63
64 // CHECKSTYLE: stop NoWhitespaceAfter check
65 /** Harris-Priester min-max density (kg/m3) vs. altitude (m) table.
66 * These data are valid for a mean solar activity. */
67 private static final double[][] ALT_RHO = {
68 { 100000.0, 4.974e-07, 4.974e-07 },
69 { 120000.0, 2.490e-08, 2.490e-08 },
70 { 130000.0, 8.377e-09, 8.710e-09 },
71 { 140000.0, 3.899e-09, 4.059e-09 },
72 { 150000.0, 2.122e-09, 2.215e-09 },
73 { 160000.0, 1.263e-09, 1.344e-09 },
74 { 170000.0, 8.008e-10, 8.758e-10 },
75 { 180000.0, 5.283e-10, 6.010e-10 },
76 { 190000.0, 3.617e-10, 4.297e-10 },
77 { 200000.0, 2.557e-10, 3.162e-10 },
78 { 210000.0, 1.839e-10, 2.396e-10 },
79 { 220000.0, 1.341e-10, 1.853e-10 },
80 { 230000.0, 9.949e-11, 1.455e-10 },
81 { 240000.0, 7.488e-11, 1.157e-10 },
82 { 250000.0, 5.709e-11, 9.308e-11 },
83 { 260000.0, 4.403e-11, 7.555e-11 },
84 { 270000.0, 3.430e-11, 6.182e-11 },
85 { 280000.0, 2.697e-11, 5.095e-11 },
86 { 290000.0, 2.139e-11, 4.226e-11 },
87 { 300000.0, 1.708e-11, 3.526e-11 },
88 { 320000.0, 1.099e-11, 2.511e-11 },
89 { 340000.0, 7.214e-12, 1.819e-11 },
90 { 360000.0, 4.824e-12, 1.337e-11 },
91 { 380000.0, 3.274e-12, 9.955e-12 },
92 { 400000.0, 2.249e-12, 7.492e-12 },
93 { 420000.0, 1.558e-12, 5.684e-12 },
94 { 440000.0, 1.091e-12, 4.355e-12 },
95 { 460000.0, 7.701e-13, 3.362e-12 },
96 { 480000.0, 5.474e-13, 2.612e-12 },
97 { 500000.0, 3.916e-13, 2.042e-12 },
98 { 520000.0, 2.819e-13, 1.605e-12 },
99 { 540000.0, 2.042e-13, 1.267e-12 },
100 { 560000.0, 1.488e-13, 1.005e-12 },
101 { 580000.0, 1.092e-13, 7.997e-13 },
102 { 600000.0, 8.070e-14, 6.390e-13 },
103 { 620000.0, 6.012e-14, 5.123e-13 },
104 { 640000.0, 4.519e-14, 4.121e-13 },
105 { 660000.0, 3.430e-14, 3.325e-13 },
106 { 680000.0, 2.632e-14, 2.691e-13 },
107 { 700000.0, 2.043e-14, 2.185e-13 },
108 { 720000.0, 1.607e-14, 1.779e-13 },
109 { 740000.0, 1.281e-14, 1.452e-13 },
110 { 760000.0, 1.036e-14, 1.190e-13 },
111 { 780000.0, 8.496e-15, 9.776e-14 },
112 { 800000.0, 7.069e-15, 8.059e-14 },
113 { 840000.0, 4.680e-15, 5.741e-14 },
114 { 880000.0, 3.200e-15, 4.210e-14 },
115 { 920000.0, 2.210e-15, 3.130e-14 },
116 { 960000.0, 1.560e-15, 2.360e-14 },
117 { 1000000.0, 1.150e-15, 1.810e-14 }
118 };
119 // CHECKSTYLE: resume NoWhitespaceAfter check
120
121 /** Cosine exponent from 2 to 6 according to inclination. */
122 private double n;
123
124 /** Earth body shape. */
125 private OneAxisEllipsoid earth;
126
127 /** Density table. */
128 private double[][] tabAltRho;
129
130 /** Simple constructor for Modified Harris-Priester atmosphere model.
131 * <p>The cosine exponent value is set to 4 by default.</p>
132 * <p>The default embedded density table is the one given in the referenced
133 * book from Montenbruck & Gill. It is given for mean solar activity and
134 * spreads over 100 to 1000 km.</p>
135 * @param sun the sun position
136 * @param earth the earth body shape
137 */
138 public HarrisPriester(final ExtendedPositionProvider sun,
139 final OneAxisEllipsoid earth) {
140 this(sun, earth, ALT_RHO, N_DEFAULT);
141 }
142
143 /** Constructor for Modified Harris-Priester atmosphere model.
144 * <p>Recommanded values for the cosine exponent spread over the range
145 * 2, for low inclination orbits, to 6, for polar orbits.</p>
146 * <p> The default embedded density table is the one given in the referenced
147 * book from Montenbruck & Gill. It is given for mean solar activity and
148 * spreads over 100 to 1000 km. </p>
149 * @param sun the sun position
150 * @param earth the earth body shape
151 * @param n the cosine exponent
152 */
153 public HarrisPriester(final ExtendedPositionProvider sun,
154 final OneAxisEllipsoid earth,
155 final double n) {
156 this(sun, earth, ALT_RHO, n);
157 }
158
159 /** Constructor for Modified Harris-Priester atmosphere model.
160 * <p>The provided density table must be an array such as:
161 * <ul>
162 * <li>tabAltRho[][0] = altitude (m)</li>
163 * <li>tabAltRho[][1] = min density (kg/m³)</li>
164 * <li>tabAltRho[][2] = max density (kg/m³)</li>
165 * </ul>
166 * <p> The altitude must be increasing without limitation in range. The
167 * internal density table is a copy of the provided one.
168 *
169 * <p>The cosine exponent value is set to 4 by default.</p>
170 * @param sun the sun position
171 * @param earth the earth body shape
172 * @param tabAltRho the density table
173 */
174 public HarrisPriester(final ExtendedPositionProvider sun,
175 final OneAxisEllipsoid earth,
176 final double[][] tabAltRho) {
177 this(sun, earth, tabAltRho, N_DEFAULT);
178 }
179
180 /** Constructor for Modified Harris-Priester atmosphere model.
181 * <p>Recommanded values for the cosine exponent spread over the range
182 * 2, for low inclination orbits, to 6, for polar orbits.</p>
183 * <p>The provided density table must be an array such as:
184 * <ul>
185 * <li>tabAltRho[][0] = altitude (m)</li>
186 * <li>tabAltRho[][1] = min density (kg/m³)</li>
187 * <li>tabAltRho[][2] = max density (kg/m³)</li>
188 * </ul>
189 * <p> The altitude must be increasing without limitation in range. The
190 * internal density table is a copy of the provided one.
191 *
192 * @param sun the sun position
193 * @param earth the earth body shape
194 * @param tabAltRho the density table
195 * @param n the cosine exponent
196 */
197 public HarrisPriester(final ExtendedPositionProvider sun,
198 final OneAxisEllipsoid earth,
199 final double[][] tabAltRho,
200 final double n) {
201 super(sun);
202 this.earth = earth;
203 setTabDensity(tabAltRho);
204 setN(n);
205 }
206
207 /** {@inheritDoc} */
208 public Frame getFrame() {
209 return earth.getBodyFrame();
210 }
211
212 /** Set parameter N, the cosine exponent.
213 * @param n the cosine exponent
214 */
215 private void setN(final double n) {
216 this.n = n;
217 }
218
219 /** Set a user define density table to deal with different solar activities.
220 * @param tab density vs. altitude table
221 */
222 private void setTabDensity(final double[][] tab) {
223 this.tabAltRho = new double[tab.length][];
224 for (int i = 0; i < tab.length; i++) {
225 this.tabAltRho[i] = tab[i].clone();
226 }
227 }
228
229 /** Get the current density table.
230 * <p>The density table is an array such as:
231 * <ul>
232 * <li>tabAltRho[][0] = altitude (m)</li>
233 * <li>tabAltRho[][1] = min density (kg/m³)</li>
234 * <li>tabAltRho[][2] = max density (kg/m³)</li>
235 * </ul>
236 * <p> The altitude must be increasing without limitation in range.
237 *
238 * <p>
239 * The returned density table is a copy of the current one.
240 * </p>
241 * @return density vs. altitude table
242 */
243 public double[][] getTabDensity() {
244 final double[][] copy = new double[tabAltRho.length][];
245 for (int i = 0; i < tabAltRho.length; i++) {
246 copy[i] = tabAltRho[i].clone();
247 }
248 return copy;
249 }
250
251 /** Get the minimal altitude for the model.
252 * <p>No computation is possible below this altitude.</p>
253 * @return the minimal altitude (m)
254 */
255 public double getMinAlt() {
256 return tabAltRho[0][0];
257 }
258
259 /** Get the maximal altitude for the model.
260 * <p>Above this altitude, density is assumed to be zero.</p>
261 * @return the maximal altitude (m)
262 */
263 public double getMaxAlt() {
264 return tabAltRho[tabAltRho.length - 1][0];
265 }
266
267 /** Get the local density.
268 * @param sunInEarth position of the Sun in Earth frame (m)
269 * @param posInEarth target position in Earth frame (m)
270 * @return the local density (kg/m³)
271 */
272 public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth) {
273
274 final double posAlt = getHeight(posInEarth);
275 // Check for height boundaries
276 if (posAlt < getMinAlt()) {
277 throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
278 }
279 if (posAlt > getMaxAlt()) {
280 return 0.;
281 }
282
283 // Diurnal bulge apex direction
284 final Vector3D sunDir = sunInEarth.normalize();
285 final Vector3D bulDir = new Vector3D(sunDir.getX() * SCLAG.cos() - sunDir.getY() * SCLAG.sin(),
286 sunDir.getX() * SCLAG.sin() + sunDir.getY() * SCLAG.cos(),
287 sunDir.getZ());
288
289 // Cosine of angle Psi between the diurnal bulge apex and the satellite
290 final double cosPsi = bulDir.normalize().dotProduct(posInEarth.normalize());
291 // (1 + cos(Psi))/2 = cos²(Psi/2)
292 final double c2Psi2 = (1. + cosPsi) / 2.;
293 final double cPsi2 = FastMath.sqrt(c2Psi2);
294 final double cosPow = (cPsi2 > MIN_COS) ? c2Psi2 * FastMath.pow(cPsi2, n - 2) : 0.;
295
296 // Search altitude index in density table
297 int ia = 0;
298 while (ia < tabAltRho.length - 2 && posAlt > tabAltRho[ia + 1][0]) {
299 ia++;
300 }
301
302 // Fractional satellite height
303 final double dH = (tabAltRho[ia][0] - posAlt) / (tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
304
305 // Min exponential density interpolation
306 final double rhoMin = tabAltRho[ia][1] * FastMath.pow(tabAltRho[ia + 1][1] / tabAltRho[ia][1], dH);
307
308 if (Precision.equals(cosPow, 0.)) {
309 return rhoMin;
310 } else {
311 // Max exponential density interpolation
312 final double rhoMax = tabAltRho[ia][2] * FastMath.pow(tabAltRho[ia + 1][2] / tabAltRho[ia][2], dH);
313 return rhoMin + (rhoMax - rhoMin) * cosPow;
314 }
315
316 }
317
318 /** Get the local density.
319 * @param sunInEarth position of the Sun in Earth frame (m)
320 * @param posInEarth target position in Earth frame (m)
321 * @return the local density (kg/m³)
322 * @param <T> instance of CalculusFieldElement<T>
323 */
324 public <T extends CalculusFieldElement<T>> T getDensity(final FieldVector3D<T> sunInEarth, final FieldVector3D<T> posInEarth) {
325 final T zero = posInEarth.getX().getField().getZero();
326 final T posAlt = getHeight(posInEarth);
327 // Check for height boundaries
328 if (posAlt.getReal() < getMinAlt()) {
329 throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
330 }
331 if (posAlt.getReal() > getMaxAlt()) {
332 return zero;
333 }
334
335 // Diurnal bulge apex direction
336 final FieldVector3D<T> sunDir = sunInEarth.normalize();
337 final FieldVector3D<T> bulDir = new FieldVector3D<>(sunDir.getX().multiply(SCLAG.cos()).subtract(sunDir.getY().multiply(SCLAG.sin())),
338 sunDir.getX().multiply(SCLAG.sin()).add(sunDir.getY().multiply(SCLAG.cos())),
339 sunDir.getZ());
340
341 // Cosine of angle Psi between the diurnal bulge apex and the satellite
342 final T cosPsi = posInEarth.normalize().dotProduct(bulDir.normalize());
343 // (1 + cos(Psi))/2 = cos²(Psi/2)
344 final T c2Psi2 = cosPsi.add(1.).divide(2);
345 final T cPsi2 = c2Psi2.sqrt();
346 final T cosPow = (cPsi2.getReal() > MIN_COS) ? c2Psi2.multiply(cPsi2.pow(n - 2)) : zero;
347
348 // Search altitude index in density table
349 int ia = 0;
350 while (ia < tabAltRho.length - 2 && posAlt.getReal() > tabAltRho[ia + 1][0]) {
351 ia++;
352 }
353
354 // Fractional satellite height
355 final T dH = posAlt.negate().add(tabAltRho[ia][0]).divide(tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
356
357 // Min exponential density interpolation
358 final T rhoMin = zero.newInstance(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);
359
360 if (Precision.equals(cosPow.getReal(), 0.)) {
361 return rhoMin;
362 } else {
363 // Max exponential density interpolation
364 final T rhoMax = zero.newInstance(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
365 return rhoMin.add(rhoMax.subtract(rhoMin).multiply(cosPow));
366 }
367
368 }
369
370 /** Get the local density at some position.
371 * @param date current date
372 * @param position current position
373 * @param frame the frame in which is defined the position
374 * @return local density (kg/m³)
375 * or if altitude is below the model minimal altitude
376 */
377 public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame) {
378
379 // Sun position in earth frame
380 final Vector3D sunInEarth = getSunPosition(date, earth.getBodyFrame());
381
382 // Target position in earth frame
383 final Vector3D posInEarth = frame
384 .getStaticTransformTo(earth.getBodyFrame(), date)
385 .transformPosition(position);
386
387 return getDensity(sunInEarth, posInEarth);
388 }
389
390 /** Get the local density at some position.
391 * @param date current date
392 * @param position current position
393 * @param <T> implements a CalculusFieldElement
394 * @param frame the frame in which is defined the position
395 * @return local density (kg/m³)
396 * or if altitude is below the model minimal altitude
397 */
398 public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
399 final FieldVector3D<T> position, final Frame frame) {
400 // Sun position in earth frame
401 final FieldVector3D<T> sunInEarth = getSunPosition(date, earth.getBodyFrame());
402
403 // Target position in earth frame
404 final FieldVector3D<T> posInEarth = frame
405 .getStaticTransformTo(earth.getBodyFrame(), date.toAbsoluteDate())
406 .transformPosition(position);
407
408 return getDensity(sunInEarth, posInEarth);
409 }
410
411 /** Get the height above the Earth for the given position.
412 * <p>
413 * The height computation is an approximation valid for the considered atmosphere.
414 * </p>
415 * @param position current position in Earth frame
416 * @return height (m)
417 */
418 private double getHeight(final Vector3D position) {
419 final double a = earth.getEquatorialRadius();
420 final double f = earth.getFlattening();
421 final double e2 = f * (2. - f);
422 final double r = position.getNorm();
423 final double sl = position.getZ() / r;
424 final double cl2 = 1. - sl * sl;
425 final double coef = FastMath.sqrt((1. - e2) / (1. - e2 * cl2));
426
427 return r - a * coef;
428 }
429
430 /** Get the height above the Earth for the given position.
431 * <p>
432 * The height computation is an approximation valid for the considered atmosphere.
433 * </p>
434 * @param position current position in Earth frame
435 * @param <T> instance of CalculusFieldElement<T>
436 * @return height (m)
437 */
438 private <T extends CalculusFieldElement<T>> T getHeight(final FieldVector3D<T> position) {
439 final double a = earth.getEquatorialRadius();
440 final double f = earth.getFlattening();
441 final double e2 = f * (2. - f);
442 final T r = position.getNorm();
443 final T sl = position.getZ().divide(r);
444 final T cl2 = sl.square().negate().add(1.);
445 final T coef = cl2.multiply(-e2).add(1.).reciprocal().multiply(1. - e2).sqrt();
446
447 return r.subtract(coef.multiply(a));
448 }
449
450 }