1 /* Copyright 2002-2023 Luc Maisonobe
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.utils;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21 import org.hipparchus.geometry.euclidean.threed.Vector3D;
22 import org.hipparchus.util.FastMath;
23 import org.orekit.bodies.OneAxisEllipsoid;
24 import org.orekit.propagation.FieldSpacecraftState;
25 import org.orekit.propagation.SpacecraftState;
26
27 /** Computation engine for occultation events.
28 * @author Luc Maisonobe
29 * @since 12.0
30 */
31 public class OccultationEngine {
32
33 /** Occulting body. */
34 private final OneAxisEllipsoid occulting;
35
36 /** Occulted body. */
37 private final ExtendedPVCoordinatesProvider occulted;
38
39 /** Occulted body radius (m). */
40 private final double occultedRadius;
41
42 /** Build a new occultation engine.
43 * @param occulted the body to be occulted
44 * @param occultedRadius the radius of the body to be occulted (m)
45 * @param occulting the occulting body
46 */
47 public OccultationEngine(final ExtendedPVCoordinatesProvider occulted, final double occultedRadius,
48 final OneAxisEllipsoid occulting) {
49 this.occulted = occulted;
50 this.occultedRadius = FastMath.abs(occultedRadius);
51 this.occulting = occulting;
52 }
53
54 /** Getter for the occulting body.
55 * @return the occulting body
56 */
57 public OneAxisEllipsoid getOcculting() {
58 return occulting;
59 }
60
61 /** Getter for the occulted body.
62 * @return the occulted body
63 */
64 public ExtendedPVCoordinatesProvider getOcculted() {
65 return occulted;
66 }
67
68 /** Getter for the occultedRadius.
69 * @return the occultedRadius
70 */
71 public double getOccultedRadius() {
72 return occultedRadius;
73 }
74
75 /** Compute the occultation angles as seen from a spacecraft.
76 * @param state the current state information: date, kinematics, attitude
77 * @return occultation angles
78 */
79 public OccultationAngles angles(final SpacecraftState state) {
80
81 final Vector3D psat = state.getPosition(occulting.getBodyFrame());
82 final Vector3D pted = occulted.getPosition(state.getDate(), occulting.getBodyFrame());
83 final Vector3D plimb = occulting.pointOnLimb(psat, pted);
84 final Vector3D ps = psat.subtract(pted);
85 final Vector3D pi = psat.subtract(plimb);
86 final double angle = Vector3D.angle(ps, psat);
87 final double rs = FastMath.asin(occultedRadius / ps.getNorm());
88 final double ro = Vector3D.angle(pi, psat);
89 if (Double.isNaN(rs)) {
90 // we are inside the occulted body…
91 // set up dummy values consistent with full lighting (assuming occulted is the Sun)
92 return new OccultationAngles(FastMath.PI, 0.0, 0.0);
93 } else {
94 // regular case, we can compute limit angles as seen from spacecraft
95 return new OccultationAngles(angle, ro, rs);
96 }
97
98 }
99
100 /** Compute the occultation angles as seen from a spacecraft.
101 * @param state the current state information: date, kinematics, attitude
102 * @param <T> the type of the field elements
103 * @return occultation angles
104 */
105 public <T extends CalculusFieldElement<T>> FieldOccultationAngles<T> angles(final FieldSpacecraftState<T> state) {
106
107 final FieldVector3D<T> psat = state.getPosition(occulting.getBodyFrame());
108 final FieldVector3D<T> pted = occulted.getPosition(state.getDate(), occulting.getBodyFrame());
109 final FieldVector3D<T> plimb = occulting.pointOnLimb(psat, pted);
110 final FieldVector3D<T> ps = psat.subtract(pted);
111 final FieldVector3D<T> pi = psat.subtract(plimb);
112 final T angle = FieldVector3D.angle(ps, psat);
113 final T rs = FastMath.asin(ps.getNorm().reciprocal().multiply(occultedRadius));
114 final T ro = FieldVector3D.angle(pi, psat);
115 if (rs.isNaN()) {
116 // we are inside the occulted body…
117 // set up dummy values consistent with full lighting (assuming occulted is the Sun)
118 final T zero = rs.getField().getZero();
119 return new FieldOccultationAngles<>(zero.newInstance(FastMath.PI), zero, zero);
120 } else {
121 // regular case, we can compute limit angles as seen from spacecraft
122 return new FieldOccultationAngles<>(angle, ro, rs);
123 }
124
125 }
126
127 /** Container for occultation angles.
128 * @since 12.0
129 */
130 public static class OccultationAngles {
131
132 /** Apparent separation between occulting and occulted directions. */
133 private final double separation;
134
135 /** Limb radius in occulting/occulted plane. */
136 private final double limbRadius;
137
138 /** Apparent radius of occulted body. */
139 private final double occultedApparentRadius;
140
141 /** Simple constructor.
142 * @param separation apparent separation between occulting and occulted directions (rad)
143 * @param limbRadius limb radius in occulting/occulted plane (rad)
144 * @param occultedApparentRadius apparent radius of occulted body (rad)
145 */
146 OccultationAngles(final double separation, final double limbRadius, final double occultedApparentRadius) {
147 this.separation = separation;
148 this.limbRadius = limbRadius;
149 this.occultedApparentRadius = occultedApparentRadius;
150 }
151
152 /** Get apparent separation between occulting and occulted directions.
153 * @return apparent separation between occulting and occulted directions (rad)
154 */
155 public double getSeparation() {
156 return separation;
157 }
158
159 /** Get limb radius in occulting/occulted plane.
160 * @return limb radius in occulting/occulted plane (rad)
161 */
162 public double getLimbRadius() {
163 return limbRadius;
164 }
165
166 /** Get apparent radius of occulted body.
167 * @return apparent radius of occulted body (rad)
168 */
169 public double getOccultedApparentRadius() {
170 return occultedApparentRadius;
171 }
172
173 }
174
175 /** Container for occultation angles.
176 * @param <T> the type of the field elements
177 * @since 12.0
178 */
179 public static class FieldOccultationAngles<T extends CalculusFieldElement<T>> {
180
181 /** Apparent separation between occulting and occulted directions. */
182 private final T separation;
183
184 /** Limb radius in occulting/occulted plane. */
185 private final T limbRadius;
186
187 /** Apparent radius of occulted body. */
188 private final T occultedApparentRadius;
189
190 /** Simple constructor.
191 * @param separation apparent separation between occulting and occulted directions (rad)
192 * @param limbRadius limb radius in occulting/occulted plane (rad)
193 * @param occultedApparentRadius apparent radius of occulted body (rad)
194 */
195 FieldOccultationAngles(final T separation, final T limbRadius, final T occultedApparentRadius) {
196 this.separation = separation;
197 this.limbRadius = limbRadius;
198 this.occultedApparentRadius = occultedApparentRadius;
199 }
200
201 /** Get apparent separation between occulting and occulted directions.
202 * @return apparent separation between occulting and occulted directions (rad)
203 */
204 public T getSeparation() {
205 return separation;
206 }
207
208 /** Get limb radius in occulting/occulted plane.
209 * @return limb radius in occulting/occulted plane (rad)
210 */
211 public T getLimbRadius() {
212 return limbRadius;
213 }
214
215 /** Get apparent radius of occulted body.
216 * @return apparent radius of occulted body (rad)
217 */
218 public T getOccultedApparentRadius() {
219 return occultedApparentRadius;
220 }
221
222 }
223
224 }