Skip to content
GitLab
Commit ee2c4123 authored by Bryan Cazabonne's avatar Bryan Cazabonne
Browse files

Added Zeis model for DSST J2-squared second order terms. 

Issue #931
parent 88b993e2
Pipeline #2682 passed with stages
in 19 minutes and 4 seconds
Expand all Hide whitespace changes
Inline Side-by-side
src/changes/changes.xml
View file @ ee2c4123
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.4" date="TBD" description="TBD">
<action dev="bryan" type="add" issue="931">
Added Zeis model for DSST J2-squared second order terms.
</action>
<action dev="vincent" type="add" issue="981">
Added ability to consider LOFType as pseudo-inertial frame.
</action>
......
src/design/dsst-propagator-class-diagram.puml
View file @ ee2c4123
......@@ -62,6 +62,7 @@
DSSTForceModel -right-> AuxiliaryElements
DSSTForceModel <|-- DSSTZonal
DSSTForceModel <|-- DSSTJ2SquaredClosedForm
DSSTForceModel <|-- DSSTTesseral
DSSTForceModel <|-- DSSTThirdBody
DSSTForceModel <|-- AbstractGaussianContribution
......
src/design/dsst-zonal-contribution-class-diagram.puml
View file @ ee2c4123
......@@ -34,9 +34,13 @@
Package forces #CBDBC8 {
interface DSSTForceModel
interface ShortPeriodTerms
interface J2SquaredModel
DSSTForceModel <|-- DSSTZonal
DSSTZonal --> ZonalShortPeriodicCoefficients
ShortPeriodTerms <|-- ZonalShortPeriodicCoefficients
DSSTForceModel <|-- DSSTJ2SquaredClosedForm
J2SquaredModel "1" <--* DSSTJ2SquaredClosedForm
J2SquaredModel <|-- ZeisModel
}
package utilities #CBDBC8 {
......@@ -51,6 +55,7 @@
package forces.gravity.potential #DDEBD8 {
interface UnnormalizedSphericalHarmonicsProvider
UnnormalizedSphericalHarmonicsProvider "1" <--* DSSTZonal
UnnormalizedSphericalHarmonicsProvider "1" <--* DSSTJ2SquaredClosedForm
}
}
......
src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTJ2SquaredClosedForm.java 0 → 100644
View file @ ee2c4123
/* 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 java.util.Collections;
import java.util.List;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.util.MathArrays;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.PropagationType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldEventDetector;
import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
import org.orekit.utils.ParameterDriver;
/**
* Second order J2-squared force model.
* <p>
* The force model implements a closed-form of the J2-squared perturbation.
* The full realization of the model is based on a gaussian quadrature.
* Even if it is very accurate, a gaussian quadrature is usually time consuming.
* A closed-form is less accurate than a gaussian quadrature, but faster.
* </p>
* @author Bryan Cazabonne
*/
public class DSSTJ2SquaredClosedForm implements DSSTForceModel {
/** Model for second order terms. */
private final J2SquaredModel j2SquaredModel;
/** Gravity field to use. */
private final UnnormalizedSphericalHarmonicsProvider provider;
/**
* Constructor.
*
* @param j2SquaredModel model for second order terms
* @param provider gravity field to use
*/
public DSSTJ2SquaredClosedForm(final J2SquaredModel j2SquaredModel,
final UnnormalizedSphericalHarmonicsProvider provider) {
// Initialize fields
this.j2SquaredModel = j2SquaredModel;
this.provider = provider;
}
/** {@inheritDoc}. */
@Override
public double[] getMeanElementRate(final SpacecraftState state,
final AuxiliaryElements auxiliaryElements,
final double[] parameters) {
// Context
final DSSTJ2SquaredClosedFormContext context = new DSSTJ2SquaredClosedFormContext(auxiliaryElements, provider);
// Second-order terms
final double[] delta = j2SquaredModel.computeMeanEquinoctialSecondOrderTerms(context);
// J2
final double J2 = -provider.onDate(state.getDate()).getUnnormalizedCnm(2, 0);
final double J2SquaredOver2 = 0.5 * J2 * J2;
// Mean elements rate
final double da = 0.0;
final double dk = J2SquaredOver2 * delta[1];
final double dh = J2SquaredOver2 * delta[2];
final double dq = J2SquaredOver2 * delta[3];
final double dp = J2SquaredOver2 * delta[4];
final double dM = J2SquaredOver2 * delta[5];
// Return
return new double[] { da, dk, dh, dq, dp, dM };
}
/** {@inheritDoc}. */
@Override
public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters) {
// Field
final Field<T> field = state.getDate().getField();
// Context
final FieldDSSTJ2SquaredClosedFormContext<T> context = new FieldDSSTJ2SquaredClosedFormContext<>(auxiliaryElements, provider);
// Second-order terms
final T[] delta = j2SquaredModel.computeMeanEquinoctialSecondOrderTerms(context);
// J2
final double J2 = -provider.onDate(state.getDate().toAbsoluteDate()).getUnnormalizedCnm(2, 0);
final double J2SquaredOver2 = 0.5 * J2 * J2;
// Mean elements rate
final T da = field.getZero();
final T dk = delta[1].multiply(J2SquaredOver2);
final T dh = delta[2].multiply(J2SquaredOver2);
final T dq = delta[3].multiply(J2SquaredOver2);
final T dp = delta[4].multiply(J2SquaredOver2);
final T dM = delta[5].multiply(J2SquaredOver2);
// Return
final T[] elements = MathArrays.buildArray(field, 6);
elements[0] = da;
elements[1] = dk;
elements[2] = dh;
elements[3] = dq;
elements[4] = dp;
elements[5] = dM;
return elements;
}
/** {@inheritDoc}. */
@Override
public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
final PropagationType type,
final double[] parameters) {
// Currently, there is no short periods for J2-squared closed-form
return Collections.emptyList();
}
/** {@inheritDoc}. */
@Override
public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
final PropagationType type,
final T[] parameters) {
// Currently, there is no short periods for J2-squared closed-form
return Collections.emptyList();
}
/** {@inheritDoc}. */
@Override
public List<ParameterDriver> getParametersDrivers() {
return Collections.emptyList();
}
/** {@inheritDoc}. */
@Override
public EventDetector[] getEventsDetectors() {
return null;
}
/** {@inheritDoc}. */
@Override
public <T extends CalculusFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
return null;
}
/** {@inheritDoc}. */
@Override
public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
// Nothing is done since this contribution is not sensitive to attitude
}
/** {@inheritDoc}. */
@Override
public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {
// Currently, there is no short periods for J2-squared closed-form
}
/** {@inheritDoc}. */
@Override
@SuppressWarnings("unchecked")
public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters, final FieldSpacecraftState<T>... meanStates) {
// Currently, there is no short periods for J2-squared closed-form
}
}
src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTJ2SquaredClosedFormContext.java 0 → 100644
View file @ ee2c4123
/* 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.util.FastMath;
import org.hipparchus.util.SinCos;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
/**
* This class is a container for the common parameters used in {@link DSSTJ2SquaredClosedForm}.
* <p>
* It performs parameters initialization at each integration step for the second-order J2-squared
* contribution to the central body gravitational perturbation.
* <p>
* @author Bryan Cazabonne
*/
public class DSSTJ2SquaredClosedFormContext extends ForceModelContext {
/** Equatorial radius of the central body to the power 4. */
private final double alpha4;
/** Semi major axis to the power 4. */
private final double a4;
/** sqrt(1 - e * e). */
private final double eta;
/** Cosine of the inclination. */
private final double c;
/** Sine of the inclination to the power 2. */
private final double s2;
/**
* Simple constructor.
*
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param provider provider for spherical harmonics
*/
public DSSTJ2SquaredClosedFormContext(final AuxiliaryElements auxiliaryElements,
final UnnormalizedSphericalHarmonicsProvider provider) {
super(auxiliaryElements);
// Sine and cosine of the inclination
final double inc = auxiliaryElements.getOrbit().getI();
final SinCos scI = FastMath.sinCos(inc);
// Other parameters
this.c = scI.cos();
this.s2 = scI.sin() * scI.sin();
final double alpha2 = provider.getAe() * provider.getAe();
this.alpha4 = alpha2 * alpha2;
this.eta = FastMath.sqrt(1.0 - auxiliaryElements.getEcc() * auxiliaryElements.getEcc());
final double a2 = auxiliaryElements.getSma() * auxiliaryElements.getSma();
this.a4 = a2 * a2;
}
/**
* Get the equatorial radius of the central body to the power 4.
* @return the equatorial radius of the central body to the power 4
*/
public double getAlpha4() {
return alpha4;
}
/**
* Get the semi major axis to the power 4.
* @return the semi major axis to the power 4
*/
public double getA4() {
return a4;
}
/**
* Get the eta value.
* @return sqrt(1 - e * e)
*/
public double getEta() {
return eta;
}
/**
* Get the cosine of the inclination.
* @return the cosine of the inclination
*/
public double getC() {
return c;
}
/**
* Get the sine of the inclination to the power 2.
* @return the sine of the inclination to the power 2
*/
public double getS2() {
return s2;
}
}
src/main/java/org/orekit/propagation/semianalytical/dsst/forces/FieldDSSTJ2SquaredClosedFormContext.java 0 → 100644
View file @ ee2c4123
/* 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.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
/**
* This class is a container for the common parameters used in {@link DSSTJ2SquaredClosedForm}.
* <p>
* It performs parameters initialization at each integration step for the second-order J2-squared
* contribution to the central body gravitational perturbation.
* <p>
* @author Bryan Cazabonne
*/
public class FieldDSSTJ2SquaredClosedFormContext<T extends CalculusFieldElement<T>> extends FieldForceModelContext<T> {
/** Equatorial radius of the central body to the power 4. */
private final double alpha4;
/** Semi major axis to the power 4. */
private final T a4;
/** sqrt(1 - e * e). */
private final T eta;
/** Cosine of the inclination. */
private final T c;
/** Sine of the inclination to the power 2. */
private final T s2;
/**
* Simple constructor.
*
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param provider provider for spherical harmonics
*/
public FieldDSSTJ2SquaredClosedFormContext(final FieldAuxiliaryElements<T> auxiliaryElements,
final UnnormalizedSphericalHarmonicsProvider provider) {
super(auxiliaryElements);
// Sine and cosine of the inclination
final T inc = auxiliaryElements.getOrbit().getI();
final FieldSinCos<T> scI = FastMath.sinCos(inc);
// Other parameters
this.c = scI.cos();
this.s2 = scI.sin().multiply(scI.sin());
final double alpha2 = provider.getAe() * provider.getAe();
this.alpha4 = alpha2 * alpha2;
this.eta = FastMath.sqrt(auxiliaryElements.getEcc().multiply(auxiliaryElements.getEcc()).negate().add(1.0));
final T a2 = auxiliaryElements.getSma().multiply(auxiliaryElements.getSma());
this.a4 = a2.multiply(a2);
}
/**
* Get the equatorial radius of the central body to the power 4.
* @return the equatorial radius of the central body to the power 4
*/
public double getAlpha4() {
return alpha4;
}
/**
* Get the semi major axis to the power 4.
* @return the semi major axis to the power 4
*/
public T getA4() {
return a4;
}
/**
* Get the eta value.
* @return sqrt(1 - e * e)
*/
public T getEta() {
return eta;
}
/**
* Get the cosine of the inclination.
* @return the cosine of the inclination
*/
public T getC() {
return c;
}
/**
* Get the sine of the inclination to the power 2.
* @return the sine of the inclination to the power 2
*/
public T getS2() {
return s2;
}
}
src/main/java/org/orekit/propagation/semianalytical/dsst/forces/J2SquaredModel.java 0 → 100644
View file @ ee2c4123
/* 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;
/**
* Semi-analytical J2-squared model.
* <p>
* This interface is implemented by models providing J2-squared
* second-order terms in equinoctial elements. These terms are
* used in the computation of the closed-form J2-squared perturbation
* in semi-analytical satellite theory.
* </p>
* @see ZeisModel
* @author Bryan Cazabonne
*/
public interface J2SquaredModel {
/**
* Compute the J2-squared second-order terms in equinoctial elements.
* @param context model context
* @return the J2-squared second-order terms in equinoctial elements.
* Order must follow: [A, K, H, Q, P, M]
*/
double[] computeMeanEquinoctialSecondOrderTerms(DSSTJ2SquaredClosedFormContext context);
/**
* Compute the J2-squared second-order terms in equinoctial elements.
* @param context model context
* @param <T> type of the elements
* @return the J2-squared second-order terms in equinoctial elements.
* Order must follow: [A, K, H, Q, P, M]
*/
<T extends CalculusFieldElement<T>> T[] computeMeanEquinoctialSecondOrderTerms(FieldDSSTJ2SquaredClosedFormContext<T> context);
}
src/main/java/org/orekit/propagation/semianalytical/dsst/forces/ZeisModel.java 0 → 100644
View file @ ee2c4123
/* 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
*/
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.