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