View Javadoc
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 fr.cs.examples.propagation;
18  
19  import java.io.BufferedWriter;
20  import java.io.File;
21  import java.io.FileInputStream;
22  import java.io.FileWriter;
23  import java.io.IOException;
24  import java.text.ParseException;
25  import java.util.ArrayList;
26  import java.util.Formatter;
27  import java.util.List;
28  import java.util.Locale;
29  import java.util.NoSuchElementException;
30  
31  import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
32  import org.apache.commons.math3.ode.AbstractIntegrator;
33  import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
34  import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator;
35  import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
36  import org.apache.commons.math3.util.FastMath;
37  import org.apache.commons.math3.util.MathUtils;
38  import org.orekit.Utils;
39  import org.orekit.bodies.CelestialBodyFactory;
40  import org.orekit.bodies.OneAxisEllipsoid;
41  import org.orekit.errors.OrekitException;
42  import org.orekit.forces.SphericalSpacecraft;
43  import org.orekit.forces.drag.Atmosphere;
44  import org.orekit.forces.drag.DragForce;
45  import org.orekit.forces.drag.HarrisPriester;
46  import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
47  import org.orekit.forces.gravity.ThirdBodyAttraction;
48  import org.orekit.forces.gravity.potential.GravityFieldFactory;
49  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
50  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
51  import org.orekit.forces.radiation.SolarRadiationPressure;
52  import org.orekit.frames.Frame;
53  import org.orekit.frames.FramesFactory;
54  import org.orekit.orbits.CartesianOrbit;
55  import org.orekit.orbits.CircularOrbit;
56  import org.orekit.orbits.EquinoctialOrbit;
57  import org.orekit.orbits.KeplerianOrbit;
58  import org.orekit.orbits.Orbit;
59  import org.orekit.orbits.PositionAngle;
60  import org.orekit.propagation.SpacecraftState;
61  import org.orekit.propagation.numerical.NumericalPropagator;
62  import org.orekit.propagation.sampling.OrekitFixedStepHandler;
63  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
64  import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
65  import org.orekit.propagation.semianalytical.dsst.forces.DSSTCentralBody;
66  import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
67  import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
68  import org.orekit.time.AbsoluteDate;
69  import org.orekit.time.TimeScale;
70  import org.orekit.time.TimeScalesFactory;
71  import org.orekit.utils.Constants;
72  import org.orekit.utils.PVCoordinates;
73  
74  import fr.cs.examples.KeyValueFileParser;
75  
76  /** Orekit tutorial for semi-analytical extrapolation using the DSST.
77   *  <p>
78   *  The parameters are read from the input file dsst-propagation.in located in the user's
79   *  home directory (see commented example at src/tutorial/ressources/dsst-propagation.in).
80   *  The results are written to the ouput file dsst-propagation.out in the same directory.
81   *  </p>
82   *  <p>
83   *  Comparison between the DSST propagator and the numerical propagator can be optionally
84   *  performed. Numerical results are  written to the ouput file numerical-propagation.out.
85   *  </p>
86   *
87   *  @author Romain Di Costanzo
88   *  @author Pascal Parraud
89   */
90  public class DSSTPropagation {
91  
92      /** Program entry point.
93       * @param args program arguments
94       */
95      public static void main(String[] args) {
96          try {
97  
98              // configure Orekit data acces
99              Utils.setDataRoot("tutorial-orekit-data");
100 
101             // input/output (in user's home directory)
102             File input  = new File(new File(System.getProperty("user.home")), "dsst-propagation.in");
103             File output = new File(input.getParentFile(), "dsst-propagation.out");
104 
105             new DSSTPropagation().run(input, output);
106 
107         } catch (IOException ioe) {
108             System.err.println(ioe.getLocalizedMessage());
109             System.exit(1);
110         } catch (IllegalArgumentException iae) {
111             System.err.println(iae.getLocalizedMessage());
112             System.exit(1);
113         } catch (OrekitException oe) {
114             System.err.println(oe.getLocalizedMessage());
115             System.exit(1);
116         } catch (ParseException pe) {
117             System.err.println(pe.getLocalizedMessage());
118             System.exit(1);
119         }
120     }
121 
122     /** Input parameter keys. */
123     private static enum ParameterKey {
124         ORBIT_DATE,
125         ORBIT_CIRCULAR_A,
126         ORBIT_CIRCULAR_EX,
127         ORBIT_CIRCULAR_EY,
128         ORBIT_CIRCULAR_I,
129         ORBIT_CIRCULAR_RAAN,
130         ORBIT_CIRCULAR_ALPHA,
131         ORBIT_EQUINOCTIAL_A,
132         ORBIT_EQUINOCTIAL_EX,
133         ORBIT_EQUINOCTIAL_EY,
134         ORBIT_EQUINOCTIAL_HX,
135         ORBIT_EQUINOCTIAL_HY,
136         ORBIT_EQUINOCTIAL_LAMBDA,
137         ORBIT_KEPLERIAN_A,
138         ORBIT_KEPLERIAN_E,
139         ORBIT_KEPLERIAN_I,
140         ORBIT_KEPLERIAN_PA,
141         ORBIT_KEPLERIAN_RAAN,
142         ORBIT_KEPLERIAN_ANOMALY,
143         ORBIT_ANGLE_TYPE,
144         ORBIT_CARTESIAN_PX,
145         ORBIT_CARTESIAN_PY,
146         ORBIT_CARTESIAN_PZ,
147         ORBIT_CARTESIAN_VX,
148         ORBIT_CARTESIAN_VY,
149         ORBIT_CARTESIAN_VZ,
150         ORBIT_IS_OSCULATING,
151         START_DATE,
152         DURATION,
153         DURATION_IN_DAYS,
154         OUTPUT_STEP,
155         FIXED_INTEGRATION_STEP,
156         NUMERICAL_COMPARISON,
157         CENTRAL_BODY_ORDER,
158         CENTRAL_BODY_DEGREE,
159         THIRD_BODY_MOON,
160         THIRD_BODY_SUN,
161         MASS,
162         DRAG,
163         DRAG_CD,
164         DRAG_SF,
165         SOLAR_RADIATION_PRESSURE,
166         SOLAR_RADIATION_PRESSURE_CR,
167         SOLAR_RADIATION_PRESSURE_SF;
168     }
169 
170     private void run(final File input, final File output)
171             throws IOException, IllegalArgumentException, OrekitException, ParseException {
172 
173         // read input parameters
174         KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
175         parser.parseInput(new FileInputStream(input));
176 
177         // check mandatory input parameters
178         if (!parser.containsKey(ParameterKey.ORBIT_DATE)) {
179             throw new IOException("Orbit date is not defined.");
180         }
181         if (!parser.containsKey(ParameterKey.DURATION) && !parser.containsKey(ParameterKey.DURATION_IN_DAYS)) {
182             throw new IOException("Propagation duration is not defined.");
183         }
184 
185         // All dates in UTC
186         final TimeScale utc = TimeScalesFactory.getUTC();
187 
188         final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE);
189         final int order  = FastMath.min(degree, parser.getInt(ParameterKey.CENTRAL_BODY_ORDER));
190 
191         // Potential coefficients providers
192         final UnnormalizedSphericalHarmonicsProvider unnormalized =
193                 GravityFieldFactory.getConstantUnnormalizedProvider(degree, order);
194         final NormalizedSphericalHarmonicsProvider normalized =
195                 GravityFieldFactory.getConstantNormalizedProvider(degree, order);
196 
197         // Central body attraction coefficient (m³/s²)
198         final double mu = unnormalized.getMu();
199 
200         // Earth frame definition (for faster computation)
201         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
202 
203         // Orbit definition (inertial frame is EME2000)
204         final Orbit orbit = createOrbit(parser, FramesFactory.getEME2000(), utc, mu);
205 
206         // DSST propagator definition
207         double mass = 1000.0;
208         if (parser.containsKey(ParameterKey.MASS)) {
209             mass = parser.getDouble(ParameterKey.MASS);
210         }
211         Boolean isOsculating = false;
212         if (parser.containsKey(ParameterKey.ORBIT_IS_OSCULATING)) {
213             isOsculating = parser.getBoolean(ParameterKey.ORBIT_IS_OSCULATING);
214         }
215         double fixedStepSize = -1.;
216         if (parser.containsKey(ParameterKey.FIXED_INTEGRATION_STEP)) {
217             fixedStepSize = parser.getDouble(ParameterKey.FIXED_INTEGRATION_STEP);
218         }
219         final DSSTPropagator dsstProp = createDSSTProp(orbit, mass, isOsculating, fixedStepSize);
220 
221         // Set Force models
222         setForceModel(parser, unnormalized, earthFrame, dsstProp);
223 
224         // Simulation properties
225         AbsoluteDate start;
226         if (parser.containsKey(ParameterKey.START_DATE)) {
227             start = parser.getDate(ParameterKey.START_DATE, utc);
228         } else {
229             start = parser.getDate(ParameterKey.ORBIT_DATE, utc);
230         }
231         double duration = 0.;
232         if (parser.containsKey(ParameterKey.DURATION)) {
233             duration = parser.getDouble(ParameterKey.DURATION);
234         }
235         if (parser.containsKey(ParameterKey.DURATION_IN_DAYS)) {
236             duration = parser.getDouble(ParameterKey.DURATION_IN_DAYS) * Constants.JULIAN_DAY;
237         }
238         double outStep = parser.getDouble(ParameterKey.OUTPUT_STEP);
239 
240         // Add orbit handler
241         OrbitHandler dsstHandler = new OrbitHandler();
242         dsstProp.setMasterMode(outStep, dsstHandler);
243 
244         // DSST Propagation
245         final double dsstOn = System.currentTimeMillis();
246         dsstProp.propagate(start, start.shiftedBy(duration));
247         final double dsstOff = System.currentTimeMillis();
248         System.out.println("DSST execution time: " + (dsstOff - dsstOn) / 1000.);
249 
250         // Print results
251         printOutput(output, dsstHandler, start);
252         System.out.println("DSST results saved as file " + output);
253 
254         // Check if we want to compare numerical to DSST propagator (default is false)
255         if (parser.containsKey(ParameterKey.NUMERICAL_COMPARISON)
256                 && parser.getBoolean(ParameterKey.NUMERICAL_COMPARISON)) {
257             
258             if ( !isOsculating ) {
259                 System.out.println("\nWARNING:");
260                 System.out.println("The DSST propagator considers a mean orbit while the numerical will consider an osculating one.");
261                 System.out.println("The comparison will be meaningless.\n");
262             }
263 
264             // output (in user's home directory)
265             File output_num = new File(input.getParentFile(), "numerical-propagation.out");
266 
267             // Numerical propagator definition
268             final NumericalPropagator numProp = createNumProp(orbit, mass);
269             
270             // Set Force models
271             setForceModel(parser, normalized, earthFrame, numProp);
272 
273             // Add orbit handler
274             OrbitHandler numHandler = new OrbitHandler();
275             numProp.setMasterMode(outStep, numHandler);
276 
277             // Numerical Propagation
278             final double numOn = System.currentTimeMillis();
279             numProp.propagate(start, start.shiftedBy(duration));
280             final double numOff = System.currentTimeMillis();
281             System.out.println("Numerical execution time: " + (numOff - numOn) / 1000.);
282             
283             // Print results
284             printOutput(output_num, numHandler, start);
285             System.out.println("Numerical results saved as file " + output_num);
286         }
287 
288     }
289 
290     /** Create an orbit from input parameters
291      * @param parser input file parser
292      * @param frame  inertial frame
293      * @param scale  time scale
294      * @param mu     central attraction coefficient
295      * @throws NoSuchElementException if input parameters are mising
296      * @throws IOException if input parameters are invalid
297      */
298     private Orbit createOrbit(final KeyValueFileParser<ParameterKey> parser,
299                               final Frame frame, final TimeScale scale, final double mu)
300         throws NoSuchElementException, IOException {
301 
302         // Orbit definition
303         Orbit orbit;
304         PositionAngle angleType = PositionAngle.MEAN;
305         if (parser.containsKey(ParameterKey.ORBIT_ANGLE_TYPE)) {
306             angleType = PositionAngle.valueOf(parser.getString(ParameterKey.ORBIT_ANGLE_TYPE).toUpperCase());
307         }
308         if (parser.containsKey(ParameterKey.ORBIT_KEPLERIAN_A)) {
309             orbit = new KeplerianOrbit(parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_A) * 1000.,
310                                        parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_E),
311                                        parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_I),
312                                        parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_PA),
313                                        parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_RAAN),
314                                        parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_ANOMALY),
315                                        angleType,
316                                        frame,
317                                        parser.getDate(ParameterKey.ORBIT_DATE, scale),
318                                        mu
319                                       );
320         } else if (parser.containsKey(ParameterKey.ORBIT_EQUINOCTIAL_A)) {
321             orbit = new EquinoctialOrbit(parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_A) * 1000.,
322                                          parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EX),
323                                          parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EY),
324                                          parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HX),
325                                          parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HY),
326                                          parser.getAngle(ParameterKey.ORBIT_EQUINOCTIAL_LAMBDA),
327                                          angleType,
328                                          frame,
329                                          parser.getDate(ParameterKey.ORBIT_DATE, scale),
330                                          mu
331                                         );
332         } else if (parser.containsKey(ParameterKey.ORBIT_CIRCULAR_A)) {
333             orbit = new CircularOrbit(parser.getDouble(ParameterKey.ORBIT_CIRCULAR_A) * 1000.,
334                                       parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EX),
335                                       parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EY),
336                                       parser.getAngle(ParameterKey.ORBIT_CIRCULAR_I),
337                                       parser.getAngle(ParameterKey.ORBIT_CIRCULAR_RAAN),
338                                       parser.getAngle(ParameterKey.ORBIT_CIRCULAR_ALPHA),
339                                       angleType,
340                                       frame,
341                                       parser.getDate(ParameterKey.ORBIT_DATE, scale),
342                                       mu
343                                      );
344         } else if (parser.containsKey(ParameterKey.ORBIT_CARTESIAN_PX)) {
345             final double[] pos = {parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PX) * 1000.,
346                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PY) * 1000.,
347                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PZ) * 1000.};
348             final double[] vel = {parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VX) * 1000.,
349                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VY) * 1000.,
350                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VZ) * 1000.};
351             orbit = new CartesianOrbit(new PVCoordinates(new Vector3D(pos), new Vector3D(vel)),
352                                        frame,
353                                        parser.getDate(ParameterKey.ORBIT_DATE, scale),
354                                        mu
355                                       );
356         } else {
357             throw new IOException("Orbit definition is incomplete.");
358         }
359 
360         return orbit;
361 
362     }
363     /** Set up the DSST Propagator
364      *
365      *  @param orbit initial orbit
366      *  @param mass S/C mass (kg)
367      *  @param isOsculating if orbital elements are osculating
368      *  @param fixedStepSize step size for fixed step integrator (s)
369      *  @throws OrekitException
370      */
371     private DSSTPropagator createDSSTProp(final Orbit orbit,
372                                           final double mass,
373                                           final boolean isOsculating,
374                                           final double fixedStepSize) throws OrekitException {
375         AbstractIntegrator integrator;
376         if (fixedStepSize > 0.) {
377             integrator = new ClassicalRungeKuttaIntegrator(fixedStepSize);
378         } else {
379             final double minStep = orbit.getKeplerianPeriod();
380             final double maxStep = minStep * 100.;
381             final double[][] tol = DSSTPropagator.tolerances(1.0, orbit);
382             integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
383             ((AdaptiveStepsizeIntegrator) integrator).setInitialStepSize(10. * minStep);
384         }
385 
386         DSSTPropagator dsstProp = new DSSTPropagator(integrator);
387         dsstProp.setInitialState(new SpacecraftState(orbit, mass), isOsculating);
388 
389         return dsstProp;
390     }
391 
392     /** Create the numerical propagator
393      *
394      *  @param orbit initial orbit
395      *  @param mass S/C mass (kg)
396      *  @throws OrekitException 
397      */
398     private NumericalPropagator createNumProp(final Orbit orbit, final double mass) throws OrekitException {
399         final double[][] tol = NumericalPropagator.tolerances(1.0, orbit, orbit.getType());
400         final double minStep = 1.e-3;
401         final double maxStep = 1.e+3;
402         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
403         integrator.setInitialStepSize(100.);
404 
405         NumericalPropagator numProp = new NumericalPropagator(integrator);
406         numProp.setInitialState(new SpacecraftState(orbit, mass));
407 
408         return numProp;
409     }
410 
411     /** Set DSST propagator force models
412      *
413      *  @param parser input file parser
414      *  @param unnormalized spherical harmonics provider
415      *  @param earthFrame Earth rotating frame
416      *  @param dsstProp DSST propagator
417      *  @throws IOException
418      *  @throws OrekitException
419      */
420     private void setForceModel(final KeyValueFileParser<ParameterKey> parser,
421                                final UnnormalizedSphericalHarmonicsProvider unnormalized,
422                                final Frame earthFrame,
423                                final DSSTPropagator dsstProp) throws IOException, OrekitException {
424 
425         final double ae = unnormalized.getAe();
426         
427         final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE);
428         final int order  = parser.getInt(ParameterKey.CENTRAL_BODY_ORDER);
429 
430         if (order > degree) {
431             throw new IOException("Potential order cannot be higher than potential degree");
432         }
433 
434         // Central Body Force Model with un-normalized coefficients
435         dsstProp.addForceModel(new DSSTCentralBody(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, unnormalized));
436 
437         // 3rd body (SUN)
438         if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
439             dsstProp.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getSun()));
440         }
441 
442         // 3rd body (MOON)
443         if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
444             dsstProp.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getMoon()));
445         }
446 
447         // Drag
448         if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
449             final OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame);
450             final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
451             dsstProp.addForceModel(new DSSTAtmosphericDrag(atm,
452                                                            parser.getDouble(ParameterKey.DRAG_CD),
453                                                            parser.getDouble(ParameterKey.DRAG_SF)));
454         }
455 
456         // Solar Radiation Pressure
457         if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) {
458             dsstProp.addForceModel(new DSSTSolarRadiationPressure(parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR),
459                                                                   parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_SF),
460                                                                   CelestialBodyFactory.getSun(), ae));
461         }
462     }
463 
464     /** Set numerical propagator force models
465      *
466      *  @param parser  input file parser
467      *  @param normalized spherical harmonics provider
468      *  @param earthFrame Earth rotating frame
469      *  @param numProp numerical propagator
470      *  @throws IOException
471      *  @throws OrekitException
472      */
473     private void setForceModel(final KeyValueFileParser<ParameterKey> parser,
474                                final NormalizedSphericalHarmonicsProvider normalized,
475                                final Frame earthFrame,
476                                final NumericalPropagator numProp) throws IOException, OrekitException {
477 
478         final double ae = normalized.getAe();
479         
480         final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE);
481         final int order  = parser.getInt(ParameterKey.CENTRAL_BODY_ORDER);
482 
483         if (order > degree) {
484             throw new IOException("Potential order cannot be higher than potential degree");
485         }
486 
487         // Central Body (normalized coefficients)
488         numProp.addForceModel(new HolmesFeatherstoneAttractionModel(earthFrame, normalized));
489 
490         // 3rd body (SUN)
491         if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
492             numProp.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
493         }
494 
495         // 3rd body (MOON)
496         if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
497             numProp.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
498         }
499 
500         // Drag
501         if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
502             final OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame);
503             final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
504             final SphericalSpacecraft ssc = new SphericalSpacecraft(parser.getDouble(ParameterKey.DRAG_SF),
505                                                                     parser.getDouble(ParameterKey.DRAG_CD),
506                                                                     0., 0.);
507             numProp.addForceModel(new DragForce(atm, ssc));
508         }
509 
510         // Solar Radiation Pressure
511         if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) {
512             final double cR = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR);
513             // cR being the DSST SRP coef and assuming a spherical spacecraft, the conversion is:
514             // cR = 1 + (1 - kA) * (1 - kR) * 4 / 9
515             // with kA arbitrary sets to 0
516             final double kR = 3.25 - 2.25 * cR;
517             final SphericalSpacecraft ssc = new SphericalSpacecraft(parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_SF),
518                                                                     0., 0., kR);
519             numProp.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), ae, ssc));
520         }
521     }
522 
523     /** Print the results in the output file
524      *
525      *  @param handler orbit handler
526      *  @param output output file
527      *  @param sart start date of propagation
528      *  @throws OrekitException 
529      *  @throws IOException 
530      */
531     private void printOutput(final File output,
532                              final OrbitHandler handler,
533                              final AbsoluteDate start) throws OrekitException, IOException {
534         // Output format:
535         // time_from_start, a, e, i, raan, pa, aM, h, k, p, q, lM, px, py, pz, vx, vy, vz
536         final String format = new String(" %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e %24.16e");
537         final BufferedWriter buffer = new BufferedWriter(new FileWriter(output));
538         buffer.write("##   time_from_start(s)            a(km)                      e                      i(deg)         ");
539         buffer.write("         raan(deg)                pa(deg)              mean_anomaly(deg)              ey/h          ");
540         buffer.write("           ex/k                    hy/p                     hx/q             mean_longitude_arg(deg)");
541         buffer.write("       Xposition(km)           Yposition(km)             Zposition(km)           Xvelocity(km/s)    ");
542         buffer.write("      Yvelocity(km/s)         Zvelocity(km/s)");
543         buffer.newLine();
544         for (Orbit o : handler.getOrbits()) {
545             final Formatter f = new Formatter(new StringBuilder(), Locale.ENGLISH);
546             // Time from start (s)
547             final double time = o.getDate().durationFrom(start);
548             // Semi-major axis (km)
549             final double a = o.getA() / 1000.;
550             // Keplerian elements
551             // Eccentricity
552             final double e = o.getE();
553             // Inclination (degrees)
554             final double i = Math.toDegrees(MathUtils.normalizeAngle(o.getI(), FastMath.PI));
555             // Right Ascension of Ascending Node (degrees)
556             KeplerianOrbit ko = new KeplerianOrbit(o);
557             final double ra = Math.toDegrees(MathUtils.normalizeAngle(ko.getRightAscensionOfAscendingNode(), FastMath.PI));
558             // Perigee Argument (degrees)
559             final double pa = Math.toDegrees(MathUtils.normalizeAngle(ko.getPerigeeArgument(), FastMath.PI));
560             // Mean Anomaly (degrees)
561             final double am = Math.toDegrees(MathUtils.normalizeAngle(ko.getAnomaly(PositionAngle.MEAN), FastMath.PI));
562             // Equinoctial elements
563             // ey/h component of eccentricity vector
564             final double h = o.getEquinoctialEy();
565             // ex/k component of eccentricity vector
566             final double k = o.getEquinoctialEx();
567             // hy/p component of inclination vector
568             final double p = o.getHy();
569             // hx/q component of inclination vector
570             final double q = o.getHx();
571             // Mean Longitude Argument (degrees)
572             final double lm = Math.toDegrees(MathUtils.normalizeAngle(o.getLM(), FastMath.PI));
573             // Cartesian elements
574             // Position along X in inertial frame (km)
575             final double px = o.getPVCoordinates().getPosition().getX() / 1000.;
576             // Position along Y in inertial frame (km)
577             final double py = o.getPVCoordinates().getPosition().getY() / 1000.;
578             // Position along Z in inertial frame (km)
579             final double pz = o.getPVCoordinates().getPosition().getZ() / 1000.;
580             // Velocity along X in inertial frame (km/s)
581             final double vx = o.getPVCoordinates().getVelocity().getX() / 1000.;
582             // Velocity along Y in inertial frame (km/s)
583             final double vy = o.getPVCoordinates().getVelocity().getY() / 1000.;
584             // Velocity along Z in inertial frame (km/s)
585             final double vz = o.getPVCoordinates().getVelocity().getZ() / 1000.;
586             buffer.write(f.format(format, time, a, e, i, ra, pa, am, h, k, p, q, lm, px, py, pz, vx, vy, vz).toString());
587             buffer.newLine();
588             f.close();
589         }
590         buffer.close();
591     }
592 
593     /** Specialized step handler catching the orbit at each step. */
594     private static class OrbitHandler implements OrekitFixedStepHandler {
595 
596         /** List of orbits. */
597         private final List<Orbit> orbits;
598 
599         private OrbitHandler() {
600             // initialise an empty list of orbit
601             orbits = new ArrayList<Orbit>();
602         }
603 
604         /** {@inheritDoc} */
605         public void init(final SpacecraftState s0, final AbsoluteDate t) {
606         }
607 
608         /** {@inheritDoc} */
609         public void handleStep(SpacecraftState currentState, boolean isLast) {
610             // fill in the list with the orbit from the current step
611             orbits.add(currentState.getOrbit());
612         }
613 
614         /** Get the list of propagated orbits.
615          * @return orbits
616          */
617         public List<Orbit> getOrbits() {
618             return orbits;
619         }
620     }
621 
622 }