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.geometry.euclidean.threed.FieldVector3D;
22  import org.orekit.propagation.FieldSpacecraftState;
23  import org.orekit.propagation.events.FieldEventDetectionSettings;
24  
25  /**
26   * Abstract class for energy cost with Cartesian coordinates and non-zero mass flow rate.
27   * An energy cost is proportional to the integral over time of the squared Euclidean norm of the control vector, often scaled with 1/2.
28   * This type of cost is not optimal in terms of mass consumption, however its solutions showcase a smoother behavior favorable for convergence in shooting techniques.
29   *
30   * @param <T> field type
31   * @author Romain Serra
32   * @see FieldCartesianCost
33   * @see CartesianEnergyConsideringMass
34   * @since 13.0
35   */
36  abstract class FieldCartesianEnergyConsideringMass<T extends CalculusFieldElement<T>> extends FieldAbstractCartesianCost<T> {
37  
38      /** Detection settings for singularity detection. */
39      private final FieldEventDetectionSettings<T> eventDetectionSettings;
40  
41      /**
42       * Constructor.
43       * @param name name
44       * @param massFlowRateFactor mass flow rate factor
45       * @param eventDetectionSettings settings for singularity detections
46       */
47      protected FieldCartesianEnergyConsideringMass(final String name, final T massFlowRateFactor,
48                                                    final FieldEventDetectionSettings<T> eventDetectionSettings) {
49          super(name, massFlowRateFactor);
50          this.eventDetectionSettings = eventDetectionSettings;
51      }
52  
53      /**
54       * Getter for event detection settings.
55       * @return detection settings.
56       */
57      public FieldEventDetectionSettings<T> getEventDetectionSettings() {
58          return eventDetectionSettings;
59      }
60  
61      /** {@inheritDoc} */
62      @Override
63      public FieldVector3D<T> getFieldThrustAccelerationVector(final T[] adjointVariables, final T mass) {
64          return getFieldThrustDirection(adjointVariables).scalarMultiply(getFieldThrustForceNorm(adjointVariables, mass).divide(mass));
65      }
66  
67      /**
68       * Computes the direction of thrust.
69       * @param adjointVariables adjoint vector
70       * @return thrust direction
71       */
72      protected FieldVector3D<T> getFieldThrustDirection(final T[] adjointVariables) {
73          return new FieldVector3D<>(adjointVariables[3], adjointVariables[4], adjointVariables[5]).normalize();
74      }
75  
76      /**
77       * Computes the Euclidean norm of the thrust force.
78       * @param adjointVariables adjoint vector
79       * @param mass mass
80       * @return norm of thrust
81       */
82      protected abstract T getFieldThrustForceNorm(T[] adjointVariables, T mass);
83  
84      /** {@inheritDoc} */
85      @Override
86      public void updateFieldAdjointDerivatives(final T[] adjointVariables, final T mass, final T[] adjointDerivatives) {
87          if (getAdjointDimension() > 6) {
88              adjointDerivatives[6] = adjointDerivatives[6].add(getFieldThrustForceNorm(adjointVariables, mass)
89                      .multiply(getFieldAdjointVelocityNorm(adjointVariables)).divide(mass.square()));
90          }
91      }
92  
93      /** {@inheritDoc} */
94      @Override
95      public T getFieldHamiltonianContribution(final T[] adjointVariables, final T mass) {
96          final FieldVector3D<T> thrustForce = getFieldThrustAccelerationVector(adjointVariables, mass).scalarMultiply(mass);
97          return thrustForce.getNormSq().multiply(-1. / 2.);
98      }
99  
100     /**
101      * Field event detector for singularities in adjoint dynamics.
102      */
103     class FieldSingularityDetector extends FieldControlSwitchDetector<T> {
104 
105         /** Value to detect. */
106         private final T detectionValue;
107 
108         /**
109          * Constructor.
110          * @param detectionSettings detection settings
111          * @param detectionValue value to detect
112          */
113         FieldSingularityDetector(final FieldEventDetectionSettings<T> detectionSettings, final T detectionValue) {
114             super(detectionSettings);
115             this.detectionValue = detectionValue;
116         }
117 
118         /** {@inheritDoc} */
119         @Override
120         public T g(final FieldSpacecraftState<T> state) {
121             final T[] adjoint = state.getAdditionalState(getAdjointName());
122             return evaluateVariablePart(adjoint, state.getMass()).subtract(detectionValue);
123         }
124 
125         /**
126          * Evaluate variable part of singularity function.
127          * @param adjointVariables adjoint vector
128          * @param mass mass
129          * @return singularity function without the constant part
130          */
131         private T evaluateVariablePart(final T[] adjointVariables, final T mass) {
132             final T adjointVelocityNorm = getFieldAdjointVelocityNorm(adjointVariables);
133             T variablePart = adjointVelocityNorm.divide(mass);
134             if (getAdjointDimension() > 6) {
135                 variablePart = variablePart.subtract(adjointVariables[6].multiply(getMassFlowRateFactor()));
136             }
137             return variablePart;
138         }
139 
140     }
141 }