非惯性参考系中的轨道传播

本教程的目标是介绍使用SingleBodyAttraction类进行轨道积分。
这个类可以替代所有不同类型的质点相互作用（ThirdBodyAttraction，NewtonianAttraction）。
使用SingleBodyAttraction和InertialForces可以实现更丰富的建模，允许用户在一个不一定以主要吸引体为中心且不一定具有惯性轴的参考系中计算卫星的运动。

初始化

我们将传播一个从地月拉格朗日点L2出发的卫星的轨迹。
运动方程将在以L2为中心的参考系中计算，其X轴不断指向地球和月球。
这个参考系的旋转意味着我们不能简单地使用动力学的基本原理（牛顿第二定律）。
由于我们只考虑地球和月球的引力作用，与我们的问题最相关的惯性参考系是一个以地月质心为中心、具有惯性轴的参考系。

让我们初始化我们的程序。首先是时间设置：积分的初始时刻、积分的长度（这里以秒为单位）以及每个输出之间的时间间隔。

final AbsoluteDate initialDate = new AbsoluteDate(2000, 01, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
double integrationTime = 600000.;
double outputStep = 600.0;

然后是卫星的初始位置和速度，位于L2点。我们还没有引入参考系，所以用户必须记住在哪个参考系中定义其卫星的初始条件。

final PVCoordinates initialConditions = new PVCoordinates(new Vector3D(0.0, 0.0, 0.0), new Vector3D(0.0, 0.0, 0.0));

参考系设置

我们需要加载一些天体，这些天体与我们的卫星在引力相互作用中考虑在内，以及地月质心，因为我们将根据其附加的参考系计算惯性力。然后我们创建与L2点相关联的参考系，以及与地月质心相关联的惯性参考系。

final CelestialBody earth = CelestialBodyFactory.getEarth();
final CelestialBody moon  = CelestialBodyFactory.getMoon();
final CelestialBody earthMoonBary = CelestialBodyFactory.getEarthMoonBarycenter();

final Frame l2Frame = new L2Frame(earth, moon);
final Frame earthMoonBaryFrame = earthMoonBary.getInertiallyOrientedFrame();

final Frame inertiaFrame = earthMoonBaryFrame;
final Frame integrationFrame = l2Frame;
final Frame outputFrame = l2Frame;

传播器准备工作

我们现在将我们的PVCoordinates转换为AbsolutePVCoordinates，定义卫星姿态，并使用所有这些来构建SpacecraftState

final AbsolutePVCoordinates initialAbsPV = new AbsolutePVCoordinates(integrationFrame, initialDate, initialConditions);
Attitude arbitraryAttitude = new Attitude(integrationFrame,
                                          new TimeStampedAngularCoordinates(initialDate,
                                          new PVCoordinates(Vector3D.PLUS_I, Vector3D.PLUS_J),
                                          new PVCoordinates(Vector3D.PLUS_I, Vector3D.PLUS_J)));
final SpacecraftState initialState = new SpacecraftState(initialAbsPV, arbitraryAttitude);

我们将使用一个变步长的8(5,3) Dormand-Prince积分器。这个积分器需要一些初始化参数。首先，让我们设置积分步长的边界，以防止计算时间过长或积分步长过长。

final double minStep = 0.001;
final double maxstep = 3600.0;

我们还需要设置积分的可接受误差。这个误差用于调整积分步长，并不代表传播的整体误差。

final double positionTolerance = 0.001;
final double velocityTolerance = 0.00001;
final double massTolerance     = 1.0e-6;
final double[] vecAbsoluteTolerances = {
        positionTolerance, positionTolerance, positionTolerance,
        velocityTolerance, velocityTolerance, velocityTolerance,
        massTolerance
        };
final double[] vecRelativeTolerances = new double[vecAbsoluteTolerances.length];

现在我们可以定义传播器的数值积分器。

AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep,
                                                                       vecAbsoluteTolerances,
                                                                       vecRelativeTolerances);

构建和使用传播器

最后，我们可以构建传播器。首先，我们使用之前定义的数值积分器来构造NumericalPropagator，然后我们指定演化模型的属性。使用SingleBodyAttraction意味着我们不限制轨道类型，并且我们需要使用setOrbitType(null)和setIgnoreCentralAttraction(true)。我们的积分参考系的非惯性方面由InertialForces处理，它将计算出现在该参考系中的惯性加速度。每个引力作用的吸引体都将通过调用SingleBodyAttraction来提供给传播器。然后，我们只需要设置初始的航天器状态和传播器模式，有关传播器模式的更多详细信息可以在传播教程中找到。

NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(null);
propagator.setIgnoreCentralAttraction(true);
propagator.addForceModel(new InertialForces(earthMoonBaryFrame));
propagator.addForceModel(new SingleBodyAbsoluteAttraction(earth));
propagator.addForceModel(new SingleBodyAbsoluteAttraction(moon));
propagator.setInitialState(initialState);
propagator.getMultiplexer().add(outputStep, new TutorialStepHandler("test.dat", outputFrame));

最后，我们可以进行传播本身，从初始日期开始，持续一定的时间。

SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(integrationTime));
final PVCoordinates pv = finalState.getPVCoordinates(outputFrame);
System.out.println("初始条件: " + initialConditions);
System.out.println("最终条件: " + pv);

步骤处理

现在，我们可以处理我们的步骤处理器，它将在每个积分输出步骤时打印位置、速度和加速度。

 private static class TutorialStepHandler implements OrekitFixedStepHandler {
 private Frame outputFrame;
 private TutorialStepHandler( final Frame frame) {
    outputFrame = frame;
    }

 public void init(final SpacecraftState s0, final AbsoluteDate t) {
        System.out.format(Locale.US,
                                  "%s %s %s %s %s %s %s %s %s %s %n",
                                  "日期", "                           X", "                 Y",
                                  "                 Z", "                 Vx", "                Vy",
                                  "                Vz", "                ax", "                ay",
                                  "                az");
    }

    public void handleStep(SpacecraftState currentState) {
        final AbsoluteDate d = currentState.getDate();
        final PVCoordinates pv = currentState.getPVCoordinates(outputFrame);
        System.out.format(Locale.US,
                          "%s %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f %18.12f%n",
                          d, pv.getPosition().getX(),
                          pv.getPosition().getY(), pv.getPosition().getZ(),
                          pv.getVelocity().getX(), pv.getVelocity().getY(),
                          pv.getVelocity().getZ(), pv.getAcceleration().getX(),
                          pv.getAcceleration().getY(),
                          pv.getAcceleration().getZ());
    }

    public void finish(final SpacecraftState finalState) {
        final PVCoordinates finalPv =
                        finalState.getPVCoordinates(outputFrame);
        System.out.println();
        System.out.format(Locale.US,
                          "%s %12.0f %12.0f %12.0f %12.0f %12.0f %12.0f%n",
                          finalState.getDate(), finalPv.getPosition().getX(),
                          finalPv.getPosition().getY(),
                          finalPv.getPosition().getZ(),
                          finalPv.getVelocity().getX(),
                          finalPv.getVelocity().getY(),
                          finalPv.getVelocity().getZ());
        System.out.println();
    }

}

结果

打印的结果如下所示：

日期 X Y Z Vx Vy Vz ax ay az 2000-01-01T00:00:00.000Z 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000007203264 -0.000023617334 0.000000988884 2000-01-01T00:10:00.000Z 1.310181730693 -4.509937975829 0.177995372228 0.004389441334 -0.015464490879 0.000593311180 0.000007424953 -0.000027930976 0.000000988817 2000-01-01T00:20:00.000Z 5.292760978714 -19.075027472324 0.711964496463 0.008906039665 -0.033517172720 0.001186577761 0.000007627122 -0.000032244627 0.000000988736 2000-01-01T00:30:00.000Z 12.020527119451 -45.248282259272 1.601878011878 0.013538105089 -0.054158285678 0.001779791113 0.000007690269 -0.000035156300 0.000000989424 2000-01-01T00:40:00.000Z 21.178259937228 -80.021547310428 2.850254956138 0.017011615221 -0.062141479011 0.002381481685 0.000005918995 -0.000015387558 0.000001002821 2000-01-01T00:50:00.000Z 32.466255913066 -120.334453824629 4.459646951922 0.020640563371 -0.072665310900 0.002983149974 0.000006174249 -0.000019691905 0.000001002738 2000-01-01T01:00:00.000Z 45.976396611751 -167.736449656047 6.430024082490 0.024416814768 -0.085771787993 0.003584764051 0.000006410007 -0.000023996366 0.000001002640 2000-01-01T01:10:00.000Z 61.793555128366 -223.777138815637 8.761351268544 0.028328673735 -0.101460960959 0.004186315403 0.000006626275 -0.000028300882 0.000001002529 2000-01-01T01:20:00.000Z 79.995588058317 -290.006145592032 11.453588310832 0.032364447630 -0.119732846726 0.004787795477 0.000006823058 -0.000032605399 0.000001002403 2000-01-01T01:30:00.000Z 100.653337431216 -367.973095343246 14.506689867546 0.036512309222 -0.140585791722 0.005389196595 0.000006837664 -0.000034976714 0.000001003336 2000-01-01T01:40:00.000Z 123.453665764425 -454.667651202580 17.923140739980 0.039513023305 -0.148793728322 0.005998979319 0.000005129687 -0.000015774366 0.000001016276 2000-01-01T01:50:00.000Z 148.100098098847 -547.040990119446 21.705450574838 0.042666722324 -0.159546931615 0.006608707351 0.000005379399 -0.000020069667 0.000001016148 2000-01-01T02:00:00.000Z 174.682530302316 -646.639412217862 25.853573367492 0.045964410184 -0.172877349174 0.007218354413 0.000005609651 -0.000024365071 0.000001016006 2000-01-01T02:10:00.000Z 203.283853139330 -755.009261136276 30.367457953699 0.049394412241 -0.188785025652 0.007827911891 0.000005820447 -0.000028660521 0.000001015850 2000-01-01T02:20:00.000Z 233.979953476934 -873.696897209172 35.247047991524 0.052945056828 -0.207269972054 0.008437371140 0.000006011793 -0.000032955961 0.000001015679 2000-01-01T02:30:00.000Z 266.839718231943 -1004.248702962090 40.492281925576 0.056604717930 -0.228332679264 0.009046723193 0.000006192961 -0.000037365335 0.000001015430 2000-01-01T02:40:00.000Z 301.553619978550 -1143.666361515510 46.105600609180 0.059130769542 -0.236774827689 0.009664348966 0.000004333782 -0.000016135290 0.000001029339 2000-01-01T02:50:00.000Z 337.827105187818 -1288.892787478942 52.089480854137 0.061805266666 -0.247741890570 0.010281901244 0.000004577970 -0.000020421606 0.000001029166 2000-01-01T03:00:00.000Z 375.748076996488 -1441.470994477480 58.443860505298 0.064620448977 -0.261280773723 0.010899345538 0.000004802733 -0.000024708015 0.000001028979 2000-01-01T03:10:00.000Z 415.397450190460 -1602.944087970172 65.168672162654 0.067564662535 -0.277391515874 0.011516673144 0.000005008075 -0.000028994460 0.000001028777 2000-01-01T03:20:00.000Z 456.849148063114 -1774.855186584263 72.263843192313 0.070626256317 -0.296074122196 0.012133875320 0.000005194002 -0.000033280886 0.000001028561 2000-01-01T03:20:00.000Z 457 -1775 72 0 -0 0 初始条件：{P(0.0, 0.0, 0.0), V(0.0, 0.0, 0.0), A(0.0, 0.0, 0.0)} 最终条件：{2000-01-01T03:20:00.000, P(456.8491480631136, -1774.8551865842628, 72.26384319231306), V(0.07062625631702495, -0.2960741221959937, 0.012133875319548914), A(5.194001612691349E-6, -3.3280886249059965E-5, 1.0285608966641585E-6)}

这个示例的完整代码可以在教程的源代码树中找到，文件名为 src/main/java/org/orekit/tutorials/propagation/PropagationInNonInertialFrame.java。