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.propagation.analytical;
18
19 import org.hipparchus.util.FastMath;
20 import org.hipparchus.util.SinCos;
21 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
22 import org.orekit.orbits.EquinoctialOrbit;
23 import org.orekit.orbits.Orbit;
24 import org.orekit.orbits.OrbitType;
25 import org.orekit.orbits.PositionAngleType;
26 import org.orekit.propagation.SpacecraftState;
27 import org.orekit.time.AbsoluteDate;
28
29 /** Analytical model for J2 effect.
30 * <p>
31 * This class computes the differential effect of J2 due to an initial orbit
32 * offset. A typical case is when an inclination maneuver changes an orbit
33 * inclination at time t₀. As ascending node drift rate depends on
34 * inclination, the change induces a time-dependent change in ascending node
35 * for later dates.
36 * </p>
37 * @see org.orekit.forces.maneuvers.SmallManeuverAnalyticalModel
38 * @author Luc Maisonobe
39 */
40 public class J2DifferentialEffect
41 implements AdapterPropagator.DifferentialEffect {
42
43 /** Reference date. */
44 private final AbsoluteDate referenceDate;
45
46 /** Differential drift on perigee argument. */
47 private final double dPaDot;
48
49 /** Differential drift on ascending node. */
50 private final double dRaanDot;
51
52 /** Indicator for applying effect before reference date. */
53 private final boolean applyBefore;
54
55 /** Simple constructor.
56 * <p>
57 * The {@code applyBefore} parameter is mainly used when the differential
58 * effect is associated with a maneuver. In this case, the parameter must be
59 * set to {@code false}.
60 * </p>
61 * @param original original state at reference date
62 * @param directEffect direct effect changing the orbit
63 * @param applyBefore if true, effect is applied both before and after
64 * reference date, if false it is only applied after reference date
65 * @param gravityField gravity field to use
66 */
67 public J2DifferentialEffect(final SpacecraftState original,
68 final AdapterPropagator.DifferentialEffect directEffect,
69 final boolean applyBefore,
70 final UnnormalizedSphericalHarmonicsProvider gravityField) {
71 this(original, directEffect, applyBefore,
72 gravityField.getAe(), gravityField.getMu(),
73 -gravityField.onDate(original.getDate()).getUnnormalizedCnm(2, 0));
74 }
75
76 /** Simple constructor.
77 * <p>
78 * The {@code applyBefore} parameter is mainly used when the differential
79 * effect is associated with a maneuver. In this case, the parameter must be
80 * set to {@code false}.
81 * </p>
82 * @param orbit0 original orbit at reference date
83 * @param orbit1 shifted orbit at reference date
84 * @param applyBefore if true, effect is applied both before and after
85 * reference date, if false it is only applied after reference date
86 * @param gravityField gravity field to use
87 */
88 public J2DifferentialEffect(final Orbit orbit0, final Orbit orbit1, final boolean applyBefore,
89 final UnnormalizedSphericalHarmonicsProvider gravityField) {
90 this(orbit0, orbit1, applyBefore,
91 gravityField.getAe(), gravityField.getMu(),
92 -gravityField.onDate(orbit0.getDate()).getUnnormalizedCnm(2, 0));
93 }
94
95 /** Simple constructor.
96 * <p>
97 * The {@code applyBefore} parameter is mainly used when the differential
98 * effect is associated with a maneuver. In this case, the parameter must be
99 * set to {@code false}.
100 * </p>
101 * @param original original state at reference date
102 * @param directEffect direct effect changing the orbit
103 * @param applyBefore if true, effect is applied both before and after
104 * reference date, if false it is only applied after reference date
105 * @param referenceRadius reference radius of the Earth for the potential model (m)
106 * @param mu central attraction coefficient (m³/s²)
107 * @param j2 un-normalized zonal coefficient (about +1.08e-3 for Earth)
108 */
109 public J2DifferentialEffect(final SpacecraftState original,
110 final AdapterPropagator.DifferentialEffect directEffect,
111 final boolean applyBefore,
112 final double referenceRadius, final double mu, final double j2) {
113 this(original.getOrbit(),
114 directEffect.apply(original.shiftedBy(0.001)).getOrbit().shiftedBy(-0.001),
115 applyBefore, referenceRadius, mu, j2);
116 }
117
118 /** Simple constructor.
119 * <p>
120 * The {@code applyBefore} parameter is mainly used when the differential
121 * effect is associated with a maneuver. In this case, the parameter must be
122 * set to {@code false}.
123 * </p>
124 * @param orbit0 original orbit at reference date
125 * @param orbit1 shifted orbit at reference date
126 * @param applyBefore if true, effect is applied both before and after
127 * reference date, if false it is only applied after reference date
128 * @param referenceRadius reference radius of the Earth for the potential model (m)
129 * @param mu central attraction coefficient (m³/s²)
130 * @param j2 un-normalized zonal coefficient (about +1.08e-3 for Earth)
131 */
132 public J2DifferentialEffect(final Orbit orbit0, final Orbit orbit1, final boolean applyBefore,
133 final double referenceRadius, final double mu, final double j2) {
134
135 this.referenceDate = orbit0.getDate();
136 this.applyBefore = applyBefore;
137
138 // extract useful parameters
139 final double a0 = orbit0.getA();
140 final double e0 = orbit0.getE();
141 final double i0 = orbit0.getI();
142 final double a1 = orbit1.getA();
143 final double e1 = orbit1.getE();
144 final double i1 = orbit1.getI();
145
146 // compute reference drifts
147 final double oMe2 = 1 - e0 * e0;
148 final double ratio = referenceRadius / (a0 * oMe2);
149 final SinCos scI = FastMath.sinCos(i0);
150 final double n = FastMath.sqrt(mu / a0) / a0;
151 final double c = ratio * ratio * n * j2;
152 final double refPaDot = 0.75 * c * (4 - 5 * scI.sin() * scI.sin());
153 final double refRaanDot = -1.5 * c * scI.cos();
154
155 // differential model on perigee argument drift
156 final double dPaDotDa = -3.5 * refPaDot / a0;
157 final double dPaDotDe = 4 * refPaDot * e0 / oMe2;
158 final double dPaDotDi = -7.5 * c * scI.sin() * scI.cos();
159 dPaDot = dPaDotDa * (a1 - a0) + dPaDotDe * (e1 - e0) + dPaDotDi * (i1 - i0);
160
161 // differential model on ascending node drift
162 final double dRaanDotDa = -3.5 * refRaanDot / a0;
163 final double dRaanDotDe = 4 * refRaanDot * e0 / oMe2;
164 final double dRaanDotDi = -refRaanDot * FastMath.tan(i0);
165 dRaanDot = dRaanDotDa * (a1 - a0) + dRaanDotDe * (e1 - e0) + dRaanDotDi * (i1 - i0);
166
167 }
168
169 /** Compute the effect of the maneuver on an orbit.
170 * @param orbit1 original orbit at t₁, without maneuver
171 * @return orbit at t₁, taking the maneuver
172 * into account if t₁ > t₀
173 * @see #apply(SpacecraftState)
174 */
175 public Orbit apply(final Orbit orbit1) {
176
177 if (orbit1.getDate().compareTo(referenceDate) <= 0 && !applyBefore) {
178 // the orbit change has not occurred yet, don't change anything
179 return orbit1;
180 }
181
182 return updateOrbit(orbit1);
183
184 }
185
186 /** {@inheritDoc} */
187 public SpacecraftState apply(final SpacecraftState state1) {
188
189 if (state1.getDate().compareTo(referenceDate) <= 0 && !applyBefore) {
190 // the orbit change has not occurred yet, don't change anything
191 return state1;
192 }
193
194 return new SpacecraftState(updateOrbit(state1.getOrbit()), state1.getAttitude()).withMass(state1.getMass());
195
196 }
197
198 /** Compute the differential effect of J2 on an orbit.
199 * @param orbit1 original orbit at t₁, without differential J2
200 * @return orbit at t₁, always taking the effect into account
201 */
202 private Orbit updateOrbit(final Orbit orbit1) {
203
204 // convert current orbital state to equinoctial elements
205 final EquinoctialOrbit original =
206 (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(orbit1);
207
208 // compute differential effect
209 final AbsoluteDate date = original.getDate();
210 final double dt = date.durationFrom(referenceDate);
211 final double dPaRaan = (dPaDot + dRaanDot) * dt;
212 final SinCos scPaRaan = FastMath.sinCos(dPaRaan);
213 final double dRaan = dRaanDot * dt;
214 final SinCos scRaan = FastMath.sinCos(dRaan);
215
216 final double ex = original.getEquinoctialEx() * scPaRaan.cos() -
217 original.getEquinoctialEy() * scPaRaan.sin();
218 final double ey = original.getEquinoctialEx() * scPaRaan.sin() +
219 original.getEquinoctialEy() * scPaRaan.cos();
220 final double hx = original.getHx() * scRaan.cos() - original.getHy() * scRaan.sin();
221 final double hy = original.getHx() * scRaan.sin() + original.getHy() * scRaan.cos();
222 final double lambda = original.getLv() + dPaRaan;
223
224 // build updated orbit
225 final EquinoctialOrbit updated =
226 new EquinoctialOrbit(original.getA(), ex, ey, hx, hy, lambda, PositionAngleType.TRUE,
227 original.getFrame(), date, original.getMu());
228
229 // convert to required type
230 return orbit1.getType().convertType(updated);
231
232 }
233
234 }