1   /* Copyright 2022-2025 Romain Serra
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.control.indirect.adjoint.cost;
18  
19  
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.Field;
22  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23  import org.orekit.propagation.FieldSpacecraftState;
24  import org.orekit.propagation.events.EventDetectionSettings;
25  import org.orekit.propagation.events.FieldEventDetectionSettings;
26  import org.orekit.propagation.events.FieldEventDetector;
27  
28  import java.util.stream.Stream;
29  
30  /**
31   * Class for fuel cost with Cartesian coordinates.
32   * It is the integral over time of the Euclidean norm of the thrust vector.
33   *
34   * @author Romain Serra
35   * @see CartesianCost
36   * @since 13.0
37   */
38  public class FieldCartesianFuelCost<T extends CalculusFieldElement<T>> extends FieldAbstractCartesianCost<T> {
39  
40      /** Maximum value of thrust force Euclidean norm. */
41      private final T maximumThrustMagnitude;
42  
43      /** Detection settings for singularity detection. */
44      private final FieldEventDetectionSettings<T> eventDetectionSettings;
45  
46      /**
47       * Constructor with default detection settings.
48       * @param name name
49       * @param massFlowRateFactor mass flow rate factor
50       * @param maximumThrustMagnitude maximum thrust magnitude
51       */
52      public FieldCartesianFuelCost(final String name, final T massFlowRateFactor, final T maximumThrustMagnitude) {
53          this(name, massFlowRateFactor, maximumThrustMagnitude, new FieldEventDetectionSettings<>(massFlowRateFactor.getField(),
54                  EventDetectionSettings.getDefaultEventDetectionSettings()));
55      }
56  
57      /**
58       * Constructor.
59       * @param name name
60       * @param massFlowRateFactor mass flow rate factor
61       * @param maximumThrustMagnitude maximum thrust magnitude
62       * @param eventDetectionSettings singularity event detection settings
63       */
64      public FieldCartesianFuelCost(final String name, final T massFlowRateFactor, final T maximumThrustMagnitude,
65                                    final FieldEventDetectionSettings<T> eventDetectionSettings) {
66          super(name, massFlowRateFactor);
67          this.maximumThrustMagnitude = maximumThrustMagnitude;
68          this.eventDetectionSettings = eventDetectionSettings;
69      }
70  
71      /**
72       * Getter for event detection settings.
73       * @return detection settings.
74       */
75      public FieldEventDetectionSettings<T> getEventDetectionSettings() {
76          return eventDetectionSettings;
77      }
78  
79      /** Getter for maximum thrust magnitude.
80       * @return maximum thrust
81       */
82      public T getMaximumThrustMagnitude() {
83          return maximumThrustMagnitude;
84      }
85  
86      /**
87       * Evaluate switching function (whose sign determines the bang-bang control profile).
88       * @param adjointVariables adjoint vector
89       * @param mass mass
90       * @return value of switch function
91       */
92      private T evaluateFieldSwitchFunction(final T[] adjointVariables, final T mass) {
93          T switchFunction = getFieldAdjointVelocityNorm(adjointVariables).divide(mass).subtract(1.);
94          if (getAdjointDimension() > 6) {
95              switchFunction = switchFunction.subtract(adjointVariables[6].multiply(getMassFlowRateFactor()));
96          }
97          return switchFunction;
98      }
99  
100     /**
101      * Computes the direction of thrust.
102      * @param adjointVariables adjoint vector
103      * @return thrust direction
104      */
105     private FieldVector3D<T> getFieldThrustDirection(final T[] adjointVariables) {
106         return new FieldVector3D<>(adjointVariables[3], adjointVariables[4], adjointVariables[5]).normalize();
107     }
108 
109     /** {@inheritDoc} */
110     @Override
111     public FieldVector3D<T> getFieldThrustAccelerationVector(final T[] adjointVariables, final T mass) {
112         final T switchFunction = evaluateFieldSwitchFunction(adjointVariables, mass);
113         if (switchFunction.getReal() > 0.) {
114             return getFieldThrustDirection(adjointVariables).scalarMultiply(mass.reciprocal().multiply(maximumThrustMagnitude));
115         } else {
116             return FieldVector3D.getZero(mass.getField());
117         }
118     }
119 
120     /** {@inheritDoc} */
121     @Override
122     public void updateFieldAdjointDerivatives(final T[] adjointVariables, final T mass,
123                                               final T[] adjointDerivatives) {
124         if (getAdjointDimension() > 6) {
125             final T switchFunction = evaluateFieldSwitchFunction(adjointVariables, mass);
126             if (switchFunction.getReal() > 0.) {
127                 adjointDerivatives[6] = adjointDerivatives[6].add(getFieldAdjointVelocityNorm(adjointVariables)
128                         .multiply(maximumThrustMagnitude).divide(mass.square()));
129             }
130         }
131     }
132 
133     /** {@inheritDoc} */
134     @Override
135     public T getFieldHamiltonianContribution(final T[] adjointVariables, final T mass) {
136         final FieldVector3D<T> thrustForce = getFieldThrustAccelerationVector(adjointVariables, mass).scalarMultiply(mass);
137         return thrustForce.getNorm().negate();
138     }
139 
140     /** {@inheritDoc} */
141     @Override
142     public Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
143         return Stream.of(new FieldSwitchDetector(getEventDetectionSettings()));
144     }
145 
146     @Override
147     public CartesianFuelCost toCartesianCost() {
148         return new CartesianFuelCost(getAdjointName(), getMassFlowRateFactor().getReal(), maximumThrustMagnitude.getReal(),
149                 getEventDetectionSettings().toEventDetectionSettings());
150     }
151 
152     /**
153      * Field event detector for bang-bang switches.
154      */
155     class FieldSwitchDetector extends FieldControlSwitchDetector<T> {
156 
157         FieldSwitchDetector(final FieldEventDetectionSettings<T> detectionSettings) {
158             super(detectionSettings);
159         }
160 
161         /** {@inheritDoc} */
162         @Override
163         public T g(final FieldSpacecraftState<T> state) {
164             final T[] adjoint = state.getAdditionalState(getAdjointName());
165             return evaluateFieldSwitchFunction(adjoint, state.getMass());
166         }
167 
168     }
169 }