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