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 IIR 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 GPSBlockIIA extends AbstractGNSSAttitudeProvider {
42  
43      /** Default yaw rates for all spacecrafts in radians per seconds (indexed by prnNumber, hence entry 0 is unused). */
44      public static final double[] DEFAULT_YAW_RATES = new double[] {
45          Double.NaN,  // unused entry 0
46          FastMath.toRadians(0.1211), FastMath.toRadians(0.1339), FastMath.toRadians(0.1230), FastMath.toRadians(0.1233),
47          FastMath.toRadians(0.1180), FastMath.toRadians(0.1266), FastMath.toRadians(0.1269), FastMath.toRadians(0.1033),
48          FastMath.toRadians(0.1278), FastMath.toRadians(0.0978), FastMath.toRadians(0.2000), FastMath.toRadians(0.1990),
49          FastMath.toRadians(0.2000), FastMath.toRadians(0.0815), FastMath.toRadians(0.1303), FastMath.toRadians(0.0838),
50          FastMath.toRadians(0.1401), FastMath.toRadians(0.1069), FastMath.toRadians(0.0980), FastMath.toRadians(0.1030),
51          FastMath.toRadians(0.1366), FastMath.toRadians(0.1025), FastMath.toRadians(0.1140), FastMath.toRadians(0.1089),
52          FastMath.toRadians(0.1001), FastMath.toRadians(0.1227), FastMath.toRadians(0.1194), FastMath.toRadians(0.1260),
53          FastMath.toRadians(0.1228), FastMath.toRadians(0.1165), FastMath.toRadians(0.0969), FastMath.toRadians(0.1140)
54      };
55  
56      /** Default yaw bias (rad). */
57      public static final double DEFAULT_YAW_BIAS = FastMath.toRadians(0.5);
58  
59      /** Serializable UID. */
60      private static final long serialVersionUID = 20171114L;
61  
62      /** Satellite-Sun angle limit for a midnight turn maneuver. */
63      private static final double NIGHT_TURN_LIMIT = FastMath.toRadians(180.0 - 13.25);
64  
65      /** Margin on turn end. */
66      private final double END_MARGIN = 1800.0;
67  
68      /** Yaw rate for current spacecraft. */
69      private final double yawRate;
70  
71      /** Yaw bias. */
72      private final double yawBias;
73  
74      /** Simple constructor.
75       * @param yawRate yaw rate to use in radians per seconds (typically {@link #DEFAULT_YAW_RATES}{@code [prnNumber]})
76       * @param yawBias yaw bias to use (rad) (typicall {@link #DEFAULT_YAW_BIAS})
77       * @param validityStart start of validity for this provider
78       * @param validityEnd end of validity for this provider
79       * @param sun provider for Sun position
80       * @param inertialFrame inertial frame where velocity are computed
81       * @param prnNumber number within the GPS constellation (between 1 and 32)
82       * @deprecated as of 9.3 replaced by {@link #GPSBlockIIA(double, double,
83       * AbsoluteDate, AbsoluteDate, ExtendedPVCoordinatesProvider, Frame)}
84       */
85      @Deprecated
86      public GPSBlockIIA(final double yawRate, final double yawBias,
87                         final AbsoluteDateDate">AbsoluteDate validityStart, final AbsoluteDate validityEnd,
88                         final ExtendedPVCoordinatesProvider sun, final Frame inertialFrame, final int prnNumber) {
89          this(yawRate, yawBias, validityStart, validityEnd, sun, inertialFrame);
90      }
91  
92      /** Simple constructor.
93       * @param yawRate yaw rate to use in radians per seconds (typically {@link #DEFAULT_YAW_RATES}{@code [prnNumber]})
94       * @param yawBias yaw bias to use (rad) (typicall {@link #DEFAULT_YAW_BIAS})
95       * @param validityStart start of validity for this provider
96       * @param validityEnd end of validity for this provider
97       * @param sun provider for Sun position
98       * @param inertialFrame inertial frame where velocity are computed
99       * @since 9.3
100      */
101     public GPSBlockIIA(final double yawRate, final double yawBias,
102                        final AbsoluteDateDate">AbsoluteDate validityStart, final AbsoluteDate validityEnd,
103                        final ExtendedPVCoordinatesProvider sun, final Frame inertialFrame) {
104         super(validityStart, validityEnd, sun, inertialFrame);
105         this.yawRate = yawRate;
106         this.yawBias = yawBias;
107     }
108 
109     /** {@inheritDoc} */
110     @Override
111     protected TimeStampedAngularCoordinates correctedYaw(final GNSSAttitudeContext context) {
112 
113         // noon beta angle limit from yaw rate
114         final double aNoon  = FastMath.atan(context.getMuRate() / yawRate);
115         final double aNight = NIGHT_TURN_LIMIT;
116         final double cNoon  = FastMath.cos(aNoon);
117         final double cNight = FastMath.cos(aNight);
118 
119         if (context.setUpTurnRegion(cNight, cNoon)) {
120 
121             final double absBeta = FastMath.abs(context.beta(context.getDate()));
122             context.setHalfSpan(context.inSunSide() ?
123                                 absBeta * FastMath.sqrt(aNoon / absBeta - 1.0) :
124                                 context.inOrbitPlaneAbsoluteAngle(aNight - FastMath.PI),
125                                 END_MARGIN);
126             if (context.inTurnTimeRange()) {
127 
128                 // we need to ensure beta sign does not change during the turn
129                 final double beta     = context.getSecuredBeta();
130                 final double phiStart = context.getYawStart(beta);
131                 final double dtStart  = context.timeSinceTurnStart();
132                 final double linearPhi;
133                 final double phiDot;
134                 if (context.inSunSide()) {
135                     // noon turn
136                     if (beta > 0 && beta < yawBias) {
137                         // noon turn problem for small positive beta in block IIA
138                         // rotation is in the wrong direction for these spacecrafts
139                         phiDot    = FastMath.copySign(yawRate, beta);
140                         linearPhi = phiStart + phiDot * dtStart;
141                     } else {
142                         // regular noon turn
143                         phiDot    = -FastMath.copySign(yawRate, beta);
144                         linearPhi = phiStart + phiDot * dtStart;
145                     }
146                 } else {
147                     // midnight turn
148                     phiDot    = yawRate;
149                     linearPhi = phiStart + phiDot * dtStart;
150                 }
151 
152                 if (context.linearModelStillActive(linearPhi, phiDot)) {
153                     // we are still in the linear model phase
154                     return context.turnCorrectedAttitude(linearPhi, phiDot);
155                 }
156 
157             }
158 
159         }
160 
161         // in nominal yaw mode
162         return context.nominalYaw(context.getDate());
163 
164     }
165 
166     /** {@inheritDoc} */
167     @Override
168     protected <T extends RealFieldElement<T>> TimeStampedFieldAngularCoordinates<T> correctedYaw(final GNSSFieldAttitudeContext<T> context) {
169 
170         final Field<T> field = context.getDate().getField();
171 
172         // noon beta angle limit from yaw rate
173         final T      aNoon  = FastMath.atan(context.getMuRate().divide(yawRate));
174         final T      aNight = field.getZero().add(NIGHT_TURN_LIMIT);
175         final double cNoon  = FastMath.cos(aNoon.getReal());
176         final double cNight = FastMath.cos(aNight.getReal());
177 
178         if (context.setUpTurnRegion(cNight, cNoon)) {
179 
180             final T absBeta = FastMath.abs(context.beta(context.getDate()));
181             context.setHalfSpan(context.inSunSide() ?
182                                 absBeta.multiply(FastMath.sqrt(aNoon.divide(absBeta).subtract(1.0))) :
183                                 context.inOrbitPlaneAbsoluteAngle(aNight.subtract(FastMath.PI)),
184                                 END_MARGIN);
185             if (context.inTurnTimeRange()) {
186 
187                 // we need to ensure beta sign does not change during the turn
188                 final T beta     = context.getSecuredBeta();
189                 final T phiStart = context.getYawStart(beta);
190                 final T dtStart  = context.timeSinceTurnStart();
191                 final T linearPhi;
192                 final T phiDot;
193                 if (context.inSunSide()) {
194                     // noon turn
195                     if (beta.getReal() > 0 && beta.getReal() < yawBias) {
196                         // noon turn problem for small positive beta in block IIA
197                         // rotation is in the wrong direction for these spacecrafts
198                         phiDot    = field.getZero().add(FastMath.copySign(yawRate, beta.getReal()));
199                         linearPhi = phiStart.add(phiDot.multiply(dtStart));
200                     } else {
201                         // regular noon turn
202                         phiDot    = field.getZero().add(-FastMath.copySign(yawRate, beta.getReal()));
203                         linearPhi = phiStart.add(phiDot.multiply(dtStart));
204                     }
205                 } else {
206                     // midnight turn
207                     phiDot    = field.getZero().add(yawRate);
208                     linearPhi = phiStart.add(phiDot.multiply(dtStart));
209                 }
210 
211                 if (context.linearModelStillActive(linearPhi, phiDot)) {
212                     // we are still in the linear model phase
213                     return context.turnCorrectedAttitude(linearPhi, phiDot);
214                 }
215 
216             }
217 
218         }
219 
220         // in nominal yaw mode
221         return context.nominalYaw(context.getDate());
222 
223     }
224 
225 }