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.models.earth.ionosphere.nequick;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.util.FastMath;
21  import org.orekit.bodies.FieldGeodeticPoint;
22  
23  /** Performs the computation of the coordinates along the integration path.
24   * @author Bryan Cazabonne
25   * @since 13.0
26   */
27  class FieldSegment<T extends CalculusFieldElement<T>> {
28  
29      /** Threshold for zenith segment. */
30      private static final double THRESHOLD = 1.0e-3;
31  
32      /** Supporting ray. */
33      private final FieldRay<T> ray;
34  
35      /** Integration start. */
36      private final T y;
37  
38      /** Odd points offset. */
39      private final T g;
40  
41      /** Integration step [m]. */
42      private final T deltaN;
43  
44      /** Number of points. */
45      private final int nbPoints;
46  
47      /**
48       * Constructor.
49       *
50       * @param n   number of intervals for integration (2 points per interval, hence 2n points will be generated)
51       * @param ray ray-perigee parameters
52       * @param s1  lower boundary of integration
53       * @param s2  upper boundary for integration
54       */
55      FieldSegment(final int n, final FieldRay<T> ray, final T s1, final T s2) {
56  
57          this.ray = ray;
58  
59          // Integration step (Eq. 195)
60          deltaN = s2.subtract(s1).divide(n);
61  
62          // Eq. 196
63          g = deltaN.multiply(0.5773502691896);
64  
65          // Eq. 197
66          y = s1.add(deltaN.subtract(g).multiply(0.5));
67  
68          nbPoints = 2 * n;
69  
70      }
71  
72      /** Get point along the ray.
73       * @param index point index (between O included and {@link #getNbPoints()} excluded)
74       * @return point on ray
75       * @since 13.0
76       */
77      public FieldGeodeticPoint<T> getPoint(final int index) {
78  
79          final int p = index / 2;
80          final T   s = y.add(deltaN.multiply(p)).add(g.multiply(index % 2));
81  
82          // Heights (Eq. 178)
83          final T height = FastMath.sqrt(s.multiply(s).add(ray.getRadius().multiply(ray.getRadius()))).
84                           subtract(NeQuickModel.RE);
85  
86          if (ray.getRadius().getReal() < THRESHOLD) {
87              // zenith segment
88              return new FieldGeodeticPoint<>(ray.getLatitude(),  ray.getLongitude(), height);
89          } else {
90              // Great circle parameters (Eq. 179 to 181)
91              final T tanDs = s.divide(ray.getRadius());
92              final T cosDs = FastMath.sqrt(tanDs.multiply(tanDs).add(1.0)).reciprocal();
93              final T sinDs = tanDs.multiply(cosDs);
94  
95              // Latitude (Eq. 182 to 183)
96              final T sinLatS =
97                  ray.getScLat().sin().multiply(cosDs).add(ray.getScLat().cos().multiply(sinDs).multiply(ray.getCosineAz()));
98              final T cosLatS = FastMath.sqrt(sinLatS.multiply(sinLatS).negate().add(1.0));
99  
100             // Longitude (Eq. 184 to 187)
101             final T sinLonS = sinDs.multiply(ray.getSineAz()).multiply(ray.getScLat().cos());
102             final T cosLonS = cosDs.subtract(ray.getScLat().sin().multiply(sinLatS));
103 
104             return new FieldGeodeticPoint<>(FastMath.atan2(sinLatS, cosLatS),
105                                             FastMath.atan2(sinLonS, cosLonS).add(ray.getLongitude()),
106                                             height);
107 
108         }
109     }
110 
111     /** Get number of points.
112      * <p>
113      * Note there are 2 points per interval, so {@code index} must be between 0 (included)
114      * and 2n (excluded) for a segment built with {@code n} intervals
115      * </p>
116      * @return number of points
117      */
118     public int getNbPoints() {
119         return nbPoints;
120     }
121 
122     /**
123      * Get the integration step.
124      *
125      * @return the integration step in meters
126      */
127     public T getInterval() {
128         return deltaN;
129     }
130 
131 }