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.gravity.potential;
18  
19  import org.hipparchus.util.FastMath;
20  import org.hipparchus.util.SinCos;
21  import org.orekit.data.BodiesElements;
22  
23  /** Container for ocen tides coefficients for one tide wave.
24   * @see org.orekit.forces.gravity.OceanTides
25   * @author Luc Maisonobe
26   * @since 6.1
27   * @see OceanTidesReader
28   */
29  public class OceanTidesWave {
30  
31      /** Waves of degree 0 and 1 do not affect spacecrafts. */
32      private static final int START_DEGREE = 2;
33  
34      /** Maximum supported degree. */
35      private final int degree;
36  
37      /** Maximum supported order. */
38      private final int order;
39  
40      /** Doodson number for the wave. */
41      private final int doodson;
42  
43      /** Coefficient for γ = GMST + π tide parameter. */
44      private final int cGamma;
45  
46      /** Coefficient for mean anomaly of the Moon. */
47      private final int cL;
48  
49      /** Coefficient for mean anomaly of the Sun. */
50      private final int cLPrime;
51  
52      /** Coefficient for L - Ω where L is the mean longitude of the Moon. */
53      private final int cF;
54  
55      /** Coefficient for mean elongation of the Moon from the Sun. */
56      private final int cD;
57  
58      /** Coefficient for mean longitude of the ascending node of the Moon. */
59      private final int cOmega;
60  
61      /** C<sub>n,m</sub><sup>+</sup> coefficients. */
62      private final double[][] cPlus;
63  
64      /** S<sub>n,m</sub><sup>+</sup> coefficients. */
65      private final double[][] sPlus;
66  
67      /** C<sub>n,m</sub><sup>-</sup> coefficients. */
68      private final double[][] cMinus;
69  
70      /** S<sub>n,m</sub><sup>-</sup> coefficients. */
71      private final double[][] sMinus;
72  
73      /** Simple constructor.
74       * @param doodson Doodson number for the wave
75       * @param degree max degree present in the coefficients array
76       * @param order max order present in the coefficients array
77       * @param coefficients C<sub>n,m</sub><sup>+</sup>, S<sub>n,m</sub><sup>+</sup>,
78       * C<sub>n,m</sub><sup>-</sup> and S<sub>n,m</sub><sup>-</sup> coefficients
79       */
80      public OceanTidesWave(final int doodson, final int degree, final int order,
81                            final double[][][] coefficients) {
82  
83          this.doodson = doodson;
84  
85          // compute Doodson arguments from Doodson number
86          final int cPs     = ( doodson           % 10) - 5;
87          final int cNPrime = ((doodson / 10)     % 10) - 5;
88          final int cP      = ((doodson / 100)    % 10) - 5;
89          final int cH      = ((doodson / 1000)   % 10) - 5;
90          final int cS      = ((doodson / 10000)  % 10) - 5;
91          final int cTau    =  (doodson / 100000) % 10;
92  
93          // compute Delaunay arguments from Doodson arguments
94          this.cGamma  =  cTau;
95          this.cL      = -cP;
96          this.cLPrime = -cPs;
97          this.cF      = -cTau + cS + cH + cP + cPs;
98          this.cD      = -cH - cPs;
99          this.cOmega  = -cTau + cS + cH + cP - cNPrime + cPs;
100 
101         this.degree   = degree;
102         this.order    = order;
103 
104         // distribute the coefficients
105         final int rows = degree + 1;
106         this.cPlus     = new double[rows][];
107         this.sPlus     = new double[rows][];
108         this.cMinus    = new double[rows][];
109         this.sMinus    = new double[rows][];
110         for (int i = 0; i <= degree; ++i) {
111             final int m = FastMath.min(i, order) + 1;
112             final double[][] row = coefficients[i];
113             cPlus[i]  = new double[m];
114             sPlus[i]  = new double[m];
115             cMinus[i] = new double[m];
116             sMinus[i] = new double[m];
117             for (int j = 0; j < m; ++j) {
118                 cPlus[i][j]  = row[j][0];
119                 sPlus[i][j]  = row[j][1];
120                 cMinus[i][j] = row[j][2];
121                 sMinus[i][j] = row[j][3];
122             }
123         }
124 
125     }
126 
127     /** Get the maximum supported degree.
128      * @return maximum supported degree
129      */
130     public int getMaxDegree() {
131         return degree;
132     }
133 
134     /** Get the maximum supported order.
135      * @return maximum supported order
136      */
137     public int getMaxOrder() {
138         return order;
139     }
140 
141     /** Get the Doodson number for the wave.
142      * @return Doodson number for the wave
143      */
144     public int getDoodson() {
145         return doodson;
146     }
147 
148     /** Add the contribution of the wave to Stokes coefficients.
149      * @param elements nutation elements
150      * @param cnm spherical harmonic cosine coefficients table to add contribution too
151      * @param snm spherical harmonic sine coefficients table to add contribution too
152      */
153     public void addContribution(final BodiesElements elements,
154                                 final double[][] cnm, final double[][] snm) {
155 
156         final double thetaF = cGamma * elements.getGamma() +
157                               cL * elements.getL() + cLPrime * elements.getLPrime() + cF * elements.getF() +
158                               cD * elements.getD() + cOmega * elements.getOmega();
159         final SinCos sc = FastMath.sinCos(thetaF);
160 
161         for (int i = START_DEGREE; i <= degree; ++i) {
162             for (int j = 0; j <= FastMath.min(i, order); ++j) {
163                 // from IERS conventions 2010, section 6.3, equation 6.15
164                 cnm[i][j] += (cPlus[i][j] + cMinus[i][j]) * sc.cos() + (sPlus[i][j] + sMinus[i][j]) * sc.sin();
165                 snm[i][j] += (sPlus[i][j] - sMinus[i][j]) * sc.cos() - (cPlus[i][j] - cMinus[i][j]) * sc.sin();
166             }
167         }
168 
169     }
170 
171 }