1   /* Copyright 2002-2022 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.attitudes;
18  
19  
20  import org.hipparchus.Field;
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.geometry.euclidean.threed.Rotation;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.util.Decimal64Field;
25  import org.hipparchus.util.FastMath;
26  import org.junit.Assert;
27  import org.junit.Before;
28  import org.junit.Test;
29  import org.orekit.Utils;
30  import org.orekit.bodies.CelestialBodyFactory;
31  import org.orekit.frames.Frame;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.orbits.FieldOrbit;
34  import org.orekit.orbits.KeplerianOrbit;
35  import org.orekit.orbits.Orbit;
36  import org.orekit.orbits.PositionAngle;
37  import org.orekit.propagation.FieldSpacecraftState;
38  import org.orekit.propagation.Propagator;
39  import org.orekit.propagation.SpacecraftState;
40  import org.orekit.propagation.analytical.KeplerianPropagator;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.DateComponents;
43  import org.orekit.time.FieldAbsoluteDate;
44  import org.orekit.time.TimeComponents;
45  import org.orekit.time.TimeScalesFactory;
46  import org.orekit.utils.AngularCoordinates;
47  import org.orekit.utils.PVCoordinates;
48  import org.orekit.utils.PVCoordinatesProvider;
49  
50  public class SpinStabilizedTest {
51  
52      @Test
53      public void testBBQMode() {
54          PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
55          AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 01, 01),
56                                               new TimeComponents(3, 25, 45.6789),
57                                               TimeScalesFactory.getTAI());
58          double rate = 2.0 * FastMath.PI / (12 * 60);
59          AttitudeProvider cbp = new CelestialBodyPointed(FramesFactory.getEME2000(), sun, Vector3D.PLUS_K,
60                                                          Vector3D.PLUS_I, Vector3D.PLUS_K);
61          SpinStabilized bbq = new SpinStabilized(cbp, date, Vector3D.PLUS_K, rate);
62          PVCoordinates pv =
63              new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
64                                new Vector3D(0, 0, 3680.853673522056));
65          KeplerianOrbit kep = new KeplerianOrbit(pv, FramesFactory.getEME2000(), date, 3.986004415e14);
66          Attitude attitude = bbq.getAttitude(kep, date, kep.getFrame());
67          checkField(Decimal64Field.getInstance(), bbq, kep, kep.getDate(), kep.getFrame());
68          Vector3D xDirection = attitude.getRotation().applyInverseTo(Vector3D.PLUS_I);
69          Assert.assertEquals(FastMath.atan(1.0 / 5000.0),
70                       Vector3D.angle(xDirection, sun.getPVCoordinates(date, FramesFactory.getEME2000()).getPosition()),
71                       2.0e-15);
72          Assert.assertEquals(rate, attitude.getSpin().getNorm(), 1.0e-6);
73          Assert.assertSame(cbp, bbq.getUnderlyingAttitudeProvider());
74  
75      }
76  
77      @Test
78      public void testSpin() {
79  
80          AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 01, 01),
81                                               new TimeComponents(3, 25, 45.6789),
82                                               TimeScalesFactory.getUTC());
83          double rate = 2.0 * FastMath.PI / (12 * 60);
84          AttitudeProvider law =
85              new SpinStabilized(new InertialProvider(Rotation.IDENTITY),
86                                 date, Vector3D.PLUS_K, rate);
87          KeplerianOrbit orbit =
88              new KeplerianOrbit(7178000.0, 1.e-4, FastMath.toRadians(50.),
89                                FastMath.toRadians(10.), FastMath.toRadians(20.),
90                                FastMath.toRadians(30.), PositionAngle.MEAN,
91                                FramesFactory.getEME2000(), date, 3.986004415e14);
92  
93          Propagator propagator = new KeplerianPropagator(orbit, law);
94  
95          double h = 10.0;
96          SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h));
97          SpacecraftState s0     = propagator.propagate(date);
98          SpacecraftState sPlus  = propagator.propagate(date.shiftedBy(h));
99          Vector3D spin0         = s0.getAttitude().getSpin();
100 
101         // check spin is consistent with attitude evolution
102         double errorAngleMinus     = Rotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(),
103                                                        s0.getAttitude().getRotation());
104         double evolutionAngleMinus = Rotation.distance(sMinus.getAttitude().getRotation(),
105                                                        s0.getAttitude().getRotation());
106         Assert.assertTrue(errorAngleMinus <= 1.0e-6 * evolutionAngleMinus);
107         double errorAnglePlus      = Rotation.distance(s0.getAttitude().getRotation(),
108                                                        sPlus.shiftedBy(-h).getAttitude().getRotation());
109         double evolutionAnglePlus  = Rotation.distance(s0.getAttitude().getRotation(),
110                                                        sPlus.getAttitude().getRotation());
111         Assert.assertTrue(errorAnglePlus <= 1.0e-6 * evolutionAnglePlus);
112 
113         // compute spin axis using finite differences
114         Rotation rM = sMinus.getAttitude().getRotation();
115         Rotation rP = sPlus.getAttitude().getRotation();
116         Vector3D reference = AngularCoordinates.estimateRate(rM, rP, 2 * h);
117 
118         Assert.assertEquals(2 * FastMath.PI / reference.getNorm(), 2 * FastMath.PI / spin0.getNorm(), 0.05);
119         Assert.assertEquals(0.0, FastMath.toDegrees(Vector3D.angle(reference, spin0)), 1.0e-10);
120         Assert.assertEquals(0.0, FastMath.toDegrees(Vector3D.angle(Vector3D.PLUS_K, spin0)), 1.0e-10);
121 
122     }
123 
124     private <T extends CalculusFieldElement<T>> void checkField(final Field<T> field, final AttitudeProvider provider,
125                                                             final Orbit orbit, final AbsoluteDate date,
126                                                             final Frame frame)
127         {
128         Attitude attitudeD = provider.getAttitude(orbit, date, frame);
129         final FieldOrbit<T> orbitF = new FieldSpacecraftState<>(field, new SpacecraftState(orbit)).getOrbit();
130         final FieldAbsoluteDate<T> dateF = new FieldAbsoluteDate<>(field, date);
131         FieldAttitude<T> attitudeF = provider.getAttitude(orbitF, dateF, frame);
132         Assert.assertEquals(0.0, Rotation.distance(attitudeD.getRotation(), attitudeF.getRotation().toRotation()), 1.0e-15);
133         Assert.assertEquals(0.0, Vector3D.distance(attitudeD.getSpin(), attitudeF.getSpin().toVector3D()), 1.0e-15);
134         Assert.assertEquals(0.0, Vector3D.distance(attitudeD.getRotationAcceleration(), attitudeF.getRotationAcceleration().toVector3D()), 1.0e-15);
135     }
136 
137     @Before
138     public void setUp() {
139         Utils.setDataRoot("regular-data");
140     }
141 
142 }
143