CunninghamAttractionModel.java

  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.forces.gravity;


  20. import java.util.Collections;

  22. import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
  23. import org.apache.commons.math3.geometry.euclidean.threed.FieldRotation;
  24. import org.apache.commons.math3.geometry.euclidean.threed.FieldVector3D;
  25. import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
  26. import org.apache.commons.math3.ode.AbstractParameterizable;
  27. import org.apache.commons.math3.util.FastMath;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.forces.ForceModel;
  31. import org.orekit.forces.gravity.potential.GravityFieldFactory;
  32. import org.orekit.forces.gravity.potential.TideSystem;
  33. import org.orekit.forces.gravity.potential.TideSystemProvider;
  34. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  35. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  36. import org.orekit.frames.Frame;
  37. import org.orekit.frames.Transform;
  38. import org.orekit.propagation.SpacecraftState;
  39. import org.orekit.propagation.events.EventDetector;
  40. import org.orekit.propagation.numerical.Jacobianizer;
  41. import org.orekit.propagation.numerical.ParameterConfiguration;
  42. import org.orekit.propagation.numerical.TimeDerivativesEquations;
  43. import org.orekit.time.AbsoluteDate;

  45. /** This class represents the gravitational field of a celestial body.
  46.  * <p>The algorithm implemented in this class has been designed by
  47.  * Leland E. Cunningham (Lockheed Missiles and Space Company, Sunnyvale
  48.  * and Astronomy Department University of California, Berkeley) in
  49.  * his 1969 paper: <em>On the computation of the spherical harmonic
  50.  * terms needed during the numerical integration of the orbital motion
  51.  * of an artificial satellite</em> (Celestial Mechanics 2, 1970).</p>
  52.  * <p>
  53.  * Note that this class can often not be used for high degrees (say
  54.  * above 90) as most modern gravity fields are provided as normalized
  55.  * coefficients and the un-normalization process to convert these
  56.  * coefficients underflows at degree and order 89. This class also
  57.  * does not provide analytical partial derivatives (it uses finite differences
  58.  * to compute them) and is much slower than {@link HolmesFeatherstoneAttractionModel}
  59.  * (even when no derivatives are computed). For all these reasons,
  60.  * it is recommended to use the {@link HolmesFeatherstoneAttractionModel
  61.  * Holmes-Featherstone model} rather than this class.
  62.  * </p>
  63.  * <p>
  64.  * As this class uses finite differences to compute derivatives, the steps for
  65.  * finite differences <strong>must</strong> be initialized by calling {@link
  66.  * #setSteps(double, double)} prior to use derivatives, otherwise an exception
  67.  * will be thrown by {@link #accelerationDerivatives(AbsoluteDate, Frame, FieldVector3D,
  68.  * FieldVector3D, FieldRotation, DerivativeStructure)} and by {@link
  69.  * #accelerationDerivatives(SpacecraftState, String)}.
  70.  * </p>
  71.  *
  72.  * @see HolmesFeatherstoneAttractionModel
  73.  * @author Fabien Maussion
  74.  * @author Luc Maisonobe
  75.  * @author V&eacute;ronique Pommier-Maurussane
  76.  */

  78. public class CunninghamAttractionModel
  79.     extends AbstractParameterizable implements ForceModel, TideSystemProvider {

  81.     /** Provider for the spherical harmonics. */
  82.     private final UnnormalizedSphericalHarmonicsProvider provider;

  84.     /** Central attraction coefficient. */
  85.     private double mu;

  87.     /** Rotating body. */
  88.     private final Frame bodyFrame;

  90.     /** Helper class computing acceleration derivatives. */
  91.     private Jacobianizer jacobianizer;

  93.     /** Creates a new instance.
  94.     * @param centralBodyFrame rotating body frame
  95.     * @param equatorialRadius reference equatorial radius of the potential
  96.     * @param mu central body attraction coefficient (m<sup>3</sup>/s<sup>2</sup>)
  97.     * @param C un-normalized coefficients array (cosine part)
  98.     * @param S un-normalized coefficients array (sine part)
  99.     * @deprecated since 6.0, replaced by {@link #CunninghamAttractionModel(Frame, UnnormalizedSphericalHarmonicsProvider)}
  100.     */
  101.     public CunninghamAttractionModel(final Frame centralBodyFrame,
  102.                                      final double equatorialRadius, final double mu,
  103.                                      final double[][] C, final double[][] S) {
  104.         this(centralBodyFrame, GravityFieldFactory.getUnnormalizedProvider(equatorialRadius, mu,
  105.                                                                            TideSystem.UNKNOWN, C, S));
  106.     }

  108.    /** Creates a new instance.
  109.    * @param centralBodyFrame rotating body frame
  110.    * @param provider provider for spherical harmonics
  111.    * @since 6.0
  112.    */
  113.     public CunninghamAttractionModel(final Frame centralBodyFrame,
  114.                                      final UnnormalizedSphericalHarmonicsProvider provider) {
  115.         super(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT);

  117.         this.provider     = provider;
  118.         this.mu           = provider.getMu();
  119.         this.bodyFrame    = centralBodyFrame;
  120.         this.jacobianizer = null;

  122.     }

  124.     /** Set the step for finite differences with respect to spacecraft position.
  125.      * @param hPosition step used for finite difference computation
  126.      * with respect to spacecraft position (m)
  127.      * @param hMu step used for finite difference computation
  128.      * with respect to central attraction coefficient (m<sup>3</sup>/s<sup>2</sup>)
  129.      */
  130.     public void setSteps(final double hPosition, final double hMu) {
  131.         final ParameterConfiguration muConfig =
  132.                 new ParameterConfiguration(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT, hMu);
  133.         jacobianizer = new Jacobianizer(this, mu, Collections.singletonList(muConfig), hPosition);
  134.     }

  136.     /** {@inheritDoc} */
  137.     public TideSystem getTideSystem() {
  138.         return provider.getTideSystem();
  139.     }

  141.     /** {@inheritDoc} */
  142.     public void addContribution(final SpacecraftState s, final TimeDerivativesEquations adder)
  143.         throws OrekitException {

  145.         // get the position in body frame
  146.         final AbsoluteDate date = s.getDate();
  147.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
  148.         final Transform fromBodyFrame = bodyFrame.getTransformTo(s.getFrame(), date);
  149.         final Transform toBodyFrame   = fromBodyFrame.getInverse();
  150.         final Vector3D relative = toBodyFrame.transformPosition(s.getPVCoordinates().getPosition());

  152.         final double x = relative.getX();
  153.         final double y = relative.getY();
  154.         final double z = relative.getZ();

  156.         final double x2 = x * x;
  157.         final double y2 = y * y;
  158.         final double z2 = z * z;
  159.         final double r2 = x2 + y2 + z2;
  160.         final double r = FastMath.sqrt(r2);
  161.         final double equatorialRadius = provider.getAe();
  162.         if (r <= equatorialRadius) {
  163.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r);
  164.         }

  166.         // define some intermediate variables
  167.         final double onR2 = 1 / r2;
  168.         final double onR3 = onR2 / r;
  169.         final double rEqOnR2  = equatorialRadius / r2;
  170.         final double rEqOnR4  = rEqOnR2 / r2;
  171.         final double rEq2OnR2 = equatorialRadius * rEqOnR2;

  173.         double cmx   = -x * rEqOnR2;
  174.         double cmy   = -y * rEqOnR2;
  175.         double cmz   = -z * rEqOnR2;

  177.         final double dx   = -2 * cmx;
  178.         final double dy   = -2 * cmy;
  179.         final double dz   = -2 * cmz;

  181.         // intermediate variables gradients
  182.         // since dcy/dx = dcx/dy, dcz/dx = dcx/dz and dcz/dy = dcy/dz,
  183.         // we reuse the existing variables

  185.         double dcmxdx = (x2 - y2 - z2) * rEqOnR4;
  186.         double dcmxdy =  dx * y * onR2;
  187.         double dcmxdz =  dx * z * onR2;
  188.         double dcmydy = (y2 - x2 - z2) * rEqOnR4;
  189.         double dcmydz =  dy * z * onR2;
  190.         double dcmzdz = (z2 - x2 - y2) * rEqOnR4;

  192.         final double ddxdx = -2 * dcmxdx;
  193.         final double ddxdy = -2 * dcmxdy;
  194.         final double ddxdz = -2 * dcmxdz;
  195.         final double ddydy = -2 * dcmydy;
  196.         final double ddydz = -2 * dcmydz;
  197.         final double ddzdz = -2 * dcmzdz;

  199.         final double donr2dx = -dx * rEqOnR2;
  200.         final double donr2dy = -dy * rEqOnR2;
  201.         final double donr2dz = -dz * rEqOnR2;

  203.         // potential coefficients (4 per matrix)
  204.         double vrn  = 0.0;
  205.         double vin  = 0.0;
  206.         double vrd  = 1.0 / (equatorialRadius * r);
  207.         double vid  = 0.0;
  208.         double vrn1 = 0.0;
  209.         double vin1 = 0.0;
  210.         double vrn2 = 0.0;
  211.         double vin2 = 0.0;

  213.         // gradient coefficients (4 per matrix)
  214.         double gradXVrn  = 0.0;
  215.         double gradXVin  = 0.0;
  216.         double gradXVrd  = -x * onR3 / equatorialRadius;
  217.         double gradXVid  = 0.0;
  218.         double gradXVrn1 = 0.0;
  219.         double gradXVin1 = 0.0;
  220.         double gradXVrn2 = 0.0;
  221.         double gradXVin2 = 0.0;

  223.         double gradYVrn  = 0.0;
  224.         double gradYVin  = 0.0;
  225.         double gradYVrd  = -y * onR3 / equatorialRadius;
  226.         double gradYVid  = 0.0;
  227.         double gradYVrn1 = 0.0;
  228.         double gradYVin1 = 0.0;
  229.         double gradYVrn2 = 0.0;
  230.         double gradYVin2 = 0.0;

  232.         double gradZVrn  = 0.0;
  233.         double gradZVin  = 0.0;
  234.         double gradZVrd  = -z * onR3 / equatorialRadius;
  235.         double gradZVid  = 0.0;
  236.         double gradZVrn1 = 0.0;
  237.         double gradZVin1 = 0.0;
  238.         double gradZVrn2 = 0.0;
  239.         double gradZVin2 = 0.0;

  241.         // acceleration coefficients
  242.         double vdX = 0.0;
  243.         double vdY = 0.0;
  244.         double vdZ = 0.0;

  246.         // start calculating
  247.         for (int m = 0; m <= provider.getMaxOrder(); m++) {

  249.             double cx = cmx;
  250.             double cy = cmy;
  251.             double cz = cmz;

  253.             double dcxdx = dcmxdx;
  254.             double dcxdy = dcmxdy;
  255.             double dcxdz = dcmxdz;
  256.             double dcydy = dcmydy;
  257.             double dcydz = dcmydz;
  258.             double dczdz = dcmzdz;

  260.             for (int n = m; n <= provider.getMaxDegree(); n++) {

  262.                 if (n == m) {
  263.                     // calculate the first element of the next column

  265.                     vrn      = equatorialRadius * vrd;
  266.                     vin      = equatorialRadius * vid;

  268.                     gradXVrn = equatorialRadius * gradXVrd;
  269.                     gradXVin = equatorialRadius * gradXVid;
  270.                     gradYVrn = equatorialRadius * gradYVrd;
  271.                     gradYVin = equatorialRadius * gradYVid;
  272.                     gradZVrn = equatorialRadius * gradZVrd;
  273.                     gradZVin = equatorialRadius * gradZVid;

  275.                     final double tmpGradXVrd = (cx + dx) * gradXVrd - (cy + dy) * gradXVid + (dcxdx + ddxdx) * vrd - (dcxdy + ddxdy) * vid;
  276.                     gradXVid = (cy + dy) * gradXVrd + (cx + dx) * gradXVid + (dcxdy + ddxdy) * vrd + (dcxdx + ddxdx) * vid;
  277.                     gradXVrd = tmpGradXVrd;

  279.                     final double tmpGradYVrd = (cx + dx) * gradYVrd - (cy + dy) * gradYVid + (dcxdy + ddxdy) * vrd - (dcydy + ddydy) * vid;
  280.                     gradYVid = (cy + dy) * gradYVrd + (cx + dx) * gradYVid + (dcydy + ddydy) * vrd + (dcxdy + ddxdy) * vid;
  281.                     gradYVrd = tmpGradYVrd;

  283.                     final double tmpGradZVrd = (cx + dx) * gradZVrd - (cy + dy) * gradZVid + (dcxdz + ddxdz) * vrd - (dcydz + ddydz) * vid;
  284.                     gradZVid = (cy + dy) * gradZVrd + (cx + dx) * gradZVid + (dcydz + ddydz) * vrd + (dcxdz + ddxdz) * vid;
  285.                     gradZVrd = tmpGradZVrd;

  287.                     final double tmpVrd = (cx + dx) * vrd - (cy + dy) * vid;
  288.                     vid = (cy + dy) * vrd + (cx + dx) * vid;
  289.                     vrd = tmpVrd;

  291.                 } else if (n == m + 1) {
  292.                     // calculate the second element of the column
  293.                     vrn = cz * vrn1;
  294.                     vin = cz * vin1;

  296.                     gradXVrn = cz * gradXVrn1 + dcxdz * vrn1;
  297.                     gradXVin = cz * gradXVin1 + dcxdz * vin1;

  299.                     gradYVrn = cz * gradYVrn1 + dcydz * vrn1;
  300.                     gradYVin = cz * gradYVin1 + dcydz * vin1;

  302.                     gradZVrn = cz * gradZVrn1 + dczdz * vrn1;
  303.                     gradZVin = cz * gradZVin1 + dczdz * vin1;

  305.                 } else {
  306.                     // calculate the other elements of the column
  307.                     final double inv   = 1.0 / (n - m);
  308.                     final double coeff = n + m - 1.0;

  310.                     vrn = (cz * vrn1 - coeff * rEq2OnR2 * vrn2) * inv;
  311.                     vin = (cz * vin1 - coeff * rEq2OnR2 * vin2) * inv;

  313.                     gradXVrn = (cz * gradXVrn1 - coeff * rEq2OnR2 * gradXVrn2 + dcxdz * vrn1 - coeff * donr2dx * vrn2) * inv;
  314.                     gradXVin = (cz * gradXVin1 - coeff * rEq2OnR2 * gradXVin2 + dcxdz * vin1 - coeff * donr2dx * vin2) * inv;
  315.                     gradYVrn = (cz * gradYVrn1 - coeff * rEq2OnR2 * gradYVrn2 + dcydz * vrn1 - coeff * donr2dy * vrn2) * inv;
  316.                     gradYVin = (cz * gradYVin1 - coeff * rEq2OnR2 * gradYVin2 + dcydz * vin1 - coeff * donr2dy * vin2) * inv;
  317.                     gradZVrn = (cz * gradZVrn1 - coeff * rEq2OnR2 * gradZVrn2 + dczdz * vrn1 - coeff * donr2dz * vrn2) * inv;
  318.                     gradZVin = (cz * gradZVin1 - coeff * rEq2OnR2 * gradZVin2 + dczdz * vin1 - coeff * donr2dz * vin2) * inv;
  319.                 }

  321.                 // increment variables
  322.                 cx += dx;
  323.                 cy += dy;
  324.                 cz += dz;

  326.                 dcxdx += ddxdx;
  327.                 dcxdy += ddxdy;
  328.                 dcxdz += ddxdz;
  329.                 dcydy += ddydy;
  330.                 dcydz += ddydz;
  331.                 dczdz += ddzdz;

  333.                 vrn2 = vrn1;
  334.                 vin2 = vin1;
  335.                 gradXVrn2 = gradXVrn1;
  336.                 gradXVin2 = gradXVin1;
  337.                 gradYVrn2 = gradYVrn1;
  338.                 gradYVin2 = gradYVin1;
  339.                 gradZVrn2 = gradZVrn1;
  340.                 gradZVin2 = gradZVin1;

  342.                 vrn1 = vrn;
  343.                 vin1 = vin;
  344.                 gradXVrn1 = gradXVrn;
  345.                 gradXVin1 = gradXVin;
  346.                 gradYVrn1 = gradYVrn;
  347.                 gradYVin1 = gradYVin;
  348.                 gradZVrn1 = gradZVrn;
  349.                 gradZVin1 = gradZVin;

  351.                 // compute the acceleration due to the Cnm and Snm coefficients
  352.                 // ignoring the central attraction
  353.                 if (n > 0) {
  354.                     final double cnm = harmonics.getUnnormalizedCnm(n, m);
  355.                     final double snm = harmonics.getUnnormalizedSnm(n, m);
  356.                     vdX += cnm * gradXVrn + snm * gradXVin;
  357.                     vdY += cnm * gradYVrn + snm * gradYVin;
  358.                     vdZ += cnm * gradZVrn + snm * gradZVin;
  359.                 }

  361.             }

  363.             // increment variables
  364.             cmx += dx;
  365.             cmy += dy;
  366.             cmz += dz;

  368.             dcmxdx += ddxdx;
  369.             dcmxdy += ddxdy;
  370.             dcmxdz += ddxdz;
  371.             dcmydy += ddydy;
  372.             dcmydz += ddydz;
  373.             dcmzdz += ddzdz;

  375.         }

  377.         // compute acceleration in inertial frame
  378.         final Vector3D acceleration =
  379.             fromBodyFrame.transformVector(new Vector3D(mu * vdX, mu * vdY, mu * vdZ));
  380.         adder.addXYZAcceleration(acceleration.getX(), acceleration.getY(), acceleration.getZ());

  382.     }

  384.     /** {@inheritDoc} */
  385.     public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final  Frame frame,
  386.                                                                       final FieldVector3D<DerivativeStructure> position,
  387.                                                                       final FieldVector3D<DerivativeStructure> velocity,
  388.                                                                       final FieldRotation<DerivativeStructure> rotation,
  389.                                                                       final DerivativeStructure mass)
  390.         throws OrekitException {
  391.         if (jacobianizer == null) {
  392.             throw new OrekitException(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES);
  393.         }
  394.         return jacobianizer.accelerationDerivatives(date, frame, position, velocity, rotation, mass);
  395.     }

  397.     /** {@inheritDoc} */
  398.     public FieldVector3D<DerivativeStructure> accelerationDerivatives(final SpacecraftState s, final String paramName)
  399.         throws OrekitException {
  400.         if (jacobianizer == null) {
  401.             throw new OrekitException(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES);
  402.         }
  403.         return jacobianizer.accelerationDerivatives(s, paramName);
  404.     }

  406.     /** {@inheritDoc} */
  407.     public EventDetector[] getEventsDetectors() {
  408.         return new EventDetector[0];
  409.     }

  411.     /** {@inheritDoc} */
  412.     public double getParameter(final String name)
  413.         throws IllegalArgumentException {
  414.         complainIfNotSupported(name);
  415.         return mu;
  416.     }

  418.     /** {@inheritDoc} */
  419.     public void setParameter(final String name, final double value)
  420.         throws IllegalArgumentException {
  421.         complainIfNotSupported(name);
  422.         mu = value;
  423.     }

  425. }
Created with JaCoCo 0.6.3.201306030806