Torque Models

// Create a new type derived from TorqueModel. This type represents a time-varying torque
// that is equal to the InitialTorque at time Epoch but varies linearly over time
// at a rate equal to TorqueRate.
public static class LinearTorque extends TorqueModel {
    public LinearTorque(Cartesian initialTorque, Cartesian torqueRate, JulianDate epoch, Axes referenceAxes) {
        m_initialTorque = initialTorque;
        m_torqueRate = torqueRate;
        m_epoch = epoch;
        m_referenceAxes = referenceAxes;
    }

    private LinearTorque(LinearTorque existingInstance, CopyContext context) {
        super(existingInstance, context);
        m_initialTorque = existingInstance.m_initialTorque;
        m_torqueRate = existingInstance.m_torqueRate;
        m_epoch = existingInstance.m_epoch;
        m_referenceAxes = context.updateReference(existingInstance.m_referenceAxes);
    }

    // This is called to make a copy of the object, which it does by calling the
    // copy constructor above.
    @Override
    public Object clone(CopyContext context) {
        return new LinearTorque(this, context);
    }

    // This method is only called by the IsSameDefinition method in the base class to
    // determine if two torque models are equivalent. Derived classes MUST override this method and
    // check all new fields introduced by the derived class for definitional equivalence. It
    // is NOT necessary to check base class fields because the base class will already have
    // done that.
    @Override
    protected final boolean checkForSameDefinition(TorqueModel other) {
        LinearTorque o = other instanceof LinearTorque ? (LinearTorque) other : null;
        return o != null 
                && areSameDefinition(m_initialTorque, o.m_initialTorque) 
                && areSameDefinition(m_torqueRate, o.m_torqueRate) 
                && areSameDefinition(m_epoch, o.m_epoch)
                && areSameDefinition(m_referenceAxes, o.m_referenceAxes) 
                && checkForSameDefinition(o);
    }

    // This method allows all types derived from LinearTorque to continue the process
    // of making unique CheckForSameDefinition methods.
    protected boolean checkForSameDefinition(LinearTorque other) {
        return other != null && other.getClass().equals(LinearTorque.class);
    }

    // Called to determine a hash code for the current configuration of this object.
    // Derived classes MUST override this method and compute a hash code that combines: 
    // a unique hash code seed, the base implementation result, and the 
    // hash codes of all new fields introduced by the derived class which are used 
    // in the CheckForSameDefinition method.
    @Override
    protected int computeCurrentDefinitionHashCode() {
        return HashCode.combine(LinearTorque.class.hashCode(),
                                super.computeCurrentDefinitionHashCode(), 
                                getDefinitionHashCode(m_initialTorque),
                                getDefinitionHashCode(m_torqueRate), 
                                getDefinitionHashCode(m_epoch),
                                getDefinitionHashCode(m_referenceAxes));
    }

    // Called to enumerate all of the other objects on which this object depends.  This
    // allows clients to navigate the graph of objects related to a computation.
    @Override
    public void enumerateDependencies(DependencyEnumerator enumerator) {
        super.enumerateDependencies(enumerator);
        enumerator.enumerate(m_referenceAxes);
    }

    // The time in which the torque is equal to the InitialTorque.
    public final JulianDate getEpoch() {
        return m_epoch;
    }

    // The time in which the torque is equal to the InitialTorque.
    public final void setEpoch(JulianDate value) {
        // Called to ensure that properties cannot be set if
        // the object is frozen.
        throwIfFrozen();
        m_epoch = value;
    }

    // The torque vector at the Epoch.
    public final Cartesian getInitialTorque() {
        return m_initialTorque;
    }

    // The torque vector at the Epoch.
    public final void setInitialTorque(Cartesian value) {
        throwIfFrozen();
        m_initialTorque = value;
    }

    // The rate at which the torque changes over time.
    public final Cartesian getTorqueRate() {
        return m_torqueRate;
    }

    // The rate at which the torque changes over time.
    public final void setTorqueRate(Cartesian value) {
        throwIfFrozen();
        m_torqueRate = value;
    }

    // The reference axes that the torque is defined relative to.
    // These are usually the body axes of the rigid body.
    public final Axes getReferenceAxes() {
        return m_referenceAxes;
    }

    // The reference axes that the torque is defined relative to.
    // These are usually the body axes of the rigid body.
    public final void setReferenceAxes(Axes value) {
        throwIfFrozen();
        m_referenceAxes = value;
    }

    // This method is responsible for returning an instance of the private
    // LinearTorqueEvaluator class. It should ensure that the properties are not null or
    // in an invalid state, and then use the evaluator group when it constructs
    // and returns an instance of the LinearTorqueEvaluator.
    @Override
    public TorqueEvaluator getTorqueEvaluator(EvaluatorGroup group) {
        if (group == null) {
            throw new ArgumentNullException("group");
        }
        if (getReferenceAxes() == null) {
            throw new PropertyInvalidException("ReferenceAxes", PropertyInvalidException.PropertyCannotBeNull);
        }

        return group.createEvaluator(getEvaluatorCallback);
    }

    // This method, which is passed to the evaluator group in the method above as a delegate,
    // will only be called by the delegate if the evaluator does not yet exist in the group
    // and needs to be created.
    private final LinearTorqueEvaluator createEvaluator(EvaluatorGroup group) {
        return new LinearTorqueEvaluator(group, m_epoch, m_initialTorque, m_torqueRate, m_referenceAxes);
    }

    // This method is used by PropagationEulerianAxes to combine all the torques
    // that are applied upon the rigid body.
    @Override
    public void buildTorqueEvaluator(ResultantTorqueBuilder builder, EvaluatorGroup group) {
        if (builder == null) {
            throw new ArgumentNullException("builder");
        }
        if (group == null) {
            throw new ArgumentNullException("group");
        }
        builder.addTorque(getTorqueEvaluator(group));
    }

    private final EvaluatorGroup.Callback0<TorqueEvaluator> getEvaluatorCallback =
            EvaluatorGroup.Callback0.of(this::createEvaluator);
    private Cartesian m_initialTorque = new Cartesian();
    private Cartesian m_torqueRate = new Cartesian();
    private JulianDate m_epoch = new JulianDate();
    private Axes m_referenceAxes;

    // This is the definition of the Evaluator that is used to actually evaluate this TorqueModel.
    // Because it is a private, nested class, it is not visible outside of the
    // LinearTorque class. This is ok, though, because once it is created users only
    // interact with it via a reference to its base class: TorqueEvaluator.
    private static final class LinearTorqueEvaluator extends TorqueEvaluator {
        public LinearTorqueEvaluator(EvaluatorGroup group, JulianDate epoch, 
                Cartesian initialTorque, Cartesian torqueRate, Axes referenceAxes) {
            super(group);
            m_initialTorque = initialTorque;
            m_torqueRate = torqueRate;
            m_epoch = epoch;
            m_referenceAxes = referenceAxes;
        }

        private LinearTorqueEvaluator(LinearTorqueEvaluator existingInstance, CopyContext context) {
            super(existingInstance, context);
            m_initialTorque = existingInstance.m_initialTorque;
            m_torqueRate = existingInstance.m_torqueRate;
            m_epoch = existingInstance.m_epoch;

            // Use the UpdateReference method on the CopyContext to perform any updates to the
            // reference that are required by the context.
            m_referenceAxes = context.updateReference(existingInstance.m_referenceAxes);
        }

        // This method is used by the EvaluatorGroup system to avoid redundant evaluations. The
        // EvaluatorGroup may call it on your evaluator in order to replace your evaluator's
        // reference to another evaluator with a reference to a version that caches its last
        // result.
        // 
        // In this example, there are no internal evaluators so the method does nothing.
        @Override
        public void updateEvaluatorReferences(CopyContext context) {}

        // The Clone method should call the copy constructor.
        @Override
        public Object clone(CopyContext context) {
            return new LinearTorqueEvaluator(this, context);
        }

        // Override the Dispose method to call the Dispose() method on any nested
        // evaluators or other disposable nested types.
        @Override
        protected void dispose(boolean disposing) {}

        // This method returns a collection of time intervals when data is 
        // available from this Evaluator.
        @Override
        public TimeIntervalCollection getAvailabilityIntervals(TimeIntervalCollection consideredIntervals) {
            return consideredIntervals;
        }

        // This method determines if there is data available from this Evaluator at 
        // the specified date.
        @Override
        public boolean isAvailable(JulianDate date) {
            return true;
        }

        // This property determines if this Evaluator result changes depending on the time at
        // which it is evaluated.
        @Override
        public boolean getIsTimeVarying() {
            // In this example, the torque is time-varying as long as the torque rate is not zero.
            return Cartesian.notEquals(m_torqueRate, Cartesian.getZero());
        }

        // This property determines if this Evaluator can safely be used from multiple threads
        // simultaneously. If the evaluator stores data during the Evaluate call, it is not thread
        // safe. Otherwise, it generally is thread safe as long as any nested evaluators it uses
        // are thread safe.
        @Override
        public boolean getIsThreadSafe() {
            return true;
        }

        // This is where we do the actual evaluation of the torque at a particular time.
        @Override
        public Cartesian evaluate(JulianDate date) {
            // The initial torque vector is added to the torque rate multiplied
            // by the time elapsed since the reference epoch to form the linear
            // torque model.
            return Cartesian.add(m_initialTorque, Cartesian.multiply(m_epoch.secondsDifference(date), m_torqueRate));
        }

        // The evaluator is defined to produce torques that are defined
        // with respect to the reference axes, which is usually the body
        // axes of the rigid body.
        @Override
        public Axes getDefinedIn() {
            return m_referenceAxes;
        }

        private Cartesian m_initialTorque = new Cartesian();
        private Cartesian m_torqueRate = new Cartesian();
        private JulianDate m_epoch = new JulianDate();
        private Axes m_referenceAxes;
    }
}

// The first step of creating a NumericalPropagator is to create and populate
// a NumericalPropagatorDefinition.
NumericalPropagatorDefinition definition = new NumericalPropagatorDefinition();

// Epoch is September 5, 2022 at 11:30:00.000 UTC.
JulianDate epoch = new GregorianDate(2022, 9, 5, 11, 30, 0.0).toJulianDate();
definition.setEpoch(epoch);

// A PropagationNewtonianPoint is defined to provide the position of a
// gravity-gradient satellite in an equatorial circular orbit with an
// altitude of 300 km.
KeplerianElements elements = new KeplerianElements(WorldGeodeticSystem1984.SemimajorAxis + 300000.0,
        0.0, 0.0, 0.0, 0.0, 0.0, WorldGeodeticSystem1984.GravitationalParameter);
Motion1<Cartesian> initialOrbitalState = elements.toCartesian();

PropagationNewtonianPoint point = new PropagationNewtonianPoint();
point.setInitialPosition(initialOrbitalState.getValue());
point.setInitialVelocity(initialOrbitalState.getFirstDerivative());
point.getAppliedForces().add(new TwoBodyGravity(point.getIntegrationPoint()));
point.setMass(Scalar.toScalar(100D)); // kg

// The attitude of the spacecraft is defined by a PropagationEulerianAxes.
Cartesian initialAngularVelocity = new Cartesian(0.05, 0.0, 0.01); // rad/s
UnitQuaternion initialAttitudeQuaternion = UnitQuaternion.getIdentity(); // The body frame and inertial frame are initially aligned.
final double transverseMomentOfInertia = 200.0; // kg*m^2
final double spinMomentOfInertia = 100.0; // kg*m^2
Matrix3By3Symmetric axisymmetricInertiaMatrix = Matrix3By3Symmetric.diagonalMatrix(transverseMomentOfInertia, transverseMomentOfInertia, spinMomentOfInertia);

PropagationEulerianAxes attitude = new PropagationEulerianAxes();
attitude.setInitialAttitudeQuaternion(initialAttitudeQuaternion);
attitude.setInitialAngularVelocity(initialAngularVelocity);
attitude.setInertiaMatrix(axisymmetricInertiaMatrix);

// Add linear torque, magnetic coil torque, and gravity-gradient torque to PropagationEulerianAxes.
LinearTorque linearTorque = new LinearTorque(Cartesian.getZero(), UnitCartesian.multiply(0.00000001, UnitCartesian.getUnitZ()), epoch, attitude.getIntegrationAxes());

// The number of loops is the number of times the magnetic coil winds around its axis.
final double numberOfLoops = 3.0;

// The electric current is set to be constant at 10 amps.
Scalar current = new ScalarFixed(10.0);

// The directed area is 1.0 square meters in the body Z-axis. This means
// that the current is flowing through the loops in the counter-clockwise
// direction in the +Z axis of the body frame.
Cartesian directedArea = Cartesian.toCartesian(UnitCartesian.getUnitZ());
Vector directedAreaVector = new VectorFixed(attitude.getIntegrationAxes(), directedArea);

// The WMM model can be used to define Earth's magnetic field.
MagneticCoilTorque magneticCoilTorque = MagneticCoilTorque.createMagneticCoilTorqueUsingWmmModel(wmmPath, point.getIntegrationPoint(), numberOfLoops, current, directedAreaVector,
        attitude.getIntegrationAxes());

// The IGRF model can also be used. It can be defined by either a .txt file or a .COF file.

GravityGradientTorque gravityGradientTorque = new GravityGradientTorque(point.getIntegrationPoint(), axisymmetricInertiaMatrix, attitude.getIntegrationAxes());
attitude.getAppliedTorques().add(linearTorque);
attitude.getAppliedTorques().add(magneticCoilTorque);
attitude.getAppliedTorques().add(gravityGradientTorque);

// Add point and attitude to NumericalPropagatorDefinition.
definition.getIntegrationElements().add(point);
definition.getIntegrationElements().add(attitude);

// Define integrator.
RungeKutta4Integrator integrator = new RungeKutta4Integrator();
integrator.setInitialStepSize(1.0);
integrator.setDirection(IntegrationSense.INCREASING);
definition.setIntegrator(integrator);

// Create a propagator and propagate for 1 day.
NumericalPropagator propagator = definition.createPropagator();
NumericalPropagationStateHistory history = propagator.propagate(new Duration(1, 0.0), 1);

// Extract results for the location and attitude of the gravity-gradient satellite.
DateMotionCollection1<Cartesian> locationHistory 
    = history.getDateMotionCollection(point.getIdentification());
DateMotionCollection1<RotationVectorAngularVelocity> rotationHistory 
    = history.getDateMotionCollection(attitude.getIdentification());

// Transform rotation history to a more usable form.
DateMotionCollection2<UnitQuaternion, Cartesian> transformedRotationHistory 
    = RotationVectorAngularVelocity.toDateMotionCollectionUnitQuaternionCartesian(rotationHistory, 1);
// Extract the final Motion value from this collection.
Motion2<UnitQuaternion, Cartesian> finalMotion 
    = transformedRotationHistory.getMotions().get(transformedRotationHistory.getCount() - 1);
UnitQuaternion finalAttitude = finalMotion.getValue();
Cartesian finalAngularVelocity = finalMotion.getFirstDerivative(); // rad/s, in inertial coordinates.