1   /* Copyright 2022-2025 Thales Alenia Space
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.models.earth.ionosphere.nequick;
18  
19  import org.hipparchus.util.FastMath;
20  import org.hipparchus.util.SinCos;
21  import org.orekit.time.DateTimeComponents;
22  
23  /**
24   * Fourier time series for the NeQuick model.
25   * @see NeQuickModel#computeFourierTimeSeries(DateTimeComponents, double)
26   * @author Luc Maisonobe
27   * @since 13.0.1
28   */
29  public class FourierTimeSeries {
30  
31      /** Date. */
32      private final DateTimeComponents dateTime;
33  
34      /** Effective ionisation level. */
35      private final double az;
36  
37      /** Effective sunspot number (Eq. 19). */
38      private final double azr;
39  
40      /** Fourier time series for foF2. */
41      private final double[] cf2;
42  
43      /** Fourier time series for M(3000)F2. */
44      private final double[] cm3;
45  
46      /**
47       * Simple constructor.
48       * @param dateTime   current date time components
49       * @param az         effective ionisation level
50       * @param flattenF2  F2 coefficients used by the F2 layer (flatten array)
51       * @param flattenFm3 Fm3 coefficients used by the M(3000)F2 layer (flatten array)
52       */
53      FourierTimeSeries(final DateTimeComponents dateTime, final double az,
54                        final double[] flattenF2, final double[] flattenFm3) {
55  
56          this.dateTime = dateTime;
57          this.az       = az;
58  
59          // Effective sunspot number (Eq. 19)
60          this.azr = FastMath.sqrt(167273.0 + (az - 63.7) * 1123.6) - 408.99;
61  
62          // Hours
63          final double hours = dateTime.getTime().getSecondsInUTCDay() / 3600.0;
64  
65          // Time argument (Eq. 49)
66          final double t = FastMath.toRadians(15 * hours) - FastMath.PI;
67  
68          // Compute Fourier time series for foF2 and M(3000)F2
69          final double[] scT = sinCos(t, 6);
70          this.cf2 = computeCF2(flattenF2, scT);
71          this.cm3 = computeCm3(flattenFm3, scT);
72  
73      }
74  
75      /** Get date time components.
76       * @return date time components
77       */
78      public DateTimeComponents getDateTime() {
79          return dateTime;
80      }
81  
82      /** Get effective ionisation level.
83       * @return effective ionisation level
84       */
85      public double getAz() {
86          return az;
87      }
88  
89      /** Get effective sunspot number.
90       * @return effective sunspot number
91       */
92      public double getAzr() {
93          return azr;
94      }
95  
96      /** Get Fourier time series for foF2.
97       * <p>
98       * Beware that for efficiency purposes, this method returns
99       * a reference to an internal array; this is the reason why
100      * this method visibility is limited to package level.
101      * </p>
102      * @return Fourier time series for foF2 (reference to an internal array)
103      */
104     double[] getCf2Reference() {
105         return cf2;
106     }
107 
108     /** Get Fourier time series for M(3000)F2.
109      * <p>
110      * Beware that for efficiency purposes, this method returns
111      * a reference to an internal array; this is the reason why
112      * this method visibility is limited to package level.
113      * </p>
114      * @return Fourier time series for M(3000)F2 (reference to an internal array)
115      */
116     double[] getCm3Reference() {
117         return cm3;
118     }
119 
120     /** Computes cf2 coefficients.
121      * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
122      * @param scT       sines/cosines array of time argument
123      * @return the cf2 coefficients array
124      */
125     private double[] computeCF2(final double[] flattenF2, final double[] scT) {
126 
127         // interpolation coefficients for effective spot number
128         final double azr01   = azr * 0.01;
129         final double omazr01 = 1 - azr01;
130 
131         // Eq. 44 and Eq. 50 merged into one loop
132         final double[] array   = new double[76];
133         int            index = 0;
134         for (int i = 0; i < array.length; i++) {
135             array[i] =
136                     omazr01 * flattenF2[index] + azr01 * flattenF2[index + 1] +
137                    (omazr01 * flattenF2[index + 2] + azr01 * flattenF2[index + 3]) * scT[0] +
138                    (omazr01 * flattenF2[index + 4] + azr01 * flattenF2[index + 5]) * scT[1] +
139                    (omazr01 * flattenF2[index + 6] + azr01 * flattenF2[index + 7]) * scT[2] +
140                    (omazr01 * flattenF2[index + 8] + azr01 * flattenF2[index + 9]) * scT[3] +
141                    (omazr01 * flattenF2[index + 10] + azr01 * flattenF2[index + 11]) * scT[4] +
142                    (omazr01 * flattenF2[index + 12] + azr01 * flattenF2[index + 13]) * scT[5] +
143                    (omazr01 * flattenF2[index + 14] + azr01 * flattenF2[index + 15]) * scT[6] +
144                    (omazr01 * flattenF2[index + 16] + azr01 * flattenF2[index + 17]) * scT[7] +
145                    (omazr01 * flattenF2[index + 18] + azr01 * flattenF2[index + 19]) * scT[8] +
146                    (omazr01 * flattenF2[index + 20] + azr01 * flattenF2[index + 21]) * scT[9] +
147                    (omazr01 * flattenF2[index + 22] + azr01 * flattenF2[index + 23]) * scT[10] +
148                    (omazr01 * flattenF2[index + 24] + azr01 * flattenF2[index + 25]) * scT[11];
149             index += 26;
150         }
151         return array;
152     }
153 
154     /** Computes Cm3 coefficients.
155      * @param flattenFm3 Fm3 coefficients used by the M(3000)F2 layer (flatten array)
156      * @param scT        sines/cosines array of time argument
157      * @return the Cm3 coefficients array
158      */
159     private double[] computeCm3(final double[] flattenFm3, final double[] scT) {
160 
161         // interpolation coefficients for effective spot number
162         final double azr01   = azr * 0.01;
163         final double omazr01 = 1 - azr01;
164 
165         // Eq. 44 and Eq. 51 merged into one loop
166         final double[] array   = new double[49];
167         int            index = 0;
168         for (int i = 0; i < array.length; i++) {
169             array[i] =
170                     omazr01 * flattenFm3[index] + azr01 * flattenFm3[index + 1] +
171                    (omazr01 * flattenFm3[index + 2] + azr01 * flattenFm3[index + 3]) * scT[0] +
172                    (omazr01 * flattenFm3[index + 4] + azr01 * flattenFm3[index + 5]) * scT[1] +
173                    (omazr01 * flattenFm3[index + 6] + azr01 * flattenFm3[index + 7]) * scT[2] +
174                    (omazr01 * flattenFm3[index + 8] + azr01 * flattenFm3[index + 9]) * scT[3] +
175                    (omazr01 * flattenFm3[index + 10] + azr01 * flattenFm3[index + 11]) * scT[4] +
176                    (omazr01 * flattenFm3[index + 12] + azr01 * flattenFm3[index + 13]) * scT[5] +
177                    (omazr01 * flattenFm3[index + 14] + azr01 * flattenFm3[index + 15]) * scT[6] +
178                    (omazr01 * flattenFm3[index + 16] + azr01 * flattenFm3[index + 17]) * scT[7];
179             index += 18;
180         }
181         return array;
182     }
183 
184     /** Compute sines and cosines.
185      * @param a argument
186      * @param n number of terms
187      * @return sin(a), cos(a), sin(2a), cos(2a) … sin(n a), cos(n a) array
188      */
189     static double[] sinCos(final double a, final int n) {
190 
191         final SinCos sc0 = FastMath.sinCos(a);
192         SinCos sci = sc0;
193         final double[] sc = new double[2 * n];
194         int isc = 0;
195         sc[isc++] = sci.sin();
196         sc[isc++] = sci.cos();
197         for (int i = 1; i < n; i++) {
198             sci = SinCos.sum(sc0, sci);
199             sc[isc++] = sci.sin();
200             sc[isc++] = sci.cos();
201         }
202 
203         return sc;
204 
205     }
206 
207 }