1   /* Copyright 2002-2021 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.estimation.measurements;
18  
19  import java.util.Arrays;
20  
21  import org.hipparchus.analysis.UnivariateFunction;
22  import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
23  import org.hipparchus.analysis.solvers.UnivariateSolver;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.util.FastMath;
26  import org.orekit.estimation.TLEContext;
27  import org.orekit.frames.Frame;
28  import org.orekit.frames.Transform;
29  import org.orekit.propagation.SpacecraftState;
30  import org.orekit.time.AbsoluteDate;
31  import org.orekit.utils.Constants;
32  import org.orekit.utils.PVCoordinates;
33  import org.orekit.utils.ParameterDriver;
34  
35  public class TLERangeRateMeasurementCreator extends MeasurementCreator {
36  
37      private final TLEContext context;
38      private final boolean twoWay;
39      private final ObservableSatellite satellite;
40  
41      public TLERangeRateMeasurementCreator(final TLEContext context, boolean twoWay) {
42          this.context   = context;
43          this.twoWay    = twoWay;
44          this.satellite = new ObservableSatellite(0);
45      }
46  
47      public void init(SpacecraftState s0, AbsoluteDate t, double step) {
48          for (final GroundStation station : context.stations) {
49              for (ParameterDriver driver : Arrays.asList(station.getClockOffsetDriver(),
50                                                          station.getNorthOffsetDriver(),
51                                                          station.getZenithOffsetDriver(),
52                                                          station.getPrimeMeridianOffsetDriver(),
53                                                          station.getPrimeMeridianDriftDriver(),
54                                                          station.getPolarOffsetXDriver(),
55                                                          station.getPolarDriftXDriver(),
56                                                          station.getPolarOffsetYDriver(),
57                                                          station.getPolarDriftYDriver())) {
58                  if (driver.getReferenceDate() == null) {
59                      driver.setReferenceDate(s0.getDate());
60                  }
61              }
62  
63          }
64      }
65  
66      public void handleStep(final SpacecraftState currentState) {
67          for (final GroundStation station : context.stations) {
68  
69              final AbsoluteDate     date      = currentState.getDate();
70              final Frame            inertial  = currentState.getFrame();
71              final Vector3D         position  = currentState.getPVCoordinates().getPosition();
72              final Vector3D         velocity  = currentState.getPVCoordinates().getVelocity();
73  
74              if (station.getBaseFrame().getElevation(position, inertial, date) > FastMath.toRadians(30.0)) {
75                  final UnivariateSolver solver = new BracketingNthOrderBrentSolver(1.0e-12, 5);
76  
77                  final double downLinkDelay  = solver.solve(1000, new UnivariateFunction() {
78                      public double value(final double x) {
79                          final Transform t = station.getOffsetToInertial(inertial, date.shiftedBy(x));
80                          final double d = Vector3D.distance(position, t.transformPosition(Vector3D.ZERO));
81                          return d - x * Constants.SPEED_OF_LIGHT;
82                      }
83                  }, -1.0, 1.0);
84                  final AbsoluteDate receptionDate  = currentState.getDate().shiftedBy(downLinkDelay);
85                  final PVCoordinates stationAtReception =
86                                  station.getOffsetToInertial(inertial, receptionDate).transformPVCoordinates(PVCoordinates.ZERO);
87  
88                  // line of sight at reception
89                  final Vector3D receptionLOS = (position.subtract(stationAtReception.getPosition())).normalize();
90  
91                  // relative velocity, spacecraft-station, at the date of reception
92                  final Vector3D deltaVr = velocity.subtract(stationAtReception.getVelocity());
93  
94                  final double upLinkDelay = solver.solve(1000, new UnivariateFunction() {
95                      public double value(final double x) {
96                          final Transform t = station.getOffsetToInertial(inertial, date.shiftedBy(-x));
97                          final double d = Vector3D.distance(position, t.transformPosition(Vector3D.ZERO));
98                          return d - x * Constants.SPEED_OF_LIGHT;
99                      }
100                 }, -1.0, 1.0);
101                 final AbsoluteDate emissionDate   = currentState.getDate().shiftedBy(-upLinkDelay);
102                 final PVCoordinates stationAtEmission  =
103                                 station.getOffsetToInertial(inertial, emissionDate).transformPVCoordinates(PVCoordinates.ZERO);
104 
105                 // line of sight at emission
106                 final Vector3D emissionLOS = (position.subtract(stationAtEmission.getPosition())).normalize();
107 
108                 // relative velocity, spacecraft-station, at the date of emission
109                 final Vector3D deltaVe = velocity.subtract(stationAtEmission.getVelocity());
110 
111                 // range rate at the date of reception
112                 final double rr = twoWay ?
113                                           0.5 * (deltaVr.dotProduct(receptionLOS) + deltaVe.dotProduct(emissionLOS)) :
114                                               deltaVr.dotProduct(receptionLOS);
115 
116                                           addMeasurement(new RangeRate(station, date,
117                                                                        rr,
118                                                                        1.0, 10, twoWay, satellite));
119             }
120 
121         }
122     }
123 
124 }