1   /* Copyright 2002-2019 CS Systèmes d'Information
2    * Licensed to CS Systèmes d'Information (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.forces;
18  
19  
20  import org.hipparchus.analysis.differentiation.DerivativeStructure;
21  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
22  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.util.Precision;
25  import org.junit.Assert;
26  import org.orekit.attitudes.AttitudeProvider;
27  import org.orekit.frames.Frame;
28  import org.orekit.propagation.FieldSpacecraftState;
29  import org.orekit.propagation.SpacecraftState;
30  import org.orekit.time.AbsoluteDate;
31  
32  
33  public abstract class AbstractLegacyForceModelTest extends AbstractForceModelTest {
34  
35      protected abstract FieldVector3D<DerivativeStructure> accelerationDerivatives(ForceModel forceModel,
36                                                                                    final AbsoluteDate date, final  Frame frame,
37                                                                                    final FieldVector3D<DerivativeStructure> position,
38                                                                                    final FieldVector3D<DerivativeStructure> velocity,
39                                                                                    final FieldRotation<DerivativeStructure> rotation,
40                                                                                    final DerivativeStructure mass);
41  
42      protected void checkStateJacobianVs80Implementation(final SpacecraftState state, final ForceModel forceModel,
43                                                          final AttitudeProvider attitudeProvider,
44                                                          final double checkTolerance, final boolean print)
45          {
46          FieldSpacecraftState<DerivativeStructure> fState = toDS(state, attitudeProvider);
47          FieldVector3D<DerivativeStructure> dsNew = forceModel.acceleration(fState,
48                                                                             forceModel.getParameters(fState.getDate().getField()));
49          FieldVector3D<DerivativeStructure> dsOld = accelerationDerivatives(forceModel, fState.getDate().toAbsoluteDate(),
50                                                                             fState.getFrame(),
51                                                                             fState.getPVCoordinates().getPosition(),
52                                                                             fState.getPVCoordinates().getVelocity(),
53                                                                             fState.getAttitude().getRotation(),
54                                                                             fState.getMass());
55          Vector3D dFdPXRef = new Vector3D(dsOld.getX().getPartialDerivative(1, 0, 0, 0, 0, 0),
56                                           dsOld.getY().getPartialDerivative(1, 0, 0, 0, 0, 0),
57                                           dsOld.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
58          Vector3D dFdPXRes = new Vector3D(dsNew.getX().getPartialDerivative(1, 0, 0, 0, 0, 0),
59                                           dsNew.getY().getPartialDerivative(1, 0, 0, 0, 0, 0),
60                                           dsNew.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
61          Vector3D dFdPYRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 1, 0, 0, 0, 0),
62                                           dsOld.getY().getPartialDerivative(0, 1, 0, 0, 0, 0),
63                                           dsOld.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
64          Vector3D dFdPYRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 1, 0, 0, 0, 0),
65                                           dsNew.getY().getPartialDerivative(0, 1, 0, 0, 0, 0),
66                                           dsNew.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
67          Vector3D dFdPZRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 1, 0, 0, 0),
68                                           dsOld.getY().getPartialDerivative(0, 0, 1, 0, 0, 0),
69                                           dsOld.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
70          Vector3D dFdPZRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 1, 0, 0, 0),
71                                           dsNew.getY().getPartialDerivative(0, 0, 1, 0, 0, 0),
72                                           dsNew.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
73          Vector3D dFdVXRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 1, 0, 0),
74                                           dsOld.getY().getPartialDerivative(0, 0, 0, 1, 0, 0),
75                                           dsOld.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
76          Vector3D dFdVXRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 1, 0, 0),
77                                           dsNew.getY().getPartialDerivative(0, 0, 0, 1, 0, 0),
78                                           dsNew.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
79          Vector3D dFdVYRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 0, 1, 0),
80                                           dsOld.getY().getPartialDerivative(0, 0, 0, 0, 1, 0),
81                                           dsOld.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
82          Vector3D dFdVYRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 0, 1, 0),
83                                           dsNew.getY().getPartialDerivative(0, 0, 0, 0, 1, 0),
84                                           dsNew.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
85          Vector3D dFdVZRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 0, 0, 1),
86                                           dsOld.getY().getPartialDerivative(0, 0, 0, 0, 0, 1),
87                                           dsOld.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
88          Vector3D dFdVZRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 0, 0, 1),
89                                           dsNew.getY().getPartialDerivative(0, 0, 0, 0, 0, 1),
90                                           dsNew.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
91          if (print) {
92              System.out.println("dF/dPX ref norm: " + dFdPXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPXRef, dFdPXRes) + ", rel error: " + (Vector3D.distance(dFdPXRef, dFdPXRes) / dFdPXRef.getNorm()));
93              System.out.println("dF/dPY ref norm: " + dFdPYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPYRef, dFdPYRes) + ", rel error: " + (Vector3D.distance(dFdPYRef, dFdPYRes) / dFdPYRef.getNorm()));
94              System.out.println("dF/dPZ ref norm: " + dFdPZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPZRef, dFdPZRes) + ", rel error: " + (Vector3D.distance(dFdPZRef, dFdPZRes) / dFdPZRef.getNorm()));
95              System.out.println("dF/dVX ref norm: " + dFdVXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVXRef, dFdVXRes) + ", rel error: " + (Vector3D.distance(dFdVXRef, dFdVXRes) / dFdVXRef.getNorm()));
96              System.out.println("dF/dVY ref norm: " + dFdVYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVYRef, dFdVYRes) + ", rel error: " + (Vector3D.distance(dFdVYRef, dFdVYRes) / dFdVYRef.getNorm()));
97              System.out.println("dF/dVZ ref norm: " + dFdVZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVZRef, dFdVZRes) + ", rel error: " + (Vector3D.distance(dFdVZRef, dFdVZRes) / dFdVZRef.getNorm()));
98          }
99          checkdFdP(dFdPXRef, dFdPXRes, checkTolerance);
100         checkdFdP(dFdPYRef, dFdPYRes, checkTolerance);
101         checkdFdP(dFdPZRef, dFdPZRes, checkTolerance);
102         checkdFdP(dFdVXRef, dFdVXRes, checkTolerance);
103         checkdFdP(dFdVYRef, dFdVYRes, checkTolerance);
104         checkdFdP(dFdVZRef, dFdVZRes, checkTolerance);
105 
106     }
107 
108     private void checkdFdP(final Vector3D reference, final Vector3D result, final double checkTolerance) {
109         if (reference.getNorm() == 0) {
110             // if dF/dP is exactly zero (i.e. no dependency between F and P),
111             // then the result should also be exactly zero
112             Assert.assertEquals(0, result.getNorm(), Precision.SAFE_MIN);
113         } else {
114             Assert.assertEquals(0, Vector3D.distance(reference, result), checkTolerance * reference.getNorm());
115         }
116     }
117 
118 }
119 
120