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