CR3BPMultipleShooter.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.List;
import java.util.Map;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.utils.AbsolutePVCoordinates;
import org.orekit.utils.AbstractMultipleShooting;
/**
* Multiple shooting method applicable for orbits, either propagation in CR3BP, or in an ephemeris model.
* @see "TRAJECTORY DESIGN AND ORBIT MAINTENANCE STRATEGIES IN MULTI-BODY DYNAMICAL REGIMES by Thomas A. Pavlak, Purdue University"
* @author William Desprats
* @author Alberto Fossà
*/
public class CR3BPMultipleShooter extends AbstractMultipleShooting {
/** Name of the derivatives. */
private static final String STM = "stmEquations";
/** Derivatives linked to the Propagators.
* @since 11.1
*/
private final List<STMEquations> stmEquations;
/** True if orbit is closed. */
private boolean isClosedOrbit;
/** Simple Constructor.
* <p> Standard constructor for multiple shooting which can be used with the CR3BP model. </p>
* @param initialGuessList initial patch points to be corrected
* @param propagatorList list of propagators associated to each patch point
* @param stmEquations list of additional derivatives providers linked to propagatorList
* @param tolerance convergence tolerance on the constraint vector
* @param maxIter maximum number of iterations
*/
public CR3BPMultipleShooter(final List<SpacecraftState> initialGuessList,
final List<NumericalPropagator> propagatorList,
final List<STMEquations> stmEquations,
final double tolerance, final int maxIter) {
super(initialGuessList, propagatorList, tolerance, maxIter, true, STM);
this.stmEquations = stmEquations;
}
/** {@inheritDoc} */
protected SpacecraftState getAugmentedInitialState(final int i) {
return stmEquations.get(i).setInitialPhi(getPatchPoint(i));
}
/** {@inheritDoc} */
@Override
public void setEpochFreedom(final int patchIndex, final boolean isFree) {
throw new OrekitException(OrekitMessages.FUNCTION_NOT_IMPLEMENTED);
}
/** {@inheritDoc} */
@Override
public void setScaleLength(final double scaleLength) {
throw new OrekitException(OrekitMessages.FUNCTION_NOT_IMPLEMENTED);
}
/** {@inheritDoc} */
@Override
public void setScaleTime(final double scaleTime) {
throw new OrekitException(OrekitMessages.FUNCTION_NOT_IMPLEMENTED);
}
/** Set the constraint of a closed orbit or not.
* @param isClosed true if orbit should be closed
*/
public void setClosedOrbitConstraint(final boolean isClosed) {
this.isClosedOrbit = isClosed;
}
/** {@inheritDoc} */
@Override
protected double[][] computeAdditionalJacobianMatrix(final List<SpacecraftState> propagatedSP) {
final Map<Integer, Double> mapConstraints = getConstraintsMap();
final boolean[] freeCompsMap = getFreeCompsMap();
// Number of additional constraints
final int nRows = mapConstraints.size() + (isClosedOrbit ? 6 : 0);
final int nCols = getNumberOfFreeComponents();
final double[][] M = new double[nRows][nCols];
int j = 0;
if (isClosedOrbit) {
// The Jacobian matrix has the following form:
//
// [-1 0 0 ... 1 0 ]
// [ 0 -1 0 0 ... 1 0 ]
// F = [ 0 -1 0 0 ... 1 0 ]
// [ 0 -1 0 0 ... 1 0 ]
// [ 0 -1 0 0 ... 1 0 ]
// [ 0 0 -1 0 ... 0 1 ]
int index = 0;
for (int i = 0; i < 6; i++) {
if (freeCompsMap[i]) {
M[i][index] = -1.0;
index++;
}
}
index = nCols - 6;
for (int i = 0; i < 6; i++) {
if (freeCompsMap[nCols - 6 + i]) {
M[i][index] = 1.0;
index++;
}
}
j = 6;
}
for (int k : mapConstraints.keySet()) {
M[j][k] = 1.0;
j++;
}
return M;
}
/** {@inheritDoc} */
@Override
protected double[] computeAdditionalConstraints(final List<SpacecraftState> propagatedSP) {
// The additional constraint vector has the following form :
// [ xni - x1i ]----
// [ yni - x1i ] |
// Fadd(X) = [ zni - x1i ] vector's component
// [vxni - vx1i] for a closed orbit
// [vyni - vy1i] |
// [vzni - vz1i]----
// [ y1i - y1d ]---- other constraints (component of
// [ ... ] | a patch point equals to a
// [vz2i - vz2d]---- desired value)
final Map<Integer, Double> mapConstraints = getConstraintsMap();
final List<SpacecraftState> patchedSpacecraftStates = getPatchedSpacecraftState();
final double[] fxAdditional = new double[mapConstraints.size() + (isClosedOrbit ? 6 : 0)];
int i = 0;
if (isClosedOrbit) {
final AbsolutePVCoordinates apv1i = patchedSpacecraftStates.get(0).getAbsPVA();
final AbsolutePVCoordinates apvni = patchedSpacecraftStates.get(patchedSpacecraftStates.size() - 1).getAbsPVA();
fxAdditional[0] = apvni.getPosition().getX() - apv1i.getPosition().getX();
fxAdditional[1] = apvni.getPosition().getY() - apv1i.getPosition().getY();
fxAdditional[2] = apvni.getPosition().getZ() - apv1i.getPosition().getZ();
fxAdditional[3] = apvni.getVelocity().getX() - apv1i.getVelocity().getX();
fxAdditional[4] = apvni.getVelocity().getY() - apv1i.getVelocity().getY();
fxAdditional[5] = apvni.getVelocity().getZ() - apv1i.getVelocity().getZ();
i = 6;
}
// Update additional constraints
updateAdditionalConstraints(i, fxAdditional);
return fxAdditional;
}
/** {@inheritDoc} */
@Override
protected int getNumberOfConstraints() {
return super.getNumberOfConstraints() + (isClosedOrbit ? 6 : 0);
}
}