DSST Polynomial Acceleration for low thrust estimation

Hi everyone,

I tried to implement a DSST version of the polynomial acceleration force model, in order to use it in the context of orbit determination and maneuver estimation. I ran a simple test, where the acceleration is constant and low (see below). The mean element rate vector looks like:

7.395132202849769E17 -5282237.267737973 -528223.7267771859 1.1221496749908623E-10 -4.1975607666515704E-11 -3.5241106052380067E-7

Some debugging showed that in the polynomial acceleration model the gravitational constant is used instead of the acceleration value to compute the amplitude.
I could not find a (clean) way of setting things up so that the parameters are in proper order for both the getMeanElementsRates() method and the PolynomialParametricAcceleration.signedAmplitude() method.

A very dirty workaround consisting in redefining the parameters vector in the constructor of the AbstractGaussianContribution.IntegrableFucntion class allows to compute the correct derivatives (I actually only checked the first value).

0.4638186158939421 -3.312990100977405E-12 -3.312990100989179E-13 3.6869891406591825E-29 -6.08131982587497E-30 1.9052155422727878E-25

I don’t know if there are errors in my test or if it is just a wrong way of using the PolynomialParametricAcceleration class in combination with AbstractGaussianConstribution.

Further, any advice of better ways of using orekit existing force models for thrust estimation during OD are of course appreciated.

In any case I’ll be grateful for your help.
Cheers,
Andrea

Test code:

public class TestPolynomialAccelerationDSST {

private static final double SPACECRAFT_MASS = 20; // kg
private static final double THRUST_MODULUS = 5E-3; // N
private static final double THRESHOLD = 1E-15;
private static final String OREKIT_DATA_MASTER_PATH = "  "; // add your path here

public class DSSTPolynomialAcceleration extends AbstractGaussianContribution {

    private ForceModel contribution;

    protected DSSTPolynomialAcceleration (String coefficientsKeyPrefix, double threshold,
            ForceModel contribution, double mu) {
        super(coefficientsKeyPrefix, threshold, contribution, mu);
        this.contribution = contribution;
    }

    @Override
    public <T extends RealFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(
            Field<T> field) {
        // TODO Auto-generated method stub
        return null;
    }

    @Override
    public EventDetector[] getEventsDetectors() {
        // TODO Auto-generated method stub
        return null;
    }

    @Override
    public ParameterDriver[] getParametersDriversWithoutMu() {

        return contribution.getParametersDrivers();
    }

    @Override
    protected <T extends RealFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
            FieldAuxiliaryElements<T> auxiliaryElements) {
        // TODO Auto-generated method stub
        return null;
    }

    @Override
    protected double[] getLLimits(SpacecraftState state, AuxiliaryElements auxiliaryElements) {
        return new double[] { 0, 2 * FastMath.PI };
    }
}

@Before
public void setup() {
    DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
    manager.addProvider(new DirectoryCrawler(new File(OREKIT_DATA_MASTER_PATH)));
}

@Test
public void testAcceleration() {

    AbsoluteDate referenceDate = new AbsoluteDate(2021, 01, 01, TimeScalesFactory.getUTC());

    double acc = THRUST_MODULUS / SPACECRAFT_MASS;

    // Initial condition
    EquinoctialOrbit orbit = new EquinoctialOrbit(7000.e3, 1e-4, 1e-5, 1e-4, 1e-7, 0,
            PositionAngle.MEAN, FramesFactory.getEME2000(), referenceDate,
            Constants.EGM96_EARTH_MU);

    // Numerical force model
    PolynomialParametricAcceleration polynomialAcceleration = new PolynomialParametricAcceleration(
            new Vector3D(new double[] { 1, 0, 0 }),
            new LofOffset(orbit.getFrame(), LOFType.TNW), "polynomial", referenceDate, 0);
    polynomialAcceleration.getParametersDrivers()[0].setValue(acc);

    DSSTPolynomialAcceleration dsstPolynomialAcc = new DSSTPolynomialAcceleration("polyDSST",
            THRESHOLD, polynomialAcceleration, Constants.EGM96_EARTH_MU);
    dsstPolynomialAcc.registerAttitudeProvider(new LofOffset(orbit.getFrame(), LOFType.TNW));
    
    // Initialize model parameters
    dsstPolynomialAcc.initialize(new AuxiliaryElements(orbit, 1), PropagationType.MEAN,
            dsstPolynomialAcc.getParameters());

    SpacecraftState initialState = new SpacecraftState(orbit, SPACECRAFT_MASS);

    double[] parameters = dsstPolynomialAcc.getParameters();
    double[] rates = dsstPolynomialAcc.getMeanElementRate(initialState,
            new AuxiliaryElements(orbit, 1), parameters);

    System.out.println(rates[0] + " " + rates[1] + " " + rates[2] + " " + rates[3] + " "
            + rates[4] + " " + rates[5]);
}

}