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.geometry.euclidean.threed.Vector3D;
21  import org.orekit.propagation.SpacecraftState;
22  import org.orekit.propagation.events.EventDetectionSettings;
23  
24  /**
25   * Abstract class for energy cost with Cartesian coordinates.
26   * An energy cost is proportional to the integral over time of the squared Euclidean norm of the control vector, often scaled with 1/2.
27   * 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.
28   *
29   * @author Romain Serra
30   * @since 12.2
31   */
32  abstract class CartesianEnergyConsideringMass extends AbstractCartesianCost {
33  
34      /** Detection settings for singularity detection. */
35      private final EventDetectionSettings eventDetectionSettings;
36  
37      /**
38       * Constructor.
39       * @param name name
40       * @param massFlowRateFactor mass flow rate factor
41       * @param eventDetectionSettings settings for singularity detections
42       */
43      protected CartesianEnergyConsideringMass(final String name, final double massFlowRateFactor,
44                                               final EventDetectionSettings eventDetectionSettings) {
45          super(name, massFlowRateFactor);
46          this.eventDetectionSettings = eventDetectionSettings;
47      }
48  
49      /**
50       * Getter for event detection settings.
51       * @return detection settings.
52       */
53      public EventDetectionSettings getEventDetectionSettings() {
54          return eventDetectionSettings;
55      }
56  
57      /** {@inheritDoc} */
58      @Override
59      public Vector3D getThrustAccelerationVector(final double[] adjointVariables, final double mass) {
60          return getThrustDirection(adjointVariables).scalarMultiply(getThrustForceNorm(adjointVariables, mass) / mass);
61      }
62  
63      /**
64       * Computes the direction of thrust.
65       * @param adjointVariables adjoint vector
66       * @return thrust direction
67       */
68      protected Vector3D getThrustDirection(final double[] adjointVariables) {
69          return new Vector3D(adjointVariables[3], adjointVariables[4], adjointVariables[5]).normalize();
70      }
71  
72      /**
73       * Computes the Euclidean norm of the thrust force.
74       * @param adjointVariables adjoint vector
75       * @param mass mass
76       * @return norm of thrust
77       */
78      protected abstract double getThrustForceNorm(double[] adjointVariables, double mass);
79  
80      /** {@inheritDoc} */
81      @Override
82      public void updateAdjointDerivatives(final double[] adjointVariables, final double mass,
83                                           final double[] adjointDerivatives) {
84          if (getAdjointDimension() > 6) {
85              adjointDerivatives[6] += getThrustForceNorm(adjointVariables, mass) * getAdjointVelocityNorm(adjointVariables) / (mass * mass);
86          }
87      }
88  
89      /** {@inheritDoc} */
90      @Override
91      public double getHamiltonianContribution(final double[] adjointVariables, final double mass) {
92          final Vector3D thrustForce = getThrustAccelerationVector(adjointVariables, mass).scalarMultiply(mass);
93          return -thrustForce.getNormSq() / 2.;
94      }
95  
96      /**
97       * Event detector for singularities in adjoint dynamics.
98       */
99      class SingularityDetector extends ControlSwitchDetector {
100 
101         /** Value to detect. */
102         private final double detectionValue;
103 
104         /**
105          * Constructor.
106          * @param detectionSettings detection settings
107          * @param detectionValue value to detect
108          */
109         SingularityDetector(final EventDetectionSettings detectionSettings, final double detectionValue) {
110             super(detectionSettings);
111             this.detectionValue = detectionValue;
112         }
113 
114         /** {@inheritDoc} */
115         @Override
116         public double g(final SpacecraftState state) {
117             final double[] adjoint = state.getAdditionalState(getAdjointName());
118             return evaluateVariablePart(adjoint, state.getMass()) - detectionValue;
119         }
120 
121         /**
122          * Evaluate variable part of singularity function.
123          * @param adjointVariables adjoint vector
124          * @param mass mass
125          * @return singularity function without the constant part
126          */
127         private double evaluateVariablePart(final double[] adjointVariables, final double mass) {
128             final double adjointVelocityNorm = getAdjointVelocityNorm(adjointVariables);
129             double variablePart = adjointVelocityNorm / mass;
130             if (getAdjointDimension() > 6) {
131                 variablePart -= getMassFlowRateFactor() * adjointVariables[6];
132             }
133             return variablePart;
134         }
135 
136     }
137 
138 }