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.gnss.attitude;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.RealFieldElement;
21  import org.hipparchus.util.FastMath;
22  import org.orekit.frames.Frame;
23  import org.orekit.time.AbsoluteDate;
24  import org.orekit.utils.ExtendedPVCoordinatesProvider;
25  import org.orekit.utils.TimeStampedAngularCoordinates;
26  import org.orekit.utils.TimeStampedFieldAngularCoordinates;
27  
28  /**
29   * Attitude providers for GPS block IIF navigation satellites.
30   * <p>
31   * This class is based on the May 2017 version of J. Kouba eclips.f
32   * subroutine available at <a href="http://acc.igs.org/orbits">IGS Analysis
33   * Center Coordinator site</a>. The eclips.f code itself is not used ; its
34   * hard-coded data are used and its low level models are used, but the
35   * structure of the code and the API have been completely rewritten.
36   * </p>
37   * @author J. Kouba original fortran routine
38   * @author Luc Maisonobe Java translation
39   * @since 9.2
40   */
41  public class GPSBlockIIF extends AbstractGNSSAttitudeProvider {
42  
43      /** Default yaw rates for all spacecrafts in radians per seconds. */
44      public static final double DEFAULT_YAW_RATE = FastMath.toRadians(0.11);
45  
46      /** Default yaw bias (rad). */
47      public static final double DEFAULT_YAW_BIAS = FastMath.toRadians(-0.7);
48  
49      /** Serializable UID. */
50      private static final long serialVersionUID = 20171114L;
51  
52      /** Satellite-Sun angle limit for a midnight turn maneuver. */
53      private static final double NIGHT_TURN_LIMIT = FastMath.toRadians(180.0 - 13.25);
54  
55      /** Margin on turn end. */
56      private final double END_MARGIN = 1800.0;
57  
58      /** Yaw rate. */
59      private final double yawRate;
60  
61      /** Yaw bias. */
62      private final double yawBias;
63  
64      /** Simple constructor.
65       * @param yawRate yaw rate to use in radians per seconds (typically {@link #DEFAULT_YAW_RATE})
66       * @param yawBias yaw bias to use (rad) (typicall {@link #DEFAULT_YAW_BIAS})
67       * @param validityStart start of validity for this provider
68       * @param validityEnd end of validity for this provider
69       * @param sun provider for Sun position
70       * @param inertialFrame inertial frame where velocity are computed
71       */
72      public GPSBlockIIF(final double yawRate, final double yawBias,
73                         final AbsoluteDateDate">AbsoluteDate validityStart, final AbsoluteDate validityEnd,
74                         final ExtendedPVCoordinatesProvider sun, final Frame inertialFrame) {
75          super(validityStart, validityEnd, sun, inertialFrame);
76          this.yawRate = yawRate;
77          this.yawBias = yawBias;
78      }
79  
80      /** {@inheritDoc} */
81      @Override
82      protected TimeStampedAngularCoordinates correctedYaw(final GNSSAttitudeContext context) {
83  
84          // noon beta angle limit from yaw rate
85          final double aNoon  = FastMath.atan(context.getMuRate() / yawRate);
86          final double aNight = NIGHT_TURN_LIMIT;
87          final double cNoon  = FastMath.cos(aNoon);
88          final double cNight = FastMath.cos(aNight);
89  
90          if (context.setUpTurnRegion(cNight, cNoon)) {
91  
92              final double absBeta = FastMath.abs(context.beta(context.getDate()));
93              context.setHalfSpan(context.inSunSide() ?
94                                  absBeta * FastMath.sqrt(aNoon / absBeta - 1.0) :
95                                  context.inOrbitPlaneAbsoluteAngle(aNight - FastMath.PI),
96                                  END_MARGIN);
97              if (context.inTurnTimeRange()) {
98  
99                  // we need to ensure beta sign does not change during the turn
100                 final double beta     = context.getSecuredBeta();
101                 final double phiStart = context.getYawStart(beta);
102                 final double dtStart  = context.timeSinceTurnStart();
103                 final double phiDot;
104                 final double linearPhi;
105                 if (context.inSunSide()) {
106                     // noon turn
107                     if (beta > yawBias && beta < 0) {
108                         // noon turn problem for small negative beta in block IIF
109                         // rotation is in the wrong direction for these spacecrafts
110                         phiDot    = FastMath.copySign(yawRate, beta);
111                         linearPhi = phiStart + phiDot * dtStart;
112                     } else {
113                         // regular noon turn
114                         phiDot    = -FastMath.copySign(yawRate, beta);
115                         linearPhi = phiStart + phiDot * dtStart;
116                     }
117                 } else {
118                     // midnight turn
119                     phiDot    = context.yawRate(beta);
120                     linearPhi = phiStart + phiDot * dtStart;
121                 }
122 
123                 if (context.linearModelStillActive(linearPhi, phiDot)) {
124                     // we are still in the linear model phase
125                     return context.turnCorrectedAttitude(linearPhi, phiDot);
126                 }
127 
128 
129             }
130 
131         }
132 
133         // in nominal yaw mode
134         return context.nominalYaw(context.getDate());
135 
136     }
137 
138     /** {@inheritDoc} */
139     @Override
140     protected <T extends RealFieldElement<T>> TimeStampedFieldAngularCoordinates<T> correctedYaw(final GNSSFieldAttitudeContext<T> context) {
141 
142         final Field<T> field = context.getDate().getField();
143 
144         // noon beta angle limit from yaw rate
145         final T      aNoon  = FastMath.atan(context.getMuRate().divide(yawRate));
146         final T      aNight = field.getZero().add(NIGHT_TURN_LIMIT);
147         final double cNoon  = FastMath.cos(aNoon.getReal());
148         final double cNight = FastMath.cos(aNight.getReal());
149 
150         if (context.setUpTurnRegion(cNight, cNoon)) {
151 
152             final T absBeta = FastMath.abs(context.beta(context.getDate()));
153             context.setHalfSpan(context.inSunSide() ?
154                                 absBeta.multiply(FastMath.sqrt(aNoon.divide(absBeta).subtract(1.0))) :
155                                 context.inOrbitPlaneAbsoluteAngle(aNight.subtract(FastMath.PI)),
156                                 END_MARGIN);
157             if (context.inTurnTimeRange()) {
158 
159                 // we need to ensure beta sign does not change during the turn
160                 final T beta     = context.getSecuredBeta();
161                 final T phiStart = context.getYawStart(beta);
162                 final T dtStart  = context.timeSinceTurnStart();
163                 final T phiDot;
164                 final T linearPhi;
165                 if (context.inSunSide()) {
166                     // noon turn
167                     if (beta.getReal() > yawBias && beta.getReal() < 0) {
168                         // noon turn problem for small negative beta in block IIF
169                         // rotation is in the wrong direction for these spacecrafts
170                         phiDot    = field.getZero().add(FastMath.copySign(yawRate, beta.getReal()));
171                         linearPhi = phiStart.add(phiDot.multiply(dtStart));
172                     } else {
173                         // regular noon turn
174                         phiDot    = field.getZero().add(-FastMath.copySign(yawRate, beta.getReal()));
175                         linearPhi = phiStart.add(phiDot.multiply(dtStart));
176                     }
177                 } else {
178                     // midnight turn
179                     phiDot    = context.yawRate(beta);
180                     linearPhi = phiStart.add(phiDot.multiply(dtStart));
181                 }
182 
183                 if (context.linearModelStillActive(linearPhi, phiDot)) {
184                     // we are still in the linear model phase
185                     return context.turnCorrectedAttitude(linearPhi, phiDot);
186                 }
187 
188             }
189 
190         }
191 
192         // in nominal yaw mode
193         return context.nominalYaw(context.getDate());
194 
195     }
196 
197 }