ZeisModel.java
/* Copyright 2022 Bryan Cazabonne
* 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.
* Bryan Cazabonne 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.semianalytical.dsst.forces;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.util.MathArrays;
import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
/**
* Zeis model for J2-squared second-order terms.
*
* @see "ZEIS, Eric and CEFOLA, P. Computerized algebraic utilities for the
* construction of nonsingular satellite theories. Journal of Guidance and
* Control, 1980, vol. 3, no 1, p. 48-54."
*
* @see "SAN-JUAN, Juan F., LĂ“PEZ, Rosario, et CEFOLA, Paul J. A Second-Order
* Closed-Form $$ J_2 $$ Model for the Draper Semi-Analytical Satellite
* Theory. The Journal of the Astronautical Sciences, 2022, p. 1-27."
*
* @author Bryan Cazabonne
* @since 12.0
*/
public class ZeisModel implements J2SquaredModel {
/**
* Retrograde factor I.
* <p>
* DSST model needs equinoctial orbit as internal representation. Classical
* equinoctial elements have discontinuities when inclination is close to zero.
* In this representation, I = +1. <br>
* To avoid this discontinuity, another representation exists and equinoctial
* elements can be expressed in a different way, called "retrograde" orbit. This
* implies I = -1. <br>
* As Orekit doesn't implement the retrograde orbit, I is always set to +1. But
* for the sake of consistency with the theory, the retrograde factor has been
* kept in the formulas.
* </p>
*/
private static final int I = 1;
/** Constructor. */
public ZeisModel() {
// Nothing to do...
}
/** {@inheritDoc}. */
@Override
public double[] computeMeanEquinoctialSecondOrderTerms(final DSSTJ2SquaredClosedFormContext context) {
// Auxiliary elements
final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
// Zeis constant
final double c2z = computeC2Z(context);
// Useful terms
final double s2mf = 19.0 * context.getS2() - 15.0;
final double s2pIcmo = context.getS2() + I * context.getC() - 1.0;
final double s4mts2 = 19.0 * context.getS2() * context.getS2() - 30.0 * context.getS2() + 12.0;
// Second-order terms (Ref [2] Eq. 37)
final double deltaA = 0.0;
final double deltaK = -c2z * auxiliaryElements.getH() * s2mf * s2pIcmo;
final double deltaH = c2z * auxiliaryElements.getK() * s2mf * s2pIcmo;
final double deltaQ = -c2z * context.getC() * auxiliaryElements.getP() * s2mf;
final double deltaP = c2z * context.getC() * auxiliaryElements.getQ() * s2mf;
final double deltaM = 0.5 * c2z * (2.0 * s2mf * s2pIcmo + 5.0 * s4mts2 * context.getEta());
// Return
return new double[] { deltaA, deltaK, deltaH, deltaQ, deltaP, deltaM };
}
/** {@inheritDoc}. */
@Override
public <T extends CalculusFieldElement<T>> T[] computeMeanEquinoctialSecondOrderTerms(final FieldDSSTJ2SquaredClosedFormContext<T> context) {
// Auxiliary elements
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
// Field
final Field<T> field = auxiliaryElements.getDate().getField();
// Zeis constant
final T c2z = computeC2Z(context);
// Useful terms
final T s2mf = context.getS2().multiply(19.0).subtract(15.0);
final T s2pIcmo = context.getS2().add(context.getC().multiply(I)).subtract(1.0);
final T s4mts2 = context.getS2().multiply(context.getS2()).multiply(19.0).subtract(context.getS2().multiply(30.0)).add(12.0);
// Second-order terms (Ref [2] Eq. 37)
final T deltaA = field.getZero();
final T deltaK = c2z.multiply(auxiliaryElements.getH()).multiply(s2mf).multiply(s2pIcmo).negate();
final T deltaH = c2z.multiply(auxiliaryElements.getK()).multiply(s2mf).multiply(s2pIcmo);
final T deltaQ = c2z.multiply(context.getC()).multiply(auxiliaryElements.getP()).multiply(s2mf).negate();
final T deltaP = c2z.multiply(context.getC()).multiply(auxiliaryElements.getQ()).multiply(s2mf);
final T deltaM = c2z.multiply(0.5).multiply(s2mf.multiply(s2pIcmo).multiply(2.0).add(s4mts2.multiply(context.getEta()).multiply(5.0)));
// Return
final T[] terms = MathArrays.buildArray(field, 6);
terms[0] = deltaA;
terms[1] = deltaK;
terms[2] = deltaH;
terms[3] = deltaQ;
terms[4] = deltaP;
terms[5] = deltaM;
return terms;
}
/**
* Get the value of the Zeis constant.
*
* @param context model context
* @return the value of the Zeis constant
*/
public double computeC2Z(final DSSTJ2SquaredClosedFormContext context) {
final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
return 0.75 * context.getAlpha4() * auxiliaryElements.getMeanMotion() / (context.getA4() * context.getEta());
}
/**
* Get the value of the Zeis constant.
*
* @param context model context
* @param <T> type of the elements
* @return the value of the Zeis constant
*/
public <T extends CalculusFieldElement<T>> T computeC2Z(final FieldDSSTJ2SquaredClosedFormContext<T> context) {
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
return auxiliaryElements.getMeanMotion().multiply(context.getAlpha4()).multiply(0.75).divide(context.getA4().multiply(context.getEta()));
}
}