1   /* Copyright 2002-2025 CS GROUP
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.forces.empirical;
18  
19  import java.util.ArrayList;
20  import java.util.Collections;
21  import java.util.List;
22  
23  import org.hipparchus.CalculusFieldElement;
24  import org.hipparchus.util.FastMath;
25  import org.hipparchus.util.MathUtils;
26  import org.orekit.propagation.FieldSpacecraftState;
27  import org.orekit.propagation.SpacecraftState;
28  import org.orekit.time.AbsoluteDate;
29  import org.orekit.utils.ParameterDriver;
30  
31  /** Harmonic acceleration model.
32   * @since 10.3
33   * @author Luc Maisonobe
34   * @author Bryan Cazabonne
35   */
36  public class HarmonicAccelerationModel implements AccelerationModel {
37  
38      /** Amplitude scaling factor.
39       * <p>
40       * 2⁻²⁰ is the order of magnitude of third body perturbing acceleration.
41       * </p>
42       * <p>
43       * We use a power of 2 to avoid numeric noise introduction
44       * in the multiplications/divisions sequences.
45       * </p>
46       */
47      private static final double AMPLITUDE_SCALE = FastMath.scalb(1.0, -20);
48  
49      /** Phase scaling factor.
50       * <p>
51       * 2⁻²³ is the order of magnitude of an angle corresponding to one meter along
52       * track for a Low Earth Orbiting satellite.
53       * </p>
54       * <p>
55       * We use a power of 2 to avoid numeric noise introduction
56       * in the multiplications/divisions sequences.
57       * </p>
58       */
59      private static final double PHASE_SCALE = FastMath.scalb(1.0, -23);
60  
61      /** Drivers for the parameters. */
62      private final List<ParameterDriver> drivers;
63  
64      /** Reference date for computing phase. */
65      private AbsoluteDate referenceDate;
66  
67      /** Angular frequency ω = 2kπ/T. */
68      private final double omega;
69  
70      /** Simple constructor.
71       * @param prefix prefix to use for parameter drivers
72       * @param referenceDate reference date for computing polynomials, if null
73       * the reference date will be automatically set at propagation start
74       * @param fundamentalPeriod fundamental period (typically set to initial orbit
75       * {@link org.orekit.orbits.Orbit#getKeplerianPeriod() Keplerian period})
76       * @param harmonicMultiplier multiplier to compute harmonic period from
77       * fundamental period)
78       */
79      public HarmonicAccelerationModel(final String prefix, final AbsoluteDate referenceDate,
80                                       final double fundamentalPeriod, final int harmonicMultiplier) {
81          this.referenceDate = referenceDate;
82          this.omega         = harmonicMultiplier * MathUtils.TWO_PI / fundamentalPeriod;
83          this.drivers       = new ArrayList<>(2);
84          drivers.add(new ParameterDriver(prefix + " γ",
85                                          0.0, AMPLITUDE_SCALE, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY));
86          drivers.add(new ParameterDriver(prefix + " φ",
87                                          0.0, PHASE_SCALE, -MathUtils.TWO_PI, MathUtils.TWO_PI));
88      }
89  
90      /** {@inheritDoc} */
91      @Override
92      public void init(final SpacecraftState initialState, final AbsoluteDate target) {
93          if (referenceDate == null) {
94              referenceDate = initialState.getDate();
95          }
96      }
97  
98      /** {@inheritDoc} */
99      @Override
100     public double signedAmplitude(final SpacecraftState state,
101                                   final double[] parameters) {
102         final double dt = state.getDate().durationFrom(referenceDate);
103         return parameters[0] * FastMath.sin(dt * omega + parameters[1]);
104     }
105 
106     /** {@inheritDoc} */
107     @Override
108     public <T extends CalculusFieldElement<T>> T signedAmplitude(final FieldSpacecraftState<T> state,
109                                                              final T[] parameters) {
110         final T dt = state.getDate().durationFrom(referenceDate);
111         return parameters[0].multiply(dt.multiply(omega).add(parameters[1]).sin());
112     }
113 
114     /** {@inheritDoc} */
115     @Override
116     public List<ParameterDriver> getParametersDrivers() {
117         return Collections.unmodifiableList(drivers);
118     }
119 
120 }