FieldDeepSDP4.java
/* Copyright 2002-2023 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.analytical.tle;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.data.DataContext;
import org.orekit.frames.Frame;
import org.orekit.time.DateTimeComponents;
import org.orekit.utils.Constants;
/** This class contains the methods that compute deep space perturbation terms.
* <p>
* The user should not bother in this class since it is handled internaly by the
* {@link TLEPropagator}.
* </p>
* <p>This implementation is largely inspired from the paper and source code <a
* href="https://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
* Report #3</a> and is fully compliant with its results and tests cases.</p>
* @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
* @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
* @author Fabien Maussion (java translation)
* @author Thomas Paulet (field translation)
* @since 11.0
* @param <T> type of the field elements
*/
public class FieldDeepSDP4<T extends CalculusFieldElement<T>> extends FieldSDP4<T> {
// CHECKSTYLE: stop JavadocVariable check
/** Integration step (seconds). */
private static final double SECULAR_INTEGRATION_STEP = 720.0;
/** Intermediate values. */
private double thgr;
private T xnq;
private T omegaq;
private double zcosil;
private double zsinil;
private double zsinhl;
private double zcoshl;
private double zmol;
private double zcosgl;
private double zsingl;
private double zmos;
private T savtsn;
private T ee2;
private T e3;
private T xi2;
private T xi3;
private T xl2;
private T xl3;
private T xl4;
private T xgh2;
private T xgh3;
private T xgh4;
private T xh2;
private T xh3;
private T d2201;
private T d2211;
private T d3210;
private T d3222;
private T d4410;
private T d4422;
private T d5220;
private T d5232;
private T d5421;
private T d5433;
private T xlamo;
private T sse;
private T ssi;
private T ssl;
private T ssh;
private T ssg;
private T se2;
private T si2;
private T sl2;
private T sgh2;
private T sh2;
private T se3;
private T si3;
private T sl3;
private T sgh3;
private T sh3;
private T sl4;
private T sgh4;
private T del1;
private T del2;
private T del3;
private T xfact;
private T xli;
private T xni;
private T atime;
private T pe;
private T pinc;
private T pl;
private T pgh;
private T ph;
private T[] derivs;
// CHECKSTYLE: resume JavadocVariable check
/** Flag for resonant orbits. */
private boolean resonant;
/** Flag for synchronous orbits. */
private boolean synchronous;
/** Flag for compliance with Dundee modifications. */
private boolean isDundeeCompliant = true;
/** Constructor for a unique initial TLE.
*
* <p>This constructor uses the {@link DataContext#getDefault() default data context}.
*
* @param initialTLE the TLE to propagate.
* @param attitudeProvider provider for attitude computation
* @param mass spacecraft mass (kg)
* @param parameters SGP4 and SDP4 model parameters
* @see #FieldDeepSDP4(FieldTLE, AttitudeProvider, CalculusFieldElement, Frame, CalculusFieldElement[])
*/
@DefaultDataContext
public FieldDeepSDP4(final FieldTLE<T> initialTLE, final AttitudeProvider attitudeProvider,
final T mass, final T[] parameters) {
this(initialTLE, attitudeProvider, mass,
DataContext.getDefault().getFrames().getTEME(), parameters);
}
/** Constructor for a unique initial TLE.
* @param initialTLE the TLE to propagate.
* @param attitudeProvider provider for attitude computation
* @param mass spacecraft mass (kg)
* @param teme the TEME frame to use for propagation.
* @param parameters SGP4 and SDP4 model parameters
*/
public FieldDeepSDP4(final FieldTLE<T> initialTLE,
final AttitudeProvider attitudeProvider,
final T mass,
final Frame teme,
final T[] parameters) {
super(initialTLE, attitudeProvider, mass, teme, parameters);
}
/** Computes luni - solar terms from initial coordinates and epoch.
*/
protected void luniSolarTermsComputation() {
final T zero = tle.getPerigeeArgument().getField().getZero();
final T pi = zero.getPi();
final FieldSinCos<T> scg = FastMath.sinCos(tle.getPerigeeArgument());
final T sing = scg.sin();
final T cosg = scg.cos();
final FieldSinCos<T> scq = FastMath.sinCos(tle.getRaan());
final T sinq = scq.sin();
final T cosq = scq.cos();
final T aqnv = a0dp.reciprocal();
// Compute julian days since 1900
final double daysSince1900 = (tle.getDate()
.getComponents(utc)
.offsetFrom(DateTimeComponents.JULIAN_EPOCH)) /
Constants.JULIAN_DAY - 2415020;
double cc = TLEConstants.C1SS;
double ze = TLEConstants.ZES;
double zn = TLEConstants.ZNS;
T zsinh = sinq;
T zcosh = cosq;
thgr = thetaG(tle.getDate());
xnq = xn0dp;
omegaq = tle.getPerigeeArgument();
final double xnodce = 4.5236020 - 9.2422029e-4 * daysSince1900;
final SinCos scTem = FastMath.sinCos(xnodce);
final double stem = scTem.sin();
final double ctem = scTem.cos();
final double c_minus_gam = 0.228027132 * daysSince1900 - 1.1151842;
final double gam = 5.8351514 + 0.0019443680 * daysSince1900;
zcosil = 0.91375164 - 0.03568096 * ctem;
zsinil = FastMath.sqrt(1.0 - zcosil * zcosil);
zsinhl = 0.089683511 * stem / zsinil;
zcoshl = FastMath.sqrt(1.0 - zsinhl * zsinhl);
zmol = MathUtils.normalizeAngle(c_minus_gam, pi.getReal());
double zx = 0.39785416 * stem / zsinil;
final double zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
zx = FastMath.atan2( zx, zy) + gam - xnodce;
final SinCos scZx = FastMath.sinCos(zx);
zcosgl = scZx.cos();
zsingl = scZx.sin();
zmos = MathUtils.normalizeAngle(6.2565837 + 0.017201977 * daysSince1900, pi.getReal());
// Do solar terms
savtsn = zero.add(1e20);
T zcosi = zero.add(0.91744867);
T zsini = zero.add(0.39785416);
T zsing = zero.add(-0.98088458);
T zcosg = zero.add(0.1945905);
T se = zero;
T sgh = zero;
T sh = zero;
T si = zero;
T sl = zero;
// There was previously some convoluted logic here, but it boils
// down to this: we compute the solar terms, then the lunar terms.
// On a second pass, we recompute the solar terms, taking advantage
// of the improved data that resulted from computing lunar terms.
for (int iteration = 0; iteration < 2; ++iteration) {
final T a1 = zcosh.multiply(zcosg).add(zsinh.multiply(zsing).multiply(zcosi));
final T a3 = zcosh.multiply(zsing.negate()).add(zsinh.multiply(zcosg).multiply(zcosi));
final T a7 = zsinh.negate().multiply(zcosg).add(zcosh.multiply(zcosi).multiply(zsing));
final T a8 = zsing.multiply(zsini);
final T a9 = zsinh.multiply(zsing).add(zcosh.multiply(zcosi).multiply(zcosg));
final T a10 = zcosg.multiply(zsini);
final T a2 = cosi0.multiply(a7).add(sini0.multiply(a8));
final T a4 = cosi0.multiply(a9).add(sini0.multiply(a10));
final T a5 = sini0.negate().multiply(a7).add(cosi0.multiply(a8));
final T a6 = sini0.negate().multiply(a9).add(cosi0.multiply(a10));
final T x1 = a1.multiply(cosg).add(a2.multiply(sing));
final T x2 = a3.multiply(cosg).add(a4.multiply(sing));
final T x3 = a1.negate().multiply(sing).add(a2.multiply(cosg));
final T x4 = a3.negate().multiply(sing).add(a4.multiply(cosg));
final T x5 = a5.multiply(sing);
final T x6 = a6.multiply(sing);
final T x7 = a5.multiply(cosg);
final T x8 = a6.multiply(cosg);
final T z31 = x1.multiply(x1).multiply(12).subtract(x3.multiply(x3).multiply(3));
final T z32 = x1.multiply(x2).multiply(24).subtract(x3.multiply(x4).multiply(6));
final T z33 = x2.multiply(x2).multiply(12).subtract(x4.multiply(x4).multiply(3));
final T z11 = a1.multiply(-6).multiply(a5).add(e0sq.multiply(x1.multiply(x7).multiply(-24).add(x3.multiply(x5).multiply(-6))));
final T z12 = a1.multiply(a6).add(a3.multiply(a5)).multiply(-6).add(
e0sq.multiply(x2.multiply(x7).add(x1.multiply(x8)).multiply(-24).add(
x3.multiply(x6).add(x4.multiply(x5)).multiply(-6))));
final T z13 = a3.multiply(a6).multiply(-6).add(e0sq.multiply(
x2.multiply(x8).multiply(-24).subtract(x4.multiply(x6).multiply(6))));
final T z21 = a2.multiply(a5).multiply(6).add(e0sq.multiply(
x1.multiply(x5).multiply(24).subtract(x3.multiply(x7).multiply(6))));
final T z22 = a4.multiply(a5).add(a2.multiply(a6)).multiply(6).add(
e0sq.multiply(x2.multiply(x5).add(x1.multiply(x6)).multiply(24).subtract(
x4.multiply(x7).add(x3.multiply(x8)).multiply(6))));
final T z23 = a4.multiply(a6).multiply(6).add(e0sq.multiply(x2.multiply(x6).multiply(24).subtract(x4.multiply(x8).multiply(6))));
final T s3 = xnq.reciprocal().multiply(cc);
final T s2 = beta0.reciprocal().multiply(s3.multiply(-0.5));
final T s4 = s3.multiply(beta0);
final T s1 = tle.getE().multiply(s4).multiply(-15);
final T s5 = x1.multiply(x3).add(x2.multiply(x4));
final T s6 = x2.multiply(x3).add(x1.multiply(x4));
final T s7 = x2.multiply(x4).subtract(x1.multiply(x3));
T z1 = a1.multiply(a1).add(a2.multiply(a2)).multiply(3).add(z31.multiply(e0sq));
T z2 = a1.multiply(a3).add(a2.multiply(a4)).multiply(6).add(z32.multiply(e0sq));
T z3 = a3.multiply(a3).add(a4.multiply(a4)).multiply(3).add(z33.multiply(e0sq));
z1 = z1.add(z1).add(beta02.multiply(z31));
z2 = z2.add(z2).add(beta02.multiply(z32));
z3 = z3.add(z3).add(beta02.multiply(z33));
se = s1.multiply(zn).multiply(s5);
si = s2.multiply(zn).multiply(z11.add(z13));
sl = s3.multiply(-zn).multiply(z1.add(z3).subtract(14).subtract(e0sq.multiply(6)));
sgh = s4.multiply(zn).multiply(z31.add(z33).subtract(6));
if (tle.getI().getReal() < pi.divide(60.0).getReal()) {
// inclination smaller than 3 degrees
sh = zero;
} else {
sh = s2.multiply(-zn).multiply(z21.add(z23));
}
ee2 = s1.multiply(s6).multiply(2);
e3 = s1.multiply(s7).multiply(2);
xi2 = s2.multiply(z12).multiply(2);
xi3 = s2.multiply(z13.subtract(z11)).multiply(2);
xl2 = s3.multiply(z2).multiply(-2);
xl3 = s3.multiply(z3.subtract(z1)).multiply(-2);
xl4 = s3.multiply(e0sq.multiply(-9).add(-21)).multiply(ze).multiply(-2);
xgh2 = s4.multiply(z32).multiply(2);
xgh3 = s4.multiply(z33.subtract(z31)).multiply(2);
xgh4 = s4.multiply(ze).multiply(-18);
xh2 = s2.multiply(z22).multiply(-2);
xh3 = s2.multiply(z23.subtract(z21)).multiply(-2);
if (iteration == 0) { // we compute lunar terms only on the first pass:
sse = se;
ssi = si;
ssl = sl;
ssh = (tle.getI().getReal() < pi.divide(60.0).getReal()) ? zero : sh.divide(sini0);
ssg = sgh.subtract(cosi0.multiply(ssh));
se2 = ee2;
si2 = xi2;
sl2 = xl2;
sgh2 = xgh2;
sh2 = xh2;
se3 = e3;
si3 = xi3;
sl3 = xl3;
sgh3 = xgh3;
sh3 = xh3;
sl4 = xl4;
sgh4 = xgh4;
zcosg = zero.add(zcosgl);
zsing = zero.add(zsingl);
zcosi = zero.add(zcosil);
zsini = zero.add(zsinil);
zcosh = cosq.multiply(zcoshl).add(sinq.multiply(zsinhl));
zsinh = sinq.multiply(zcoshl).subtract(cosq.multiply(zsinhl));
zn = TLEConstants.ZNL;
cc = TLEConstants.C1L;
ze = TLEConstants.ZEL;
}
} // end of solar - lunar - solar terms computation
sse = sse.add(se);
ssi = ssi.add(si);
ssl = ssl.add(sl);
ssg = ssg.add(sgh).subtract((tle.getI().getReal() < pi.divide(60.0).getReal()) ? zero : (cosi0.divide(sini0).multiply(sh)));
ssh = ssh.add((tle.getI().getReal() < pi.divide(60.0).getReal()) ? zero : sh.divide(sini0));
// Start the resonant-synchronous tests and initialization
T bfact = zero;
// if mean motion is 1.893053 to 2.117652 revs/day, and eccentricity >= 0.5,
// start of the 12-hour orbit, e > 0.5 section
if (xnq.getReal() >= 0.00826 && xnq.getReal() <= 0.00924 && tle.getE().getReal() >= 0.5) {
final T g201 = tle.getE().subtract(0.64).negate().multiply(0.440).add(-0.306);
final T eoc = tle.getE().multiply(e0sq);
final T sini2 = sini0.multiply(sini0);
final T f220 = cosi0.multiply(2).add(theta2).add(1).multiply(0.75);
final T f221 = sini2.multiply(1.5);
final T f321 = sini0.multiply(1.875).multiply(cosi0.multiply(2).negate().subtract(theta2.multiply(3)).add(1));
final T f322 = sini0.multiply(-1.875).multiply(cosi0.multiply(2).subtract(theta2.multiply(3)).add(1));
final T f441 = sini2.multiply(f220).multiply(35);
final T f442 = sini2.multiply(sini2).multiply(39.3750);
final T f522 = sini0.multiply(9.84375).multiply(sini2.multiply(cosi0.multiply(-2).add(theta2.multiply(-5)).add(1.0)).add(
cosi0.multiply(4.0).add(theta2.multiply(6.0)).add(-2).multiply(0.33333333)));
final T f523 = sini0.multiply(sini2.multiply(cosi0.multiply(-4).add(theta2.multiply(10)).add(-2)).multiply(4.92187512).add(
cosi0.multiply(2).subtract(theta2.multiply(3)).add(1).multiply(6.56250012)));
final T f542 = sini0.multiply(29.53125).multiply(cosi0.multiply(-8).add(2).add(
theta2.multiply(cosi0.multiply(8).add(theta2.multiply(10)).add(-12))));
final T f543 = sini0.multiply(29.53125).multiply(cosi0.multiply(-8).add(-2).add(
theta2.multiply(cosi0.multiply(8).subtract(theta2.multiply(10)).add(12))));
final T g211;
final T g310;
final T g322;
final T g410;
final T g422;
final T g520;
resonant = true; // it is resonant...
synchronous = false; // but it's not synchronous
// Geopotential resonance initialization for 12 hour orbits :
if (tle.getE().getReal() <= 0.65) {
g211 = tle.getE().multiply( -13.247).add( e0sq.multiply( 16.290)).add( 3.616);
g310 = tle.getE().multiply( 117.390).add( e0sq.multiply( -228.419)).add( eoc.multiply( 156.591)).add( -19.302);
g322 = tle.getE().multiply(109.7927).add( e0sq.multiply(-214.6334)).add( eoc.multiply(146.5816)).add( -18.9068);
g410 = tle.getE().multiply( 242.694).add( e0sq.multiply( -471.094)).add( eoc.multiply( 313.953)).add( -41.122);
g422 = tle.getE().multiply( 841.880).add( e0sq.multiply(-1629.014)).add( eoc.multiply(1083.435)).add( -146.407);
g520 = tle.getE().multiply(3017.977).add( e0sq.multiply(-5740.032)).add( eoc.multiply(3708.276)).add( -532.114);
} else {
g211 = tle.getE().multiply( 331.819).add( e0sq.multiply( -508.738)).add( eoc.multiply( 266.724)).add( -72.099);
g310 = tle.getE().multiply(1582.851).add( e0sq.multiply(-2415.925)).add( eoc.multiply(1246.113)).add( -346.844);
g322 = tle.getE().multiply(1554.908).add( e0sq.multiply(-2366.899)).add( eoc.multiply(1215.972)).add( -342.585);
g410 = tle.getE().multiply(4758.686).add( e0sq.multiply(-7193.992)).add( eoc.multiply(3651.957)).add(-1052.797);
g422 = tle.getE().multiply(16178.11).add( e0sq.multiply(-24462.77)).add( eoc.multiply(12422.52)).add( -3581.69);
if (tle.getE().getReal() <= 0.715) {
g520 = tle.getE().multiply(-4664.75).add( e0sq.multiply( 3763.64)).add( 1464.74);
} else {
g520 = tle.getE().multiply(29936.92).add( e0sq.multiply(-54087.36)).add( eoc.multiply(31324.56)).add( -5149.66);
}
}
final T g533;
final T g521;
final T g532;
if (tle.getE().getReal() < 0.7) {
g533 = tle.getE().multiply( 4988.61).add( e0sq.multiply( -9064.77)).add( eoc.multiply( 5542.21)).add( -919.2277);
g521 = tle.getE().multiply(4568.6173).add( e0sq.multiply(-8491.4146)).add( eoc.multiply( 5337.524)).add( -822.71072);
g532 = tle.getE().multiply( 4690.25).add( e0sq.multiply( -8624.77)).add( eoc.multiply( 5341.4)).add( -853.666);
} else {
g533 = tle.getE().multiply(161616.52).add( e0sq.multiply( -229838.2)).add( eoc.multiply(109377.94)).add( -37995.78);
g521 = tle.getE().multiply(218913.95).add( e0sq.multiply(-309468.16)).add( eoc.multiply(146349.42)).add( -51752.104);
g532 = tle.getE().multiply(170470.89).add( e0sq.multiply(-242699.48)).add( eoc.multiply(115605.82)).add( -40023.88);
}
T temp1 = xnq.multiply(xnq).multiply(aqnv).multiply(aqnv).multiply(3);
T temp = temp1.multiply(TLEConstants.ROOT22);
d2201 = temp.multiply(f220).multiply(g201);
d2211 = temp.multiply(f221).multiply(g211);
temp1 = temp1.multiply(aqnv);
temp = temp1.multiply(TLEConstants.ROOT32);
d3210 = temp.multiply(f321).multiply(g310);
d3222 = temp.multiply(f322).multiply(g322);
temp1 = temp1.multiply(aqnv);
temp = temp1.multiply(2 * TLEConstants.ROOT44);
d4410 = temp.multiply(f441).multiply(g410);
d4422 = temp.multiply(f442).multiply(g422);
temp1 = temp1.multiply(aqnv);
temp = temp1.multiply(TLEConstants.ROOT52);
d5220 = temp.multiply(f522).multiply(g520);
d5232 = temp.multiply(f523).multiply(g532);
temp = temp1.multiply(2 * TLEConstants.ROOT54);
d5421 = temp.multiply(f542).multiply(g521);
d5433 = temp.multiply(f543).multiply(g533);
xlamo = tle.getMeanAnomaly().add(tle.getRaan()).add(tle.getRaan()).subtract(thgr + thgr);
bfact = xmdot.add(xnodot).add(xnodot).subtract(TLEConstants.THDT + TLEConstants.THDT);
bfact = bfact.add(ssl).add(ssh).add(ssh);
} else if (xnq.getReal() < 0.0052359877 && xnq.getReal() > 0.0034906585) {
// if mean motion is .8 to 1.2 revs/day : (geosynch)
final T cosio_plus_1 = cosi0.add(1.0);
final T g200 = e0sq.multiply(e0sq.multiply(0.8125).add(-2.5)).add(1);
final T g300 = e0sq.multiply(e0sq.multiply(6.60937).add(-6)).add(1);
final T f311 = sini0.multiply(0.9375).multiply(sini0.multiply(cosi0.multiply(3).add(1))).subtract(cosio_plus_1.multiply(0.75));
final T g310 = e0sq.multiply(2).add(1);
final T f220 = cosio_plus_1.multiply(cosio_plus_1).multiply(0.75);
final T f330 = f220.multiply(cosio_plus_1).multiply(2.5);
resonant = true;
synchronous = true;
// Synchronous resonance terms initialization
del1 = xnq.multiply(xnq).multiply(aqnv).multiply(aqnv).multiply(3);
del2 = del1.multiply(f220).multiply(g200).multiply(2 * TLEConstants.Q22);
del3 = del1.multiply(f330).multiply(g300).multiply(aqnv).multiply(3 * TLEConstants.Q33);
del1 = del1.multiply(f311).multiply(g310).multiply(TLEConstants.Q31).multiply(aqnv);
xlamo = tle.getMeanAnomaly().add(tle.getRaan()).add(tle.getPerigeeArgument()).subtract(thgr);
bfact = xmdot.add(omgdot).add(xnodot).subtract(TLEConstants.THDT);
bfact = bfact.add(ssl).add(ssg).add(ssh);
} else {
// it's neither a high-e 12-hours orbit nor a geosynchronous:
resonant = false;
synchronous = false;
}
if (resonant) {
xfact = bfact.subtract(xnq);
// Initialize integrator
xli = xlamo;
xni = xnq;
atime = zero;
}
derivs = MathArrays.buildArray(xnq.getField(), 2);
}
/** Computes secular terms from current coordinates and epoch.
* @param t offset from initial epoch (minutes)
*/
protected void deepSecularEffects(final T t) {
xll = xll.add(ssl.multiply(t));
omgadf = omgadf.add(ssg.multiply(t));
xnode = xnode.add(ssh.multiply(t));
em = tle.getE().add(sse.multiply(t));
xinc = tle.getI().add(ssi.multiply(t));
if (resonant) {
// If we're closer to t = 0 than to the currently-stored data
// from the previous call to this function, then we're
// better off "restarting", going back to the initial data.
// The Dundee code rigs things up to _always_ take 720-minute
// steps from epoch to end time, except for the final step.
// Easiest way to arrange similar behavior in this code is
// just to always do a restart, if we're in Dundee-compliant
// mode.
if (FastMath.abs(t).getReal() < FastMath.abs(t.subtract(atime)).getReal() || isDundeeCompliant) {
// Epoch restart
atime = t.getField().getZero();
xni = xnq;
xli = xlamo;
}
boolean lastIntegrationStep = false;
// if |step|>|step max| then do one step at step max
while (!lastIntegrationStep) {
double delt = t.subtract(atime).getReal();
if (delt > SECULAR_INTEGRATION_STEP) {
delt = SECULAR_INTEGRATION_STEP;
} else if (delt < -SECULAR_INTEGRATION_STEP) {
delt = -SECULAR_INTEGRATION_STEP;
} else {
lastIntegrationStep = true;
}
computeSecularDerivs();
final T xldot = xni.add(xfact);
T xlpow = t.getField().getZero().add(1.);
xli = xli.add(xldot.multiply(delt));
xni = xni.add(derivs[0].multiply(delt));
double delt_factor = delt;
xlpow = xlpow.multiply(xldot);
derivs[1] = derivs[1].multiply(xlpow);
delt_factor *= delt / 2;
xli = xli.add(derivs[0].multiply(delt_factor));
xni = xni.add(derivs[1].multiply(delt_factor));
atime = atime.add(delt);
}
xn = xni;
final T temp = xnode.negate().add(thgr).add(t.multiply(TLEConstants.THDT));
xll = xli.add(temp).add(synchronous ? omgadf.negate() : temp);
}
}
/** Computes periodic terms from current coordinates and epoch.
* @param t offset from initial epoch (min)
*/
protected void deepPeriodicEffects(final T t) {
// If the time didn't change by more than 30 minutes,
// there's no good reason to recompute the perturbations;
// they don't change enough over so short a time span.
// However, the Dundee code _always_ recomputes, so if
// we're attempting to replicate its results, we've gotta
// recompute everything, too.
if (FastMath.abs(savtsn.subtract(t).getReal()) >= 30.0 || isDundeeCompliant) {
savtsn = t;
// Update solar perturbations for time T
T zm = t.multiply(TLEConstants.ZNS).add(zmos);
T zf = zm.add(FastMath.sin(zm).multiply(2 * TLEConstants.ZES));
FieldSinCos<T> sczf = FastMath.sinCos(zf);
T sinzf = sczf.sin();
T f2 = sinzf.multiply(sinzf).multiply(0.5).subtract(0.25);
T f3 = sinzf.multiply(sczf.cos()).multiply(-0.5);
final T ses = se2.multiply(f2).add(se3.multiply(f3));
final T sis = si2.multiply(f2).add(si3.multiply(f3));
final T sls = sl2.multiply(f2).add(sl3.multiply(f3)).add(sl4.multiply(sinzf));
final T sghs = sgh2.multiply(f2).add(sgh3.multiply(f3)).add(sgh4.multiply(sinzf));
final T shs = sh2.multiply(f2).add(sh3.multiply(f3));
// Update lunar perturbations for time T
zm = t.multiply(TLEConstants.ZNL).add(zmol);
zf = zm.add(FastMath.sin(zm).multiply(2 * TLEConstants.ZEL));
sczf = FastMath.sinCos(zf);
sinzf = sczf.sin();
f2 = sinzf.multiply(sinzf).multiply(0.5).subtract(0.25);
f3 = sinzf.multiply(sczf.cos()).multiply(-0.5);
final T sel = ee2.multiply(f2).add(e3.multiply(f3));
final T sil = xi2.multiply(f2).add(xi3.multiply(f3));
final T sll = xl2.multiply(f2).add(xl3.multiply(f3)).add(xl4.multiply(sinzf));
final T sghl = xgh2.multiply(f2).add(xgh3.multiply(f3)).add(xgh4.multiply(sinzf));
final T sh1 = xh2.multiply(f2).add(xh3.multiply(f3));
// Sum the solar and lunar contributions
pe = ses.add(sel);
pinc = sis.add(sil);
pl = sls.add(sll);
pgh = sghs.add(sghl);
ph = shs.add(sh1);
}
xinc = xinc.add(pinc);
final FieldSinCos<T> scis = FastMath.sinCos(xinc);
final T sinis = scis.sin();
final T cosis = scis.cos();
/* Add solar/lunar perturbation correction to eccentricity: */
em = em.add(pe);
xll = xll.add(pl);
omgadf = omgadf.add(pgh);
xinc = MathUtils.normalizeAngle(xinc, t.getField().getZero());
if (FastMath.abs(xinc).getReal() >= 0.2) {
// Apply periodics directly
final T temp_val = ph.divide(sinis);
omgadf = omgadf.subtract(cosis.multiply(temp_val));
xnode = xnode.add(temp_val);
} else {
// Apply periodics with Lyddane modification
final FieldSinCos<T> scok = FastMath.sinCos(xnode);
final T sinok = scok.sin();
final T cosok = scok.cos();
final T alfdp = ph.multiply(cosok).add((pinc.multiply(cosis).add(sinis)).multiply(sinok));
final T betdp = ph.negate().multiply(sinok).add((pinc.multiply(cosis).add(sinis)).multiply(cosok));
final T delta_xnode = MathUtils.normalizeAngle(FastMath.atan2(alfdp, betdp).subtract(xnode), t.getField().getZero());
final T dls = xnode.negate().multiply(sinis).multiply(pinc);
omgadf = omgadf.add(dls.subtract(cosis.multiply(delta_xnode)));
xnode = xnode.add(delta_xnode);
}
}
/** Computes internal secular derivs. */
private void computeSecularDerivs() {
final FieldSinCos<T> sc_li = FastMath.sinCos(xli);
final T sin_li = sc_li.sin();
final T cos_li = sc_li.cos();
final T sin_2li = sin_li.multiply(cos_li).multiply(2.);
final T cos_2li = cos_li.multiply(cos_li).multiply(2.).subtract(1.);
// Dot terms calculated :
if (synchronous) {
final T sin_3li = sin_2li.multiply(cos_li).add(cos_2li.multiply(sin_li));
final T cos_3li = cos_2li.multiply(cos_li).subtract(sin_2li.multiply(sin_li));
final T term1a = del1.multiply(sin_li .multiply(TLEConstants.C_FASX2) .subtract(cos_li .multiply(TLEConstants.S_FASX2 )));
final T term2a = del2.multiply(sin_2li.multiply(TLEConstants.C_2FASX4).subtract(cos_2li.multiply(TLEConstants.S_2FASX4)));
final T term3a = del3.multiply(sin_3li.multiply(TLEConstants.C_3FASX6).subtract(cos_3li.multiply(TLEConstants.S_3FASX6)));
final T term1b = del1.multiply(cos_li .multiply(TLEConstants.C_FASX2) .add(sin_li .multiply(TLEConstants.S_FASX2 )));
final T term2b = del2.multiply(cos_2li.multiply(TLEConstants.C_2FASX4) .add(sin_2li.multiply(TLEConstants.S_2FASX4))).multiply(2.0);
final T term3b = del3.multiply(cos_3li.multiply(TLEConstants.C_3FASX6) .add(sin_3li.multiply(TLEConstants.S_3FASX6))).multiply(3.0);
derivs[0] = term1a.add(term2a).add(term3a);
derivs[1] = term1b.add(term2b).add(term3b);
} else {
// orbit is a 12-hour resonant one
final T xomi = omegaq.add(omgdot.multiply(atime));
final FieldSinCos<T> sc_omi = FastMath.sinCos(xomi);
final T sin_omi = sc_omi.sin();
final T cos_omi = sc_omi.cos();
final T sin_li_m_omi = sin_li.multiply(cos_omi).subtract(sin_omi.multiply(cos_li));
final T sin_li_p_omi = sin_li.multiply(cos_omi).add( sin_omi.multiply(cos_li));
final T cos_li_m_omi = cos_li.multiply(cos_omi).add( sin_omi.multiply(sin_li));
final T cos_li_p_omi = cos_li.multiply(cos_omi).subtract(sin_omi.multiply(sin_li));
final T sin_2omi = sin_omi.multiply(cos_omi).multiply(2.0);
final T cos_2omi = cos_omi.multiply(cos_omi).multiply(2.0).subtract(1.0);
final T sin_2li_m_omi = sin_2li.multiply(cos_omi ).subtract(sin_omi .multiply(cos_2li));
final T sin_2li_p_omi = sin_2li.multiply(cos_omi ).add( sin_omi .multiply(cos_2li));
final T cos_2li_m_omi = cos_2li.multiply(cos_omi ).add( sin_omi .multiply(sin_2li));
final T cos_2li_p_omi = cos_2li.multiply(cos_omi ).subtract(sin_omi .multiply(sin_2li));
final T sin_2li_p_2omi = sin_2li.multiply(cos_2omi).add( sin_2omi.multiply(cos_2li));
final T cos_2li_p_2omi = cos_2li.multiply(cos_2omi).subtract(sin_2omi.multiply(sin_2li));
final T sin_2omi_p_li = sin_li .multiply(cos_2omi).add( sin_2omi.multiply(cos_li ));
final T cos_2omi_p_li = cos_li .multiply(cos_2omi).subtract(sin_2omi.multiply(sin_li ));
final T term1a = d2201.multiply(sin_2omi_p_li .multiply(TLEConstants.C_G22).subtract(cos_2omi_p_li .multiply(TLEConstants.S_G22))) .add(
d2211.multiply(sin_li .multiply(TLEConstants.C_G22).subtract(cos_li .multiply(TLEConstants.S_G22)))).add(
d3210.multiply(sin_li_p_omi .multiply(TLEConstants.C_G32).subtract(cos_li_p_omi .multiply(TLEConstants.S_G32)))).add(
d3222.multiply(sin_li_m_omi .multiply(TLEConstants.C_G32).subtract(cos_li_m_omi .multiply(TLEConstants.S_G32)))).add(
d5220.multiply(sin_li_p_omi .multiply(TLEConstants.C_G52).subtract(cos_li_p_omi .multiply(TLEConstants.S_G52)))).add(
d5232.multiply(sin_li_m_omi .multiply(TLEConstants.C_G52).subtract(cos_li_m_omi .multiply(TLEConstants.S_G52))));
final T term2a = d4410.multiply(sin_2li_p_2omi.multiply(TLEConstants.C_G44).subtract(cos_2li_p_2omi.multiply(TLEConstants.S_G44))) .add(
d4422.multiply(sin_2li .multiply(TLEConstants.C_G44).subtract(cos_2li .multiply(TLEConstants.S_G44)))).add(
d5421.multiply(sin_2li_p_omi .multiply(TLEConstants.C_G54).subtract(cos_2li_p_omi .multiply(TLEConstants.S_G54)))).add(
d5433.multiply(sin_2li_m_omi .multiply(TLEConstants.C_G54).subtract(cos_2li_m_omi .multiply(TLEConstants.S_G54))));
final T term1b = d2201.multiply(cos_2omi_p_li .multiply(TLEConstants.C_G22) .add(sin_2omi_p_li .multiply(TLEConstants.S_G22))) .add(
d2211.multiply(cos_li .multiply(TLEConstants.C_G22) .add(sin_li .multiply(TLEConstants.S_G22)))).add(
d3210.multiply(cos_li_p_omi .multiply(TLEConstants.C_G32) .add(sin_li_p_omi .multiply(TLEConstants.S_G32)))).add(
d3222.multiply(cos_li_m_omi .multiply(TLEConstants.C_G32) .add(sin_li_m_omi .multiply(TLEConstants.S_G32)))).add(
d5220.multiply(cos_li_p_omi .multiply(TLEConstants.C_G52) .add(sin_li_p_omi .multiply(TLEConstants.S_G52)))).add(
d5232.multiply(cos_li_m_omi .multiply(TLEConstants.C_G52) .add(sin_li_m_omi .multiply(TLEConstants.S_G52))));
final T term2b = d4410.multiply(cos_2li_p_2omi.multiply(TLEConstants.C_G44) .add(sin_2li_p_2omi.multiply(TLEConstants.S_G44))) .add(
d4422.multiply(cos_2li .multiply(TLEConstants.C_G44) .add(sin_2li .multiply(TLEConstants.S_G44)))).add(
d5421.multiply(cos_2li_p_omi .multiply(TLEConstants.C_G54) .add(sin_2li_p_omi .multiply(TLEConstants.S_G54)))).add(
d5433.multiply(cos_2li_m_omi .multiply(TLEConstants.C_G54) .add(sin_2li_m_omi .multiply(TLEConstants.S_G54)))).multiply(2.0);
derivs[0] = term1a.add(term2a);
derivs[1] = term1b.add(term2b);
}
}
}