PosVelChebyshev.java

  1. /* Copyright 2002-2016 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.bodies;

  18. import java.io.Serializable;

  19. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.orekit.time.AbsoluteDate;
  21. import org.orekit.time.TimeScale;
  22. import org.orekit.time.TimeStamped;
  23. import org.orekit.utils.PVCoordinates;


  24. /** Position-Velocity model based on Chebyshev polynomials.
  25.  * <p>This class represent the most basic element of the piecewise ephemerides
  26.  * for solar system bodies like JPL DE 405 ephemerides.</p>
  27.  * @see JPLEphemeridesLoader
  28.  * @author Luc Maisonobe
  29.  */
  30. class PosVelChebyshev implements TimeStamped, Serializable {

  31.     /** Serializable UID. */
  32.     private static final long serialVersionUID = 20151023L;

  33.     /** Time scale in which the ephemeris is defined. */
  34.     private final TimeScale timeScale;

  35.     /** Start of the validity range of the instance. */
  36.     private final AbsoluteDate start;

  37.     /** Duration of validity range of the instance. */
  38.     private final double duration;

  39.     /** Chebyshev polynomials coefficients for the X component. */
  40.     private final double[] xCoeffs;

  41.     /** Chebyshev polynomials coefficients for the Y component. */
  42.     private final double[] yCoeffs;

  43.     /** Chebyshev polynomials coefficients for the Z component. */
  44.     private final double[] zCoeffs;

  45.     /** Simple constructor.
  46.      * @param start start of the validity range of the instance
  47.      * @param timeScale time scale in which the ephemeris is defined
  48.      * @param duration duration of the validity range of the instance
  49.      * @param xCoeffs Chebyshev polynomials coefficients for the X component
  50.      * (a reference to the array will be stored in the instance)
  51.      * @param yCoeffs Chebyshev polynomials coefficients for the Y component
  52.      * (a reference to the array will be stored in the instance)
  53.      * @param zCoeffs Chebyshev polynomials coefficients for the Z component
  54.      * (a reference to the array will be stored in the instance)
  55.      */
  56.     PosVelChebyshev(final AbsoluteDate start, final TimeScale timeScale, final double duration,
  57.                     final double[] xCoeffs, final double[] yCoeffs, final double[] zCoeffs) {
  58.         this.start     = start;
  59.         this.timeScale = timeScale;
  60.         this.duration  = duration;
  61.         this.xCoeffs   = xCoeffs;
  62.         this.yCoeffs   = yCoeffs;
  63.         this.zCoeffs   = zCoeffs;
  64.     }

  65.     /** {@inheritDoc} */
  66.     public AbsoluteDate getDate() {
  67.         return start;
  68.     }

  69.     /** Check if a date is in validity range.
  70.      * @param date date to check
  71.      * @return true if date is in validity range
  72.      */
  73.     public boolean inRange(final AbsoluteDate date) {
  74.         final double dt = date.offsetFrom(start, timeScale);
  75.         return (dt >= -0.001) && (dt <= duration + 0.001);
  76.     }

  77.     /** Get the position-velocity-acceleration at a specified date.
  78.      * @param date date at which position-velocity-acceleration is requested
  79.      * @return position-velocity-acceleration at specified date
  80.      */
  81.     public PVCoordinates getPositionVelocityAcceleration(final AbsoluteDate date) {

  82.         // normalize date
  83.         final double t = (2 * date.offsetFrom(start, timeScale) - duration) / duration;
  84.         final double twoT = 2 * t;

  85.         // initialize Chebyshev polynomials recursion
  86.         double pKm1 = 1;
  87.         double pK   = t;
  88.         double xP   = xCoeffs[0];
  89.         double yP   = yCoeffs[0];
  90.         double zP   = zCoeffs[0];

  91.         // initialize Chebyshev polynomials derivatives recursion
  92.         double qKm1 = 0;
  93.         double qK   = 1;
  94.         double xV   = 0;
  95.         double yV   = 0;
  96.         double zV   = 0;

  97.         // initialize Chebyshev polynomials second derivatives recursion
  98.         double rKm1 = 0;
  99.         double rK   = 0;
  100.         double xA   = 0;
  101.         double yA   = 0;
  102.         double zA   = 0;

  103.         // combine polynomials by applying coefficients
  104.         for (int k = 1; k < xCoeffs.length; ++k) {

  105.             // consider last computed polynomials on position
  106.             xP += xCoeffs[k] * pK;
  107.             yP += yCoeffs[k] * pK;
  108.             zP += zCoeffs[k] * pK;

  109.             // consider last computed polynomials on velocity
  110.             xV += xCoeffs[k] * qK;
  111.             yV += yCoeffs[k] * qK;
  112.             zV += zCoeffs[k] * qK;

  113.             // consider last computed polynomials on acceleration
  114.             xA += xCoeffs[k] * rK;
  115.             yA += yCoeffs[k] * rK;
  116.             zA += zCoeffs[k] * rK;

  117.             // compute next Chebyshev polynomial value
  118.             final double pKm2 = pKm1;
  119.             pKm1 = pK;
  120.             pK   = twoT * pKm1 - pKm2;

  121.             // compute next Chebyshev polynomial derivative
  122.             final double qKm2 = qKm1;
  123.             qKm1 = qK;
  124.             qK   = twoT * qKm1 + 2 * pKm1 - qKm2;

  125.             // compute next Chebyshev polynomial second derivative
  126.             final double rKm2 = rKm1;
  127.             rKm1 = rK;
  128.             rK   = twoT * rKm1 + 4 * qKm1 - rKm2;

  129.         }

  130.         final double vScale = 2 / duration;
  131.         final double aScale = vScale * vScale;
  132.         return new PVCoordinates(new Vector3D(xP, yP, zP),
  133.                                  new Vector3D(xV * vScale, yV * vScale, zV * vScale),
  134.                                  new Vector3D(xA * aScale, yA * aScale, zA * aScale));

  135.     }

  136. }