CunninghamAttractionModel.java
/* Copyright 2002-2013 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.forces.gravity;
import java.util.Collections;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.geometry.euclidean.threed.FieldRotation;
import org.apache.commons.math3.geometry.euclidean.threed.FieldVector3D;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.ode.AbstractParameterizable;
import org.apache.commons.math3.util.FastMath;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.forces.ForceModel;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.TideSystem;
import org.orekit.forces.gravity.potential.TideSystemProvider;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
import org.orekit.frames.Frame;
import org.orekit.frames.Transform;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.numerical.Jacobianizer;
import org.orekit.propagation.numerical.ParameterConfiguration;
import org.orekit.propagation.numerical.TimeDerivativesEquations;
import org.orekit.time.AbsoluteDate;
/** This class represents the gravitational field of a celestial body.
* <p>The algorithm implemented in this class has been designed by
* Leland E. Cunningham (Lockheed Missiles and Space Company, Sunnyvale
* and Astronomy Department University of California, Berkeley) in
* his 1969 paper: <em>On the computation of the spherical harmonic
* terms needed during the numerical integration of the orbital motion
* of an artificial satellite</em> (Celestial Mechanics 2, 1970).</p>
* <p>
* Note that this class can often not be used for high degrees (say
* above 90) as most modern gravity fields are provided as normalized
* coefficients and the un-normalization process to convert these
* coefficients underflows at degree and order 89. This class also
* does not provide analytical partial derivatives (it uses finite differences
* to compute them) and is much slower than {@link HolmesFeatherstoneAttractionModel}
* (even when no derivatives are computed). For all these reasons,
* it is recommended to use the {@link HolmesFeatherstoneAttractionModel
* Holmes-Featherstone model} rather than this class.
* </p>
* <p>
* As this class uses finite differences to compute derivatives, the steps for
* finite differences <strong>must</strong> be initialized by calling {@link
* #setSteps(double, double)} prior to use derivatives, otherwise an exception
* will be thrown by {@link #accelerationDerivatives(AbsoluteDate, Frame, FieldVector3D,
* FieldVector3D, FieldRotation, DerivativeStructure)} and by {@link
* #accelerationDerivatives(SpacecraftState, String)}.
* </p>
*
* @see HolmesFeatherstoneAttractionModel
* @author Fabien Maussion
* @author Luc Maisonobe
* @author Véronique Pommier-Maurussane
*/
public class CunninghamAttractionModel
extends AbstractParameterizable implements ForceModel, TideSystemProvider {
/** Provider for the spherical harmonics. */
private final UnnormalizedSphericalHarmonicsProvider provider;
/** Central attraction coefficient. */
private double mu;
/** Rotating body. */
private final Frame bodyFrame;
/** Helper class computing acceleration derivatives. */
private Jacobianizer jacobianizer;
/** Creates a new instance.
* @param centralBodyFrame rotating body frame
* @param equatorialRadius reference equatorial radius of the potential
* @param mu central body attraction coefficient (m<sup>3</sup>/s<sup>2</sup>)
* @param C un-normalized coefficients array (cosine part)
* @param S un-normalized coefficients array (sine part)
* @deprecated since 6.0, replaced by {@link #CunninghamAttractionModel(Frame, UnnormalizedSphericalHarmonicsProvider)}
*/
public CunninghamAttractionModel(final Frame centralBodyFrame,
final double equatorialRadius, final double mu,
final double[][] C, final double[][] S) {
this(centralBodyFrame, GravityFieldFactory.getUnnormalizedProvider(equatorialRadius, mu,
TideSystem.UNKNOWN, C, S));
}
/** Creates a new instance.
* @param centralBodyFrame rotating body frame
* @param provider provider for spherical harmonics
* @since 6.0
*/
public CunninghamAttractionModel(final Frame centralBodyFrame,
final UnnormalizedSphericalHarmonicsProvider provider) {
super(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT);
this.provider = provider;
this.mu = provider.getMu();
this.bodyFrame = centralBodyFrame;
this.jacobianizer = null;
}
/** Set the step for finite differences with respect to spacecraft position.
* @param hPosition step used for finite difference computation
* with respect to spacecraft position (m)
* @param hMu step used for finite difference computation
* with respect to central attraction coefficient (m<sup>3</sup>/s<sup>2</sup>)
*/
public void setSteps(final double hPosition, final double hMu) {
final ParameterConfiguration muConfig =
new ParameterConfiguration(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT, hMu);
jacobianizer = new Jacobianizer(this, mu, Collections.singletonList(muConfig), hPosition);
}
/** {@inheritDoc} */
public TideSystem getTideSystem() {
return provider.getTideSystem();
}
/** {@inheritDoc} */
public void addContribution(final SpacecraftState s, final TimeDerivativesEquations adder)
throws OrekitException {
// get the position in body frame
final AbsoluteDate date = s.getDate();
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
final Transform fromBodyFrame = bodyFrame.getTransformTo(s.getFrame(), date);
final Transform toBodyFrame = fromBodyFrame.getInverse();
final Vector3D relative = toBodyFrame.transformPosition(s.getPVCoordinates().getPosition());
final double x = relative.getX();
final double y = relative.getY();
final double z = relative.getZ();
final double x2 = x * x;
final double y2 = y * y;
final double z2 = z * z;
final double r2 = x2 + y2 + z2;
final double r = FastMath.sqrt(r2);
final double equatorialRadius = provider.getAe();
if (r <= equatorialRadius) {
throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r);
}
// define some intermediate variables
final double onR2 = 1 / r2;
final double onR3 = onR2 / r;
final double rEqOnR2 = equatorialRadius / r2;
final double rEqOnR4 = rEqOnR2 / r2;
final double rEq2OnR2 = equatorialRadius * rEqOnR2;
double cmx = -x * rEqOnR2;
double cmy = -y * rEqOnR2;
double cmz = -z * rEqOnR2;
final double dx = -2 * cmx;
final double dy = -2 * cmy;
final double dz = -2 * cmz;
// intermediate variables gradients
// since dcy/dx = dcx/dy, dcz/dx = dcx/dz and dcz/dy = dcy/dz,
// we reuse the existing variables
double dcmxdx = (x2 - y2 - z2) * rEqOnR4;
double dcmxdy = dx * y * onR2;
double dcmxdz = dx * z * onR2;
double dcmydy = (y2 - x2 - z2) * rEqOnR4;
double dcmydz = dy * z * onR2;
double dcmzdz = (z2 - x2 - y2) * rEqOnR4;
final double ddxdx = -2 * dcmxdx;
final double ddxdy = -2 * dcmxdy;
final double ddxdz = -2 * dcmxdz;
final double ddydy = -2 * dcmydy;
final double ddydz = -2 * dcmydz;
final double ddzdz = -2 * dcmzdz;
final double donr2dx = -dx * rEqOnR2;
final double donr2dy = -dy * rEqOnR2;
final double donr2dz = -dz * rEqOnR2;
// potential coefficients (4 per matrix)
double vrn = 0.0;
double vin = 0.0;
double vrd = 1.0 / (equatorialRadius * r);
double vid = 0.0;
double vrn1 = 0.0;
double vin1 = 0.0;
double vrn2 = 0.0;
double vin2 = 0.0;
// gradient coefficients (4 per matrix)
double gradXVrn = 0.0;
double gradXVin = 0.0;
double gradXVrd = -x * onR3 / equatorialRadius;
double gradXVid = 0.0;
double gradXVrn1 = 0.0;
double gradXVin1 = 0.0;
double gradXVrn2 = 0.0;
double gradXVin2 = 0.0;
double gradYVrn = 0.0;
double gradYVin = 0.0;
double gradYVrd = -y * onR3 / equatorialRadius;
double gradYVid = 0.0;
double gradYVrn1 = 0.0;
double gradYVin1 = 0.0;
double gradYVrn2 = 0.0;
double gradYVin2 = 0.0;
double gradZVrn = 0.0;
double gradZVin = 0.0;
double gradZVrd = -z * onR3 / equatorialRadius;
double gradZVid = 0.0;
double gradZVrn1 = 0.0;
double gradZVin1 = 0.0;
double gradZVrn2 = 0.0;
double gradZVin2 = 0.0;
// acceleration coefficients
double vdX = 0.0;
double vdY = 0.0;
double vdZ = 0.0;
// start calculating
for (int m = 0; m <= provider.getMaxOrder(); m++) {
double cx = cmx;
double cy = cmy;
double cz = cmz;
double dcxdx = dcmxdx;
double dcxdy = dcmxdy;
double dcxdz = dcmxdz;
double dcydy = dcmydy;
double dcydz = dcmydz;
double dczdz = dcmzdz;
for (int n = m; n <= provider.getMaxDegree(); n++) {
if (n == m) {
// calculate the first element of the next column
vrn = equatorialRadius * vrd;
vin = equatorialRadius * vid;
gradXVrn = equatorialRadius * gradXVrd;
gradXVin = equatorialRadius * gradXVid;
gradYVrn = equatorialRadius * gradYVrd;
gradYVin = equatorialRadius * gradYVid;
gradZVrn = equatorialRadius * gradZVrd;
gradZVin = equatorialRadius * gradZVid;
final double tmpGradXVrd = (cx + dx) * gradXVrd - (cy + dy) * gradXVid + (dcxdx + ddxdx) * vrd - (dcxdy + ddxdy) * vid;
gradXVid = (cy + dy) * gradXVrd + (cx + dx) * gradXVid + (dcxdy + ddxdy) * vrd + (dcxdx + ddxdx) * vid;
gradXVrd = tmpGradXVrd;
final double tmpGradYVrd = (cx + dx) * gradYVrd - (cy + dy) * gradYVid + (dcxdy + ddxdy) * vrd - (dcydy + ddydy) * vid;
gradYVid = (cy + dy) * gradYVrd + (cx + dx) * gradYVid + (dcydy + ddydy) * vrd + (dcxdy + ddxdy) * vid;
gradYVrd = tmpGradYVrd;
final double tmpGradZVrd = (cx + dx) * gradZVrd - (cy + dy) * gradZVid + (dcxdz + ddxdz) * vrd - (dcydz + ddydz) * vid;
gradZVid = (cy + dy) * gradZVrd + (cx + dx) * gradZVid + (dcydz + ddydz) * vrd + (dcxdz + ddxdz) * vid;
gradZVrd = tmpGradZVrd;
final double tmpVrd = (cx + dx) * vrd - (cy + dy) * vid;
vid = (cy + dy) * vrd + (cx + dx) * vid;
vrd = tmpVrd;
} else if (n == m + 1) {
// calculate the second element of the column
vrn = cz * vrn1;
vin = cz * vin1;
gradXVrn = cz * gradXVrn1 + dcxdz * vrn1;
gradXVin = cz * gradXVin1 + dcxdz * vin1;
gradYVrn = cz * gradYVrn1 + dcydz * vrn1;
gradYVin = cz * gradYVin1 + dcydz * vin1;
gradZVrn = cz * gradZVrn1 + dczdz * vrn1;
gradZVin = cz * gradZVin1 + dczdz * vin1;
} else {
// calculate the other elements of the column
final double inv = 1.0 / (n - m);
final double coeff = n + m - 1.0;
vrn = (cz * vrn1 - coeff * rEq2OnR2 * vrn2) * inv;
vin = (cz * vin1 - coeff * rEq2OnR2 * vin2) * inv;
gradXVrn = (cz * gradXVrn1 - coeff * rEq2OnR2 * gradXVrn2 + dcxdz * vrn1 - coeff * donr2dx * vrn2) * inv;
gradXVin = (cz * gradXVin1 - coeff * rEq2OnR2 * gradXVin2 + dcxdz * vin1 - coeff * donr2dx * vin2) * inv;
gradYVrn = (cz * gradYVrn1 - coeff * rEq2OnR2 * gradYVrn2 + dcydz * vrn1 - coeff * donr2dy * vrn2) * inv;
gradYVin = (cz * gradYVin1 - coeff * rEq2OnR2 * gradYVin2 + dcydz * vin1 - coeff * donr2dy * vin2) * inv;
gradZVrn = (cz * gradZVrn1 - coeff * rEq2OnR2 * gradZVrn2 + dczdz * vrn1 - coeff * donr2dz * vrn2) * inv;
gradZVin = (cz * gradZVin1 - coeff * rEq2OnR2 * gradZVin2 + dczdz * vin1 - coeff * donr2dz * vin2) * inv;
}
// increment variables
cx += dx;
cy += dy;
cz += dz;
dcxdx += ddxdx;
dcxdy += ddxdy;
dcxdz += ddxdz;
dcydy += ddydy;
dcydz += ddydz;
dczdz += ddzdz;
vrn2 = vrn1;
vin2 = vin1;
gradXVrn2 = gradXVrn1;
gradXVin2 = gradXVin1;
gradYVrn2 = gradYVrn1;
gradYVin2 = gradYVin1;
gradZVrn2 = gradZVrn1;
gradZVin2 = gradZVin1;
vrn1 = vrn;
vin1 = vin;
gradXVrn1 = gradXVrn;
gradXVin1 = gradXVin;
gradYVrn1 = gradYVrn;
gradYVin1 = gradYVin;
gradZVrn1 = gradZVrn;
gradZVin1 = gradZVin;
// compute the acceleration due to the Cnm and Snm coefficients
// ignoring the central attraction
if (n > 0) {
final double cnm = harmonics.getUnnormalizedCnm(n, m);
final double snm = harmonics.getUnnormalizedSnm(n, m);
vdX += cnm * gradXVrn + snm * gradXVin;
vdY += cnm * gradYVrn + snm * gradYVin;
vdZ += cnm * gradZVrn + snm * gradZVin;
}
}
// increment variables
cmx += dx;
cmy += dy;
cmz += dz;
dcmxdx += ddxdx;
dcmxdy += ddxdy;
dcmxdz += ddxdz;
dcmydy += ddydy;
dcmydz += ddydz;
dcmzdz += ddzdz;
}
// compute acceleration in inertial frame
final Vector3D acceleration =
fromBodyFrame.transformVector(new Vector3D(mu * vdX, mu * vdY, mu * vdZ));
adder.addXYZAcceleration(acceleration.getX(), acceleration.getY(), acceleration.getZ());
}
/** {@inheritDoc} */
public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame,
final FieldVector3D<DerivativeStructure> position,
final FieldVector3D<DerivativeStructure> velocity,
final FieldRotation<DerivativeStructure> rotation,
final DerivativeStructure mass)
throws OrekitException {
if (jacobianizer == null) {
throw new OrekitException(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES);
}
return jacobianizer.accelerationDerivatives(date, frame, position, velocity, rotation, mass);
}
/** {@inheritDoc} */
public FieldVector3D<DerivativeStructure> accelerationDerivatives(final SpacecraftState s, final String paramName)
throws OrekitException {
if (jacobianizer == null) {
throw new OrekitException(OrekitMessages.STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES);
}
return jacobianizer.accelerationDerivatives(s, paramName);
}
/** {@inheritDoc} */
public EventDetector[] getEventsDetectors() {
return new EventDetector[0];
}
/** {@inheritDoc} */
public double getParameter(final String name)
throws IllegalArgumentException {
complainIfNotSupported(name);
return mu;
}
/** {@inheritDoc} */
public void setParameter(final String name, final double value)
throws IllegalArgumentException {
complainIfNotSupported(name);
mu = value;
}
}