1 /* Copyright 2002-2015 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (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.forces;
18
19 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
20 import org.apache.commons.math3.util.FastMath;
21 import org.apache.commons.math3.util.MathUtils;
22 import org.apache.commons.math3.util.Precision;
23 import org.orekit.errors.OrekitException;
24 import org.orekit.forces.SphericalSpacecraft;
25 import org.orekit.forces.radiation.RadiationSensitive;
26 import org.orekit.forces.radiation.SolarRadiationPressure;
27 import org.orekit.propagation.SpacecraftState;
28 import org.orekit.propagation.events.EventDetector;
29 import org.orekit.utils.PVCoordinatesProvider;
30
31 /** Solar radiation pressure contribution to the
32 * {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
33 * <p>
34 * The solar radiation pressure acceleration is computed through the
35 * acceleration model of
36 * {@link org.orekit.forces.radiation.SolarRadiationPressure SolarRadiationPressure}.
37 * </p>
38 *
39 * @author Pascal Parraud
40 */
41 public class DSSTSolarRadiationPressure extends AbstractGaussianContribution {
42
43 /** Reference distance for the solar radiation pressure (m). */
44 private static final double D_REF = 149597870000.0;
45
46 /** Reference solar radiation pressure at D_REF (N/m²). */
47 private static final double P_REF = 4.56e-6;
48
49 /** Threshold for the choice of the Gauss quadrature order. */
50 private static final double GAUSS_THRESHOLD = 1.0e-15;
51
52 /** Threshold for shadow equation. */
53 private static final double S_ZERO = 1.0e-6;
54
55 /** Sun model. */
56 private final PVCoordinatesProvider sun;
57
58 /** Central Body radius. */
59 private final double ae;
60
61 /** Spacecraft model for radiation acceleration computation. */
62 private final RadiationSensitive spacecraft;
63
64
65 /**
66 * Simple constructor with default reference values and spherical spacecraft.
67 * <p>
68 * When this constructor is used, the reference values are:
69 * </p>
70 * <ul>
71 * <li>d<sub>ref</sub> = 149597870000.0 m</li>
72 * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
73 * </ul>
74 * <p>
75 * The spacecraft has a spherical shape.
76 * </p>
77 *
78 * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
79 * @param area cross sectionnal area of satellite
80 * @param sun Sun model
81 * @param equatorialRadius central body equatorial radius (for shadow computation)
82 */
83 public DSSTSolarRadiationPressure(final double cr, final double area,
84 final PVCoordinatesProvider sun,
85 final double equatorialRadius) {
86 this(D_REF, P_REF, cr, area, sun, equatorialRadius);
87 }
88
89 /**
90 * Simple constructor with default reference values, but custom spacecraft.
91 * <p>
92 * When this constructor is used, the reference values are:
93 * </p>
94 * <ul>
95 * <li>d<sub>ref</sub> = 149597870000.0 m</li>
96 * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
97 * </ul>
98 *
99 * @param sun Sun model
100 * @param equatorialRadius central body equatorial radius (for shadow computation)
101 * @param spacecraft spacecraft model
102 */
103 public DSSTSolarRadiationPressure(final PVCoordinatesProvider sun,
104 final double equatorialRadius,
105 final RadiationSensitive spacecraft) {
106 this(D_REF, P_REF, sun, equatorialRadius, spacecraft);
107 }
108
109 /**
110 * Constructor with customizable reference values but spherical spacecraft.
111 * <p>
112 * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
113 * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
114 * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
115 * N/m² solar radiation pressure.
116 * </p>
117 *
118 * @param dRef reference distance for the solar radiation pressure (m)
119 * @param pRef reference solar radiation pressure at dRef (N/m²)
120 * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
121 * @param area cross sectionnal area of satellite
122 * @param sun Sun model
123 * @param equatorialRadius central body equatrial radius (for shadow computation)
124 */
125 public DSSTSolarRadiationPressure(final double dRef, final double pRef,
126 final double cr, final double area,
127 final PVCoordinatesProvider sun,
128 final double equatorialRadius) {
129
130 // cR being the DSST SRP coef and assuming a spherical spacecraft,
131 // the conversion is:
132 // cR = 1 + (1 - kA) * (1 - kR) * 4 / 9
133 // with kA arbitrary sets to 0
134 this(dRef, pRef, sun, equatorialRadius, new SphericalSpacecraft(
135 area, 0.0, 0.0, 3.25 - 2.25 * cr));
136 }
137
138 /**
139 * Complete constructor.
140 * <p>
141 * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
142 * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
143 * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
144 * N/m² solar radiation pressure.
145 * </p>
146 *
147 * @param dRef reference distance for the solar radiation pressure (m)
148 * @param pRef reference solar radiation pressure at dRef (N/m²)
149 * @param sun Sun model
150 * @param equatorialRadius central body equatrial radius (for shadow computation)
151 * @param spacecraft spacecraft model
152 */
153 public DSSTSolarRadiationPressure(final double dRef, final double pRef,
154 final PVCoordinatesProvider sun, final double equatorialRadius,
155 final RadiationSensitive spacecraft) {
156
157 //Call to the constructor from superclass using the numerical SRP model as ForceModel
158 super(GAUSS_THRESHOLD, new SolarRadiationPressure(
159 dRef, pRef, sun, equatorialRadius, spacecraft));
160
161 this.sun = sun;
162 this.ae = equatorialRadius;
163 this.spacecraft = spacecraft;
164 }
165
166 /** Get spacecraft shape.
167 * @return the spacecraft shape.
168 */
169 public RadiationSensitive getSpacecraft() {
170 return spacecraft;
171 }
172
173 /** {@inheritDoc} */
174 public EventDetector[] getEventsDetectors() {
175 // TODO: add eclipses handling
176 return null;
177 }
178
179 /** {@inheritDoc} */
180 protected double[] getLLimits(final SpacecraftState state) throws OrekitException {
181 // Default bounds without shadow [-PI, PI]
182 final double[] ll = {-FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0),
183 FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0)};
184
185 // Direction cosines of the Sun in the equinoctial frame
186 final Vector3D sunDir = sun.getPVCoordinates(state.getDate(), state.getFrame()).getPosition().normalize();
187 final double alpha = sunDir.dotProduct(f);
188 final double beta = sunDir.dotProduct(g);
189 final double gamma = sunDir.dotProduct(w);
190
191 // Compute limits only if the perigee is close enough from the central body to be in the shadow
192 if (FastMath.abs(gamma * a * (1. - ecc)) < ae) {
193
194 // Compute the coefficients of the quartic equation in cos(L) 3.5-(2)
195 final double bet2 = beta * beta;
196 final double h2 = h * h;
197 final double k2 = k * k;
198 final double m = ae / (a * B);
199 final double m2 = m * m;
200 final double m4 = m2 * m2;
201 final double bb = alpha * beta + m2 * h * k;
202 final double b2 = bb * bb;
203 final double cc = alpha * alpha - bet2 + m2 * (k2 - h2);
204 final double dd = 1. - bet2 - m2 * (1. + h2);
205 final double[] a = new double[5];
206 a[0] = 4. * b2 + cc * cc;
207 a[1] = 8. * bb * m2 * h + 4. * cc * m2 * k;
208 a[2] = -4. * b2 + 4. * m4 * h2 - 2. * cc * dd + 4. * m4 * k2;
209 a[3] = -8. * bb * m2 * h - 4. * dd * m2 * k;
210 a[4] = -4. * m4 * h2 + dd * dd;
211 // Compute the real roots of the quartic equation 3.5-2
212 final double[] roots = new double[4];
213 final int nbRoots = realQuarticRoots(a, roots);
214 if (nbRoots > 0) {
215 // Check for consistency
216 boolean entryFound = false;
217 boolean exitFound = false;
218 // Eliminate spurious roots
219 for (int i = 0; i < nbRoots; i++) {
220 final double cosL = roots[i];
221 final double sL = FastMath.sqrt((1. - cosL) * (1. + cosL));
222 // Check both angles: L and -L
223 for (int j = -1; j <= 1; j += 2) {
224 final double sinL = j * sL;
225 final double cPhi = alpha * cosL + beta * sinL;
226 // Is the angle on the shadow side of the central body (eq. 3.5-3) ?
227 if (cPhi < 0.) {
228 final double range = 1. + k * cosL + h * sinL;
229 final double S = 1. - m2 * range * range - cPhi * cPhi;
230 // Is the shadow equation 3.5-1 satisfied ?
231 if (FastMath.abs(S) < S_ZERO) {
232 // Is this the entry or exit angle ?
233 final double dSdL = m2 * range * (k * sinL - h * cosL) + cPhi * (alpha * sinL - beta * cosL);
234 if (dSdL > 0.) {
235 // Exit from shadow: 3.5-4
236 exitFound = true;
237 ll[0] = FastMath.atan2(sinL, cosL);
238 } else {
239 // Entry into shadow: 3.5-5
240 entryFound = true;
241 ll[1] = FastMath.atan2(sinL, cosL);
242 }
243 }
244 }
245 }
246 }
247 // Must be one entry and one exit or none
248 if (!(entryFound == exitFound)) {
249 // entry or exit found but not both ! In this case, consider there is no eclipse...
250 ll[0] = -FastMath.PI;
251 ll[1] = FastMath.PI;
252 }
253 // Quadrature between L at exit and L at entry so Lexit must be lower than Lentry
254 if (ll[0] > ll[1]) {
255 // Keep the angles between [-2PI, 2PI]
256 if (ll[1] < 0.) {
257 ll[1] += 2. * FastMath.PI;
258 } else {
259 ll[0] -= 2. * FastMath.PI;
260 }
261 }
262 }
263 }
264 return ll;
265 }
266
267 /** Get the central body equatorial radius.
268 * @return central body equatorial radius (m)
269 */
270 public double getEquatorialRadius() {
271 return ae;
272 }
273
274 /**
275 * Compute the real roots of a quartic equation.
276 *
277 * <pre>
278 * a[0] * x⁴ + a[1] * x³ + a[2] * x² + a[3] * x + a[4] = 0.
279 * </pre>
280 *
281 * @param a the 5 coefficients
282 * @param y the real roots
283 * @return the number of real roots
284 */
285 private int realQuarticRoots(final double[] a, final double[] y) {
286 /* Treat the degenerate quartic as cubic */
287 if (Precision.equals(a[0], 0.)) {
288 final double[] aa = new double[a.length - 1];
289 System.arraycopy(a, 1, aa, 0, aa.length);
290 return realCubicRoots(aa, y);
291 }
292
293 // Transform coefficients
294 final double b = a[1] / a[0];
295 final double c = a[2] / a[0];
296 final double d = a[3] / a[0];
297 final double e = a[4] / a[0];
298 final double bh = b * 0.5;
299
300 // Solve resolvant cubic
301 final double[] z3 = new double[3];
302 final int i3 = realCubicRoots(new double[] {1.0, -c, b * d - 4. * e, e * (4. * c - b * b) - d * d}, z3);
303 if (i3 == 0) {
304 return 0;
305 }
306
307 // Largest real root of resolvant cubic
308 final double z = z3[0];
309
310 // Compute auxiliary quantities
311 final double zh = z * 0.5;
312 final double p = FastMath.max(z + bh * bh - c, 0.);
313 final double q = FastMath.max(zh * zh - e, 0.);
314 final double r = bh * z - d;
315 final double pp = FastMath.sqrt(p);
316 final double qq = FastMath.copySign(FastMath.sqrt(q), r);
317
318 // Solve quadratic factors of quartic equation
319 final double[] y1 = new double[2];
320 final int n1 = realQuadraticRoots(new double[] {1.0, bh - pp, zh - qq}, y1);
321 final double[] y2 = new double[2];
322 final int n2 = realQuadraticRoots(new double[] {1.0, bh + pp, zh + qq}, y2);
323
324 if (n1 == 2) {
325 if (n2 == 2) {
326 y[0] = y1[0];
327 y[1] = y1[1];
328 y[2] = y2[0];
329 y[3] = y2[1];
330 return 4;
331 } else {
332 y[0] = y1[0];
333 y[1] = y1[1];
334 return 2;
335 }
336 } else {
337 if (n2 == 2) {
338 y[0] = y2[0];
339 y[1] = y2[1];
340 return 2;
341 } else {
342 return 0;
343 }
344 }
345 }
346
347 /**
348 * Compute the real roots of a cubic equation.
349 *
350 * <pre>
351 * a[0] * x³ + a[1] * x² + a[2] * x + a[3] = 0.
352 * </pre>
353 *
354 * @param a the 4 coefficients
355 * @param y the real roots sorted in descending order
356 * @return the number of real roots
357 */
358 private int realCubicRoots(final double[] a, final double[] y) {
359 if (Precision.equals(a[0], 0.)) {
360 // Treat the degenerate cubic as quadratic
361 final double[] aa = new double[a.length - 1];
362 System.arraycopy(a, 1, aa, 0, aa.length);
363 return realQuadraticRoots(aa, y);
364 }
365
366 // Transform coefficients
367 final double b = -a[1] / (3. * a[0]);
368 final double c = a[2] / a[0];
369 final double d = a[3] / a[0];
370 final double b2 = b * b;
371 final double p = b2 - c / 3.;
372 final double q = b * (b2 - c * 0.5) - d * 0.5;
373
374 // Compute discriminant
375 final double disc = p * p * p - q * q;
376
377 if (disc < 0.) {
378 // One real root
379 final double alpha = q + FastMath.copySign(FastMath.sqrt(-disc), q);
380 final double cbrtAl = FastMath.cbrt(alpha);
381 final double cbrtBe = p / cbrtAl;
382
383 if (p < 0.) {
384 y[0] = b + 2. * q / (cbrtAl * cbrtAl + cbrtBe * cbrtBe - p);
385 } else if (p > 0.) {
386 y[0] = b + cbrtAl + cbrtBe;
387 } else {
388 y[0] = b + cbrtAl;
389 }
390
391 return 1;
392
393 } else if (disc > 0.) {
394 // Three distinct real roots
395 final double phi = FastMath.atan2(FastMath.sqrt(disc), q) / 3.;
396 final double sqP = 2.0 * FastMath.sqrt(p);
397
398 y[0] = b + sqP * FastMath.cos(phi);
399 y[1] = b - sqP * FastMath.cos(FastMath.PI / 3. + phi);
400 y[2] = b - sqP * FastMath.cos(FastMath.PI / 3. - phi);
401
402 return 3;
403
404 } else {
405 // One distinct and two equals real roots
406 final double cbrtQ = FastMath.cbrt(q);
407 final double root1 = b + 2. * cbrtQ;
408 final double root2 = b - cbrtQ;
409 if (q < 0.) {
410 y[0] = root2;
411 y[1] = root2;
412 y[2] = root1;
413 } else {
414 y[0] = root1;
415 y[1] = root2;
416 y[2] = root2;
417 }
418
419 return 3;
420 }
421 }
422
423 /**
424 * Compute the real roots of a quadratic equation.
425 *
426 * <pre>
427 * a[0] * x² + a[1] * x + a[2] = 0.
428 * </pre>
429 *
430 * @param a the 3 coefficients
431 * @param y the real roots sorted in descending order
432 * @return the number of real roots
433 */
434 private int realQuadraticRoots(final double[] a, final double[] y) {
435 if (Precision.equals(a[0], 0.)) {
436 // Degenerate quadratic
437 if (Precision.equals(a[1], 0.)) {
438 // Degenerate linear equation: no real roots
439 return 0;
440 }
441 // Linear equation: one real root
442 y[0] = -a[2] / a[1];
443 return 1;
444 }
445
446 // Transform coefficients
447 final double b = -0.5 * a[1] / a[0];
448 final double c = a[2] / a[0];
449
450 // Compute discriminant
451 final double d = b * b - c;
452
453 if (d < 0.) {
454 // No real roots
455 return 0;
456 } else if (d > 0.) {
457 // Two distinct real roots
458 final double y0 = b + FastMath.copySign(FastMath.sqrt(d), b);
459 final double y1 = c / y0;
460 y[0] = FastMath.max(y0, y1);
461 y[1] = FastMath.min(y0, y1);
462 return 2;
463 } else {
464 // Discriminant is zero: two equal real roots
465 y[0] = b;
466 y[1] = b;
467 return 2;
468 }
469 }
470 }