CR3BPForceModel.java
/* Copyright 2002-2024 CS GROUP
* Licensed to CS GROUP (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.propagation.numerical.cr3bp;
import java.util.Collections;
import java.util.List;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.analysis.differentiation.DSFactory;
import org.hipparchus.analysis.differentiation.DerivativeStructure;
import org.hipparchus.analysis.differentiation.FDSFactory;
import org.hipparchus.analysis.differentiation.FieldDerivativeStructure;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.CR3BPSystem;
import org.orekit.forces.ForceModel;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.utils.ParameterDriver;
/** Class calculating the acceleration induced by CR3BP model.
* @see "Dynamical systems, the three-body problem, and space mission design, Koon, Lo, Marsden, Ross"
* @author Vincent Mouraux
* @since 10.2
*/
public class CR3BPForceModel implements ForceModel {
/** Suffix for parameter name for Mass Ratio enabling Jacobian processing. */
public static final String MASS_RATIO_SUFFIX = "CR3BP System Mass Ratio";
/**
* Central attraction scaling factor.
* <p>
* We use a power of 2 to avoid numeric noise introduction in the
* multiplications/divisions sequences.
* </p>
*/
private static final double MU_SCALE = FastMath.scalb(1.0, 32);
/** Driver for gravitational parameter. */
private final ParameterDriver muParameterDriver;
/** Simple constructor.
* @param cr3bp Name of the CR3BP System
*/
public CR3BPForceModel(final CR3BPSystem cr3bp) {
muParameterDriver = new ParameterDriver(cr3bp.getName() + MASS_RATIO_SUFFIX,
cr3bp.getMassRatio(), MU_SCALE, 0.0,
Double.POSITIVE_INFINITY);
}
/** {@inheritDoc} */
public Vector3D acceleration(final SpacecraftState s,
final double[] parameters) {
// Spacecraft Velocity
final double vx = s.getPVCoordinates().getVelocity().getX();
final double vy = s.getPVCoordinates().getVelocity().getY();
// Spacecraft Potential
final DerivativeStructure potential = getPotential(s);
// Potential derivatives
final double[] dU = potential.getAllDerivatives();
// first order derivatives index
final int idX = potential.getFactory().getCompiler().getPartialDerivativeIndex(1, 0, 0);
final int idY = potential.getFactory().getCompiler().getPartialDerivativeIndex(0, 1, 0);
final int idZ = potential.getFactory().getCompiler().getPartialDerivativeIndex(0, 0, 1);
// Acceleration calculation according to CR3BP Analytical Model
final double accx = dU[idX] + 2.0 * vy;
final double accy = dU[idY] - 2.0 * vx;
final double accz = dU[idZ];
// compute absolute acceleration
return new Vector3D(accx, accy, accz);
}
/** {@inheritDoc} */
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
final T[] parameters) {
// Spacecraft Velocity
final T vx = s.getPVCoordinates().getVelocity().getX();
final T vy = s.getPVCoordinates().getVelocity().getY();
// Spacecraft Potential
final FieldDerivativeStructure<T> fieldPotential = getPotential(s);
// Potential derivatives
final T[] dU = fieldPotential.getAllDerivatives();
// first order derivatives index
final int idX = fieldPotential.getFactory().getCompiler().getPartialDerivativeIndex(1, 0, 0);
final int idY = fieldPotential.getFactory().getCompiler().getPartialDerivativeIndex(0, 1, 0);
final int idZ = fieldPotential.getFactory().getCompiler().getPartialDerivativeIndex(0, 0, 1);
// Acceleration calculation according to CR3BP Analytical Model
final T accx = dU[idX].add(vy.multiply(2.0));
final T accy = dU[idY].subtract(vx.multiply(2.0));
final T accz = dU[idZ];
// compute absolute acceleration
return new FieldVector3D<>(accx, accy, accz);
}
/**
* Calculate spacecraft potential.
* @param s SpacecraftState
* @return Spacecraft Potential
*/
public DerivativeStructure getPotential(final SpacecraftState s) {
// Spacecraft Position
final double x = s.getPosition().getX();
final double y = s.getPosition().getY();
final double z = s.getPosition().getZ();
final DSFactory factoryP = new DSFactory(3, 2);
final DerivativeStructure fpx = factoryP.variable(0, x);
final DerivativeStructure fpy = factoryP.variable(1, y);
final DerivativeStructure fpz = factoryP.variable(2, z);
final DerivativeStructure zero = fpx.getField().getZero();
// Get CR3BP System mass ratio
// By construction, mudriver has 1 value for the all time period that is why
// the getValue can be called with any date argument or null argument
final DerivativeStructure mu = zero.newInstance(muParameterDriver.getValue(s.getDate()));
// Normalized distances between primaries and barycenter in CR3BP
final DerivativeStructure d1 = mu;
final DerivativeStructure d2 = mu.negate().add(1.0);
// Norm of the Spacecraft position relative to the primary body
final DerivativeStructure r1 =
FastMath.sqrt((fpx.add(d1)).multiply(fpx.add(d1)).add(fpy.square())
.add(fpz.square()));
// Norm of the Spacecraft position relative to the secondary body
final DerivativeStructure r2 =
FastMath.sqrt((fpx.subtract(d2)).multiply(fpx.subtract(d2))
.add(fpy.square()).add(fpz.square()));
// Potential of the Spacecraft
return (mu.negate().add(1.0).divide(r1)).add(mu.divide(r2))
.add(fpx.square().add(fpy.square()).multiply(0.5)).add(d1.multiply(d2).multiply(0.5));
}
/**
* Calculate spacecraft potential.
* @param <T> Field element
* @param s SpacecraftState
* @return Spacecraft Potential
*/
public <T extends CalculusFieldElement<T>> FieldDerivativeStructure<T> getPotential(final FieldSpacecraftState<T> s) {
// Spacecraft Position
final T x = s.getPosition().getX();
final T y = s.getPosition().getY();
final T z = s.getPosition().getZ();
final FDSFactory<T> factoryP = new FDSFactory<>(s.getDate().getField(), 3, 2);
final FieldDerivativeStructure<T> fpx = factoryP.variable(0, x);
final FieldDerivativeStructure<T> fpy = factoryP.variable(1, y);
final FieldDerivativeStructure<T> fpz = factoryP.variable(2, z);
final FieldDerivativeStructure<T> zero = fpx.getField().getZero();
// Get CR3BP System mass ratio
// By construction, mudriver has 1 value for the all time period that is why
// the getValue can be called with any date argument or null argument
final FieldDerivativeStructure<T> mu = zero.newInstance(muParameterDriver.getValue(s.getDate().toAbsoluteDate()));
// Normalized distances between primaries and barycenter in CR3BP
final FieldDerivativeStructure<T> d1 = mu;
final FieldDerivativeStructure<T> d2 = mu.negate().add(1.0);
// Norm of the Spacecraft position relative to the primary body
final FieldDerivativeStructure<T> r1 =
FastMath.sqrt((fpx.add(d1)).multiply(fpx.add(d1)).add(fpy.square())
.add(fpz.square()));
// Norm of the Spacecraft position relative to the secondary body
final FieldDerivativeStructure<T> r2 =
FastMath.sqrt((fpx.subtract(d2)).multiply(fpx.subtract(d2))
.add(fpy.square()).add(fpz.square()));
// Potential of the Spacecraft
return (mu.negate().add(1.0).divide(r1)).add(mu.divide(r2))
.add(fpx.square().add(fpy.square()).multiply(0.5)).add(d1.multiply(d2).multiply(0.5));
}
/** {@inheritDoc} */
public List<ParameterDriver> getParametersDrivers() {
return Collections.singletonList(muParameterDriver);
}
/** {@inheritDoc} */
public boolean dependsOnPositionOnly() {
return true;
}
}