Numerical Integration

// Define a function that returns the array of derivatives of vector y, given y and
// its dependent variable x. Since our system is one-dimensional, y has length 1, and
// so does our return value.

private double[] derivative(double x, double[] y) {
    return new double[] { Math.cos(x) };
}

// Instantiate a system by using the derivative function and specifying the order of the dependent variable (y) in the differential equation.
OrdinaryDifferentialEquationFunction function = OrdinaryDifferentialEquationFunction.of(this::derivative);
DependentVariableDifferentialEquation systemOfEquations = new DependentVariableDifferentialEquation(new OrdinaryDifferentialEquationSystem(function, 1));

// Instantiate a numerical integrator.
RungeKutta4Integrator rk4Integrator = new RungeKutta4Integrator(systemOfEquations);

RungeKuttaCashKarp45Integrator rkckIntegrator = new RungeKuttaCashKarp45Integrator(systemOfEquations, 0.1);

// Specify that we want to adapt step size to mitigate error.
rkckIntegrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);

// Specify the tolerance we want to use.
rkckIntegrator.setAbsoluteTolerance(Constants.Epsilon10);

rkckIntegrator.initialize(0.1, new double[] { 0.0 });

while (rkckIntegrator.getFinalIndependentVariableValue() < 4.7) {
    rkckIntegrator.integrate();
}

// Since the while loop guarantees that we reached or overshot our desired endpoint (x = 4.7),
// we reintegrate with a specific step size to reach exactly 4.7. Since the step size
// here is guaranteed to be smaller than our last step (which overshot 4.7),
// our computed value is at least as accurate as the value computed for our last step.
rkckIntegrator.reintegrate(4.7 - rkckIntegrator.getInitialIndependentVariableValue());
y = rkckIntegrator.getFinalDependentVariableValues();

private double[] piecewiseDerivative(double x, double[] y) {
    if ((int) x % 2 == 0) {
        return new double[] { 10.0 };
    }
    return new double[] { -10.0 };
}

// Define a function that returns the array of derivatives of vector y, given y and
// its dependent variable x.

private double[] derivatives(double x, double[] y) {
    return new double[] { y[0] * 1e-5 * x, y[1] * 1e-5 * x, y[2] * 1e-5 * x };
}

// Here we specify the order of each dependent variable (i.e. each component of y) in the differential equation.
function = OrdinaryDifferentialEquationFunction.of(this::derivatives);
systemOfEquations = new DependentVariableDifferentialEquation(new OrdinaryDifferentialEquationSystem(function, 1, 1, 1));

// Instantiate the integrator.
RungeKuttaFehlberg78Integrator rkfIntegrator = new RungeKuttaFehlberg78Integrator();
rkfIntegrator.setInitialStepSize(1);
rkfIntegrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);
rkfIntegrator.setSystemOfEquations(systemOfEquations);

// Specify the initial conditions for the differential equation.
rkfIntegrator.initialize(0.1, new double[] { -1, 0.1, 3.124 });

while (rkfIntegrator.getFinalIndependentVariableValue() < 61.1) {
    rkfIntegrator.integrate();
}
rkfIntegrator.reintegrate(61.1 - rkfIntegrator.getInitialIndependentVariableValue());
y = rkfIntegrator.getFinalDependentVariableValues();

// Define a function that returns the first derivative of the first component of y and the
// second derivative of the second component of y, given y and its dependent variable x.
// Note the format of y: y1, y1', y2. In general, a dependent variable's derivatives' values
// immediately follow its own value.

private double[] secondOrderDerivatives(double x, double[] y) {
    return new double[] { -g, y[2] * 1e-5 * x };
}

// Here we specify that y's first component has order 2, and its second component has order 1.
function = OrdinaryDifferentialEquationFunction.of(this::secondOrderDerivatives);
systemOfEquations = new DependentVariableDifferentialEquation(new OrdinaryDifferentialEquationSystem(function, 2, 1));

// Instantiate the integrator.
rk4Integrator = new RungeKutta4Integrator(systemOfEquations);

// Specify initial conditions.
rk4Integrator.initialize(0.0, new double[] { 0.0, g, 0.1 });
rk4Integrator.integrate();
rk4Integrator.integrate();

y = rk4Integrator.getFinalDependentVariableValues();

private static final JulianDate s_startDate = GregorianDate.parse("20000831T124242.00").toJulianDate();
private PointEvaluator m_pointEvaluator;
private ForceEvaluator m_gravityEvaluator;
private ForceEvaluator m_thrustEvaluator;
private ScalarEvaluator m_massEvaluator;
private static final double StepSize = 0.1;

private final void generateEvaluatorsFromPropagation() {
    CentralBody earth = CentralBodiesFacet.getFromContext().getEarth();

    // Create a newtonian point to store force models for the numerical propagator.
    PropagationNewtonianPoint newtonianPoint = new PropagationNewtonianPoint();
    newtonianPoint.setIdentification("myPoint");

    // Initialize the point approximately at the Earth's surface, at rest.
    newtonianPoint.setInitialPosition(new Cartesian(0.0, 0.0, 6.37e6));
    newtonianPoint.setInitialVelocity(new Cartesian(0.0, 0.0, 0.0));

    // Set the point's mass to a custom scalar that decreases with time, representing the rocket's fuel depletion.
    newtonianPoint.setMass(new MassScalar());

    ForceModel gravityModel = new TwoBodyGravity(newtonianPoint.getIntegrationPoint());
    ForceModel thrustModel = new ConstantForce(earth.getFixedFrame().getAxes(), new Cartesian(0.0, 0.0, 30e6));

    newtonianPoint.getAppliedForces().add(gravityModel);
    newtonianPoint.getAppliedForces().add(thrustModel);

    // Create a propagator and get a motion collection on a 20 second interval.
    NumericalPropagatorDefinition def = new NumericalPropagatorDefinition();
    final RungeKutta4Integrator rk4 = new RungeKutta4Integrator();
    rk4.setInitialStepSize(StepSize);
    def.setIntegrator(rk4);

    def.setEpoch(s_startDate);
    def.getIntegrationElements().add(newtonianPoint);

    NumericalPropagator propagator = def.createPropagator();
    NumericalPropagationStateHistory history = propagator.propagate(new Duration(0, 20.0), 1);
    DateMotionCollection1<Cartesian> motions = history.getDateMotionCollection("myPoint");

    // Get a point representing the propagation that we can use for later analysis.
    Point point = new PointInterpolator(earth.getFixedFrame(), InterpolationAlgorithmType.LINEAR, 1, motions);
    m_pointEvaluator = point.getEvaluator();

    // Get evaluators for later analysis.
    m_gravityEvaluator = new TwoBodyGravity(point).getForceEvaluator();
    m_thrustEvaluator = thrustModel.getForceEvaluator();
    m_massEvaluator = new MassScalar().getEvaluator();
}

// Recall that the derivative of work is power, i.e. force times velocity.
private double[] thrustPower(double t, double[] w) {
    JulianDate evalDate = s_startDate.add(new Duration(0, t));
    return new double[] { m_thrustEvaluator.evaluate(evalDate).dot(m_pointEvaluator.evaluate(evalDate, 1).getFirstDerivative()) };
}

generateEvaluatorsFromPropagation();

OrdinaryDifferentialEquationFunction function = OrdinaryDifferentialEquationFunction.of(this::thrustPower);
DependentVariableDifferentialEquation systemOfEquations = new DependentVariableDifferentialEquation(new OrdinaryDifferentialEquationSystem(function, 1));

RungeKutta4Integrator rk4Integrator = new RungeKutta4Integrator(systemOfEquations);
rk4Integrator.setInitialStepSize(StepSize);

// At t=0, the rocket hasn't done any work yet.
rk4Integrator.initialize(0.0, new double[] { 0.0 });

// Populate a list of work values along a 10 second interval.
ArrayList<Double> thrustWorks = new ArrayList<Double>();
for (int i = 0; i < 100; i++) {
    rk4Integrator.integrate();
    double[] w = rk4Integrator.getFinalDependentVariableValues();
    thrustWorks.add(w[0]);
}

private double[] gravityPower(double t, double[] w) {
    JulianDate evalDate = s_startDate.add(new Duration(0, t));
    // TwoBodyGravity is defined as a specific force, so we must multiply by mass to get the force in Newtons.
    return new double[] { m_massEvaluator.evaluate(evalDate) * m_gravityEvaluator.evaluate(evalDate).dot(m_pointEvaluator.evaluate(evalDate, 1).getFirstDerivative()) };
}

function = OrdinaryDifferentialEquationFunction.of(this::gravityPower);
systemOfEquations = new DependentVariableDifferentialEquation(new OrdinaryDifferentialEquationSystem(function, 1));
rk4Integrator = new RungeKutta4Integrator(systemOfEquations);
rk4Integrator.setInitialStepSize(StepSize);

rk4Integrator.initialize(0.0, new double[] { 0.0 });

ArrayList<Double> gravityWorks = new ArrayList<Double>();
for (int i = 0; i < 100; i++) {
    rk4Integrator.integrate();
    double[] w = rk4Integrator.getFinalDependentVariableValues();
    gravityWorks.add(w[0]);
}