Station keeping using low thrust

Hello @dominic-mt,

I figured i would give a full Java example instead. It is slightly different from what @Emilien_mrlt did but the result is equivalent.

To answer your question, you do not really need to implement these methods as they are only used when performing operation with CalculusFieldElement (advanced usage of Orekit i would say). Please see the code below for a concrete example in Java without implementing these functions.

Otherwise you would have to implement these methods to convert input start/stop detectors into their equivalent Field class.

Logic

In my version, i decided to implement some kind of TriggeredEvent which is a wrapping event taking one raw EventDetector to be triggered and two other EventDetector to start and stop checking the g function of the raw detector. Perhaps the same things could have been done with a EventEnablingPredicateFilter but i wanted to write something simple to understand.

By doing so, the wrapped eclipse detector that is used to trigger the maneuvers is only activated when the lower SMA threshold is triggered. Once the target upper SMA threshold is reached, the wrapped eclipse detector is deactivated until lower SMA threshold is reached again.

Java code sample

BEWARE 1: Check the loadOrekitData method as your Orekit data could be located elsewhere on your system

BEWARE 2: I used a ridiculously high drag coefficient again to avoid having to propagate for weeks/months

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.events.Action;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.LofOffset;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvider;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.DragSensitive;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.forces.maneuvers.Maneuver;
import org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel;
import org.orekit.forces.maneuvers.propulsion.PropulsionModel;
import org.orekit.forces.maneuvers.trigger.ManeuverTriggers;
import org.orekit.forces.maneuvers.trigger.StartStopEventsTrigger;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.LOFType;
import org.orekit.models.earth.ReferenceEllipsoid;
import org.orekit.models.earth.atmosphere.Atmosphere;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherData;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.AbstractDetector;
import org.orekit.propagation.events.AdaptableInterval;
import org.orekit.propagation.events.BooleanDetector;
import org.orekit.propagation.events.EclipseDetector;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldAbstractDetector;
import org.orekit.propagation.events.handlers.ContinueOnEvent;
import org.orekit.propagation.events.handlers.EventHandler;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.ParameterDriver;

import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Collections;
import java.util.List;

public class StationKeepingNotInEclipse {

    public static void main(String[] args) {
        // Load orekit data
        loadOrekitData();

        // Initialize constants
        final Frame            GCRF  = FramesFactory.getGCRF();
        final Frame            ITRF  = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
        final OneAxisEllipsoid EARTH = ReferenceEllipsoid.getIers2010(ITRF);

        // Create some initial orbit
        final AbsoluteDate initialDate = new AbsoluteDate();

        final double a     = 6950e3;
        final double e     = 0.001;
        final double i     = FastMath.toRadians(30.0);
        final double omega = FastMath.toRadians(87.0);
        final double raan  = FastMath.toRadians(45.0);
        final double nu    = FastMath.toRadians(0);

        final KeplerianOrbit initialOrbit = new KeplerianOrbit(a,
                                                               e,
                                                               i,
                                                               omega,
                                                               raan,
                                                               nu,
                                                               PositionAngleType.TRUE,
                                                               FramesFactory.getGCRF(),
                                                               initialDate,
                                                               Constants.EIGEN5C_EARTH_MU);

        // Create propagator
        final NumericalPropagator numericalPropagator = getNumericalPropagator(initialOrbit);

        // Create initial state
        final SpacecraftState initialState = new SpacecraftState(initialOrbit, 150);
        numericalPropagator.setInitialState(initialState);

        // Add drag force
        final double    crossSection = 1;
        final double    dragCoef     = 20000; // Ridiculously high drag coefficient
        final DragForce dragForce    = createDragForce(crossSection, dragCoef, EARTH);
        numericalPropagator.addForceModel(dragForce);

        // Add station keeping maneuver
        // Create attitude provider override
        final AttitudeProvider attitudeProvider = new LofOffset(GCRF, LOFType.TNW);

        // Create maneuver trigger
        final double           lowSMAThreshold  = 6900e3;
        final double           highSMAThreshold = 7000e3;
        final ManeuverTriggers maneuverTriggers = getManeuverTriggers(lowSMAThreshold, highSMAThreshold, EARTH);

        // Create propulsion model
        final double thrust = 1;
        final double isp    = 300;
        final PropulsionModel propulsionModel =
                new BasicConstantThrustPropulsionModel(thrust, isp, Vector3D.PLUS_I, "station_keeping");

        final Maneuver stationKeepingManeuver = new Maneuver(attitudeProvider, maneuverTriggers, propulsionModel);
        numericalPropagator.addForceModel(stationKeepingManeuver);

        // Add eclipse with handler for better understanding of displayed values
        final EclipseDetector eclipseDetector = getEclipseDetector(EARTH).withHandler((s, detector, increasing) -> {
            if (increasing) {
                System.out.println("Exiting eclipse");
            } else {
                System.out.println("Entering eclipse");
            }
            return Action.CONTINUE;
        });
        numericalPropagator.addEventDetector(eclipseDetector);

        // Add step handler for semi major axis display & output during propagation
        final StringBuilder dataBuilder = new StringBuilder();
        dataBuilder.append("timeFromStart,a\n");
        numericalPropagator.setStepHandler(60, (s) -> {

            // Display date and current SMA
            System.out.format("Date: %s | Semi-major axis: %3.3f km\n", s.getDate(), s.getA() / 1e3);

            // Add data for output
            dataBuilder.append(s.getDate().durationFrom(initialDate));
            dataBuilder.append(",");
            dataBuilder.append(s.getA());
            dataBuilder.append("\n");
        });

        // Propagate
        numericalPropagator.propagate(initialDate.shiftedBy(2 * Constants.JULIAN_DAY));

        // Store results in csv file
        writeResults(dataBuilder);
    }

    private static void loadOrekitData() {

        //Loading Data
        final File rootFile   = new File(System.getProperty("user.home"));
        final File orekitData = new File(rootFile, "orekit-data");
        if (!orekitData.exists()) {
            System.out.format("File %s not found.%n", orekitData.getAbsolutePath());
        }
        final DataProvider dirCrawler = new DirectoryCrawler(orekitData);
        DataContext.getDefault().getDataProvidersManager().addProvider(dirCrawler);
    }

    private static NumericalPropagator getNumericalPropagator(Orbit initialOrbit) throws OrekitException {

        // Create integrator
        double minStep           = 0.001;
        double maxStep           = 60.;
        double positionTolerance = 1e-3;
        double[][] tolerances =
                NumericalPropagator.tolerances(positionTolerance, initialOrbit, initialOrbit.getType());
        DormandPrince853Integrator integrator =
                new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);

        // Create propagator
        return new NumericalPropagator(integrator);
    }

    private static DragForce createDragForce(final double crossSection,
                                             final double dragCoef,
                                             final BodyShape earth) {

        // Create atmosphere
        final Atmosphere atmosphere =
                new NRLMSISE00(new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES),
                               CelestialBodyFactory.getSun(), earth);

        // Create drag sensitive
        final DragSensitive dragSensitive = new IsotropicDrag(crossSection, dragCoef);


        return new DragForce(atmosphere, dragSensitive);
    }

    private static ManeuverTriggers getManeuverTriggers(final double lowSMAThreshold,
                                                        final double highSMAThreshold,
                                                        final OneAxisEllipsoid earth) {
        // Create low and high sma thresholds
        final EventDetector smaLowCrossingDetector =
                BooleanDetector.notCombine(new SemiMajorAxisCrossingDetector(lowSMAThreshold));
        final EventDetector smaHighCrossingDetector = new SemiMajorAxisCrossingDetector(highSMAThreshold);

        // Create eclipse detector
        final EventDetector rawEclipseDetector = getEclipseDetector(earth);

        // Create triggered eclipse detector
        final TriggeredEvent triggeredEclipseDetector =
                new TriggeredEvent(rawEclipseDetector, smaLowCrossingDetector, smaHighCrossingDetector);

        // Create stop event trigger
        final EventDetector reversedEclipseDetector =
                BooleanDetector.notCombine(getEclipseDetector(earth));

        final BooleanDetector stopEventTrigger =
                BooleanDetector.orCombine(smaHighCrossingDetector, reversedEclipseDetector);

        return new StationKeepingTrigger(triggeredEclipseDetector, stopEventTrigger);
    }

    private static EclipseDetector getEclipseDetector(OneAxisEllipsoid earth) {
        return new EclipseDetector(CelestialBodyFactory.getSun(),
                                   Constants.SUN_RADIUS,
                                   earth);
    }

    private static void writeResults(StringBuilder dataBuilder) {
        final File root   = new File(System.getProperty("user.home"));
        final File output = new File(root, "output.csv");
        try {
            final FileWriter fileWriter = new FileWriter(output);
            fileWriter.write(dataBuilder.toString());
            fileWriter.close();
        } catch (IOException e) {
            throw new RuntimeException(e);
        }
    }

    private static class SemiMajorAxisCrossingDetector extends AbstractDetector<SemiMajorAxisCrossingDetector> {

        private final double threshold;

        public SemiMajorAxisCrossingDetector(double threshold) {
            super(AbstractDetector.DEFAULT_MAXCHECK,
                  AbstractDetector.DEFAULT_THRESHOLD,
                  AbstractDetector.DEFAULT_MAX_ITER,
                  new ContinueOnEvent());
            this.threshold = threshold;
        }

        public SemiMajorAxisCrossingDetector(double maxCheck,
                                             double threshold,
                                             int maxIter,
                                             EventHandler handler,
                                             double smaThreshold) {
            super(maxCheck, threshold, maxIter, handler);
            this.threshold = smaThreshold;
        }

        @Override
        protected SemiMajorAxisCrossingDetector create(AdaptableInterval newMaxCheck,
                                                       double newThreshold,
                                                       int newMaxIter,
                                                       EventHandler newHandler) {
            return new SemiMajorAxisCrossingDetector(newMaxCheck.currentInterval(null),
                                                     newThreshold, newMaxIter, newHandler, threshold);
        }

        @Override
        public double g(SpacecraftState s) {
            return s.getA() - threshold;
        }
    }

    private static class TriggeredEvent extends AbstractDetector<TriggeredEvent> {

        private final EventDetector triggeredDetector;
        private final EventDetector startTriggerDetector;
        private final EventDetector stopTriggerDetector;

        private boolean isTriggered = false;

        public TriggeredEvent(EventDetector triggeredDetector,
                              EventDetector startTriggerDetector,
                              EventDetector stopTriggerDetector) {
            super(AbstractDetector.DEFAULT_MAXCHECK,
                  AbstractDetector.DEFAULT_THRESHOLD,
                  AbstractDetector.DEFAULT_MAX_ITER,
                  new ContinueOnEvent());
            this.triggeredDetector    = triggeredDetector;
            this.startTriggerDetector = startTriggerDetector;
            this.stopTriggerDetector  = stopTriggerDetector;
        }

        public TriggeredEvent(double maxCheck,
                              double threshold,
                              int maxIter,
                              EventHandler handler,
                              EventDetector triggeredDetector,
                              EventDetector startTriggerDetector,
                              EventDetector stopTriggerDetector) {
            super(maxCheck, threshold, maxIter, handler);
            this.triggeredDetector    = triggeredDetector;
            this.startTriggerDetector = startTriggerDetector;
            this.stopTriggerDetector  = stopTriggerDetector;
        }

        @Override
        protected TriggeredEvent create(AdaptableInterval newMaxCheck,
                                        double newThreshold,
                                        int newMaxIter,
                                        EventHandler newHandler) {
            return new TriggeredEvent(newMaxCheck.currentInterval(null),
                                      newThreshold,
                                      newMaxIter,
                                      newHandler,
                                      triggeredDetector,
                                      startTriggerDetector,
                                      stopTriggerDetector);
        }

        @Override
        public double g(SpacecraftState s) {

            if (startTriggerDetector.g(s) > 0) {
                isTriggered = true;
            } else if (stopTriggerDetector.g(s) > 0) {
                isTriggered = false;
            }

            if (isTriggered) {
                return triggeredDetector.g(s);
            }

            return -1;
        }
    }

    private static class StationKeepingTrigger extends StartStopEventsTrigger<TriggeredEvent, BooleanDetector> {

        public StationKeepingTrigger(TriggeredEvent prototypeStartDetector,
                                     BooleanDetector prototypeStopDetector) {
            super(prototypeStartDetector, prototypeStopDetector);
        }

        @Override
        public List<ParameterDriver> getParametersDrivers() {
            return Collections.emptyList();
        }

        @Override
        protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>> FieldAbstractDetector<D, S> convertStartDetector(
                Field<S> field,
                TriggeredEvent detector) {
            return null;
        }

        @Override
        protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>> FieldAbstractDetector<D, S> convertStopDetector(
                Field<S> field,
                BooleanDetector detector) {
            return null;
        }
    }
}

Example plot

2 Likes