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