Search in sources :

Example 1 with OrekitExceptionWrapper

use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.

the class TurnAroundRangeMeasurementCreator method handleStep.

/**
 * Function handling the steps of the propagator
 * A turn-around measurement needs 2 stations, a master and a slave
 * The measurement is a signal:
 * - Emitted from the master ground station
 * - Reflected on the spacecraft
 * - Reflected on the slave ground station
 * - Reflected on the spacecraft again
 * - Received on the master ground station
 * Its value is the elapsed time between emission and reception
 * divided by 2c were c is the speed of light.
 *
 * The path of the signal is divided into 2 legs:
 *  - The 1st leg goes from emission by the master station to reception by the slave station
 *  - The 2nd leg goes from emission by the slave station to reception by the master station
 *
 * The spacecraft state date should, after a few iterations of the estimation process, be
 * set to the date of arrival/departure of the signal to/from the slave station.
 * It is guaranteed by implementation of the estimated measurement.
 * This is done to avoid big shifts in time to compute the transit states.
 * See TurnAroundRange.java for more
 * Thus the spacecraft date is the date when the 1st leg of the path ends and the 2nd leg begins
 */
public void handleStep(final SpacecraftState currentState, final boolean isLast) throws OrekitException {
    try {
        for (Map.Entry<GroundStation, GroundStation> entry : context.TARstations.entrySet()) {
            final GroundStation masterStation = entry.getKey();
            final GroundStation slaveStation = entry.getValue();
            final AbsoluteDate date = currentState.getDate();
            final Frame inertial = currentState.getFrame();
            final Vector3D position = currentState.toTransform().getInverse().transformPosition(antennaPhaseCenter);
            // Create a TAR measurement only if elevation for both stations is higher than elevationMin°
            if ((masterStation.getBaseFrame().getElevation(position, inertial, date) > FastMath.toRadians(30.0)) && (slaveStation.getBaseFrame().getElevation(position, inertial, date) > FastMath.toRadians(30.0))) {
                // The solver used
                final UnivariateSolver solver = new BracketingNthOrderBrentSolver(1.0e-12, 5);
                // Spacecraft date t = date of arrival/departure of the signal to/from from the slave station
                // Slave station position in inertial frame at t
                final Vector3D slaveStationPosition = slaveStation.getOffsetToInertial(inertial, date).transformPosition(Vector3D.ZERO);
                // Downlink time of flight to slave station
                // The date of arrival/departure of the signal to/from the slave station is known and
                // equal to spacecraft date t.
                // Therefore we can use the function "downlinkTimeOfFlight" from GroundStation class
                // final double slaveTauD = slaveStation.downlinkTimeOfFlight(currentState, date);
                final double slaveTauD = solver.solve(1000, new UnivariateFunction() {

                    public double value(final double x) throws OrekitExceptionWrapper {
                        final SpacecraftState transitState = currentState.shiftedBy(-x);
                        final double d = Vector3D.distance(transitState.toTransform().getInverse().transformPosition(antennaPhaseCenter), slaveStationPosition);
                        return d - x * Constants.SPEED_OF_LIGHT;
                    }
                }, -1.0, 1.0);
                // Uplink time of flight from slave station
                // A solver is used to know where the satellite is when it receives the signal
                // back from the slave station
                final double slaveTauU = solver.solve(1000, new UnivariateFunction() {

                    public double value(final double x) throws OrekitExceptionWrapper {
                        final SpacecraftState transitState = currentState.shiftedBy(+x);
                        final double d = Vector3D.distance(transitState.toTransform().getInverse().transformPosition(antennaPhaseCenter), slaveStationPosition);
                        return d - x * Constants.SPEED_OF_LIGHT;
                    }
                }, -1.0, 1.0);
                // Find the position of the master station at signal departure and arrival
                // ----
                // Transit state position & date for the 1st leg of the signal path
                final SpacecraftState S1 = currentState.shiftedBy(-slaveTauD);
                final Vector3D P1 = S1.toTransform().getInverse().transformPosition(antennaPhaseCenter);
                final AbsoluteDate T1 = date.shiftedBy(-slaveTauD);
                // Transit state position & date for the 2nd leg of the signal path
                final Vector3D P2 = currentState.shiftedBy(+slaveTauU).toTransform().getInverse().transformPosition(antennaPhaseCenter);
                final AbsoluteDate T2 = date.shiftedBy(+slaveTauU);
                // Master station downlink delay - from P2 to master station
                // We use a solver to know where the master station is when it receives
                // the signal back from the satellite on the 2nd leg of the path
                final double masterTauD = solver.solve(1000, new UnivariateFunction() {

                    public double value(final double x) throws OrekitExceptionWrapper {
                        try {
                            final Transform t = masterStation.getOffsetToInertial(inertial, T2.shiftedBy(+x));
                            final double d = Vector3D.distance(P2, t.transformPosition(Vector3D.ZERO));
                            return d - x * Constants.SPEED_OF_LIGHT;
                        } catch (OrekitException oe) {
                            throw new OrekitExceptionWrapper(oe);
                        }
                    }
                }, -1.0, 1.0);
                final AbsoluteDate masterReceptionDate = T2.shiftedBy(+masterTauD);
                final TimeStampedPVCoordinates masterStationAtReception = masterStation.getOffsetToInertial(inertial, masterReceptionDate).transformPVCoordinates(new TimeStampedPVCoordinates(masterReceptionDate, PVCoordinates.ZERO));
                // Master station uplink delay - from master station to P1
                // Here the state date is known. Thus we can use the function "signalTimeOfFlight"
                // of the AbstractMeasurement class
                final double masterTauU = AbstractMeasurement.signalTimeOfFlight(masterStationAtReception, P1, T1);
                final AbsoluteDate masterEmissionDate = T1.shiftedBy(-masterTauU);
                final Vector3D masterStationAtEmission = masterStation.getOffsetToInertial(inertial, masterEmissionDate).transformPosition(Vector3D.ZERO);
                // Uplink/downlink distance from/to slave station
                final double slaveDownLinkDistance = Vector3D.distance(P1, slaveStationPosition);
                final double slaveUpLinkDistance = Vector3D.distance(P2, slaveStationPosition);
                // Uplink/downlink distance from/to master station
                final double masterUpLinkDistance = Vector3D.distance(P1, masterStationAtEmission);
                final double masterDownLinkDistance = Vector3D.distance(P2, masterStationAtReception.getPosition());
                addMeasurement(new TurnAroundRange(masterStation, slaveStation, masterReceptionDate, 0.5 * (masterUpLinkDistance + slaveDownLinkDistance + slaveUpLinkDistance + masterDownLinkDistance), 1.0, 10));
            }
        }
    } catch (OrekitExceptionWrapper oew) {
        throw new OrekitException(oew.getException());
    } catch (OrekitException oe) {
        throw new OrekitException(oe);
    }
}
Also used : Frame(org.orekit.frames.Frame) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) UnivariateSolver(org.hipparchus.analysis.solvers.UnivariateSolver) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) AbsoluteDate(org.orekit.time.AbsoluteDate) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) OrekitException(org.orekit.errors.OrekitException) Transform(org.orekit.frames.Transform) Map(java.util.Map)

Example 2 with OrekitExceptionWrapper

use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.

the class SpacecraftState method interpolate.

/**
 * {@inheritDoc}
 * <p>
 * The additional states that are interpolated are the ones already present
 * in the instance. The sample instances must therefore have at least the same
 * additional states has the instance. They may have more additional states,
 * but the extra ones will be ignored.
 * </p>
 * <p>
 * As this implementation of interpolation is polynomial, it should be used only
 * with small samples (about 10-20 points) in order to avoid <a
 * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
 * and numerical problems (including NaN appearing).
 * </p>
 */
public SpacecraftState interpolate(final AbsoluteDate date, final Stream<SpacecraftState> sample) throws OrekitException {
    // prepare interpolators
    final List<Orbit> orbits = new ArrayList<>();
    final List<Attitude> attitudes = new ArrayList<>();
    final HermiteInterpolator massInterpolator = new HermiteInterpolator();
    final Map<String, HermiteInterpolator> additionalInterpolators = new HashMap<String, HermiteInterpolator>(additional.size());
    for (final String name : additional.keySet()) {
        additionalInterpolators.put(name, new HermiteInterpolator());
    }
    // extract sample data
    try {
        sample.forEach(state -> {
            try {
                final double deltaT = state.getDate().durationFrom(date);
                orbits.add(state.getOrbit());
                attitudes.add(state.getAttitude());
                massInterpolator.addSamplePoint(deltaT, new double[] { state.getMass() });
                for (final Map.Entry<String, HermiteInterpolator> entry : additionalInterpolators.entrySet()) {
                    entry.getValue().addSamplePoint(deltaT, state.getAdditionalState(entry.getKey()));
                }
            } catch (OrekitException oe) {
                throw new OrekitExceptionWrapper(oe);
            }
        });
    } catch (OrekitExceptionWrapper oew) {
        throw oew.getException();
    }
    // perform interpolations
    final Orbit interpolatedOrbit = orbit.interpolate(date, orbits);
    final Attitude interpolatedAttitude = attitude.interpolate(date, attitudes);
    final double interpolatedMass = massInterpolator.value(0)[0];
    final Map<String, double[]> interpolatedAdditional;
    if (additional.isEmpty()) {
        interpolatedAdditional = null;
    } else {
        interpolatedAdditional = new HashMap<String, double[]>(additional.size());
        for (final Map.Entry<String, HermiteInterpolator> entry : additionalInterpolators.entrySet()) {
            interpolatedAdditional.put(entry.getKey(), entry.getValue().value(0));
        }
    }
    // create the complete interpolated state
    return new SpacecraftState(interpolatedOrbit, interpolatedAttitude, interpolatedMass, interpolatedAdditional);
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) Orbit(org.orekit.orbits.Orbit) Attitude(org.orekit.attitudes.Attitude) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) HermiteInterpolator(org.hipparchus.analysis.interpolation.HermiteInterpolator) OrekitException(org.orekit.errors.OrekitException) HashMap(java.util.HashMap) Map(java.util.Map)

Example 3 with OrekitExceptionWrapper

use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.

the class FieldEventState method findRoot.

/**
 * Find a root in a bracketing interval.
 *
 * <p> When calling this method one of the following must be true. Either ga == 0, gb
 * == 0, (ga < 0  and gb > 0), or (ga > 0 and gb < 0).
 *
 * @param interpolator that covers the interval.
 * @param ta           earliest possible time for root.
 * @param ga           g(ta).
 * @param tb           latest possible time for root.
 * @param gb           g(tb).
 * @return if a zero crossing was found.
 * @throws OrekitException if the event detector throws one
 */
private boolean findRoot(final FieldOrekitStepInterpolator<T> interpolator, final FieldAbsoluteDate<T> ta, final T ga, final FieldAbsoluteDate<T> tb, final T gb) throws OrekitException {
    final T zero = ga.getField().getZero();
    // check there appears to be a root in [ta, tb]
    check(ga.getReal() == 0.0 || gb.getReal() == 0.0 || (ga.getReal() > 0.0 && gb.getReal() < 0.0) || (ga.getReal() < 0.0 && gb.getReal() > 0.0));
    final T convergence = detector.getThreshold();
    final int maxIterationCount = detector.getMaxIterationCount();
    final BracketedUnivariateSolver<UnivariateFunction> solver = new BracketingNthOrderBrentSolver(0, convergence.getReal(), 0, 5);
    // event time, just at or before the actual root.
    FieldAbsoluteDate<T> beforeRootT = null;
    T beforeRootG = zero.add(Double.NaN);
    // time on the other side of the root.
    // Initialized the the loop below executes once.
    FieldAbsoluteDate<T> afterRootT = ta;
    T afterRootG = zero;
    // the ga == 0.0 case is handled by the loop below
    if (ta.equals(tb)) {
        // both non-zero but times are the same. Probably due to reset state
        beforeRootT = ta;
        beforeRootG = ga;
        afterRootT = shiftedBy(beforeRootT, convergence);
        afterRootG = g(interpolator.getInterpolatedState(afterRootT));
    } else if (ga.getReal() != 0.0 && gb.getReal() == 0.0) {
        // hard: ga != 0.0 and gb == 0.0
        // look past gb by up to convergence to find next sign
        // throw an exception if g(t) = 0.0 in [tb, tb + convergence]
        beforeRootT = tb;
        beforeRootG = gb;
        afterRootT = shiftedBy(beforeRootT, convergence);
        afterRootG = g(interpolator.getInterpolatedState(afterRootT));
    } else if (ga.getReal() != 0.0) {
        final T newGa = g(interpolator.getInterpolatedState(ta));
        if (ga.getReal() > 0 != newGa.getReal() > 0) {
            // both non-zero, step sign change at ta, possibly due to reset state
            beforeRootT = ta;
            beforeRootG = newGa;
            afterRootT = minTime(shiftedBy(beforeRootT, convergence), tb);
            afterRootG = g(interpolator.getInterpolatedState(afterRootT));
        }
    }
    // loop to skip through "fake" roots, i.e. where g(t) = g'(t) = 0.0
    // executed once if we didn't hit a special case above
    FieldAbsoluteDate<T> loopT = ta;
    T loopG = ga;
    while ((afterRootG.getReal() == 0.0 || afterRootG.getReal() > 0.0 == g0Positive) && strictlyAfter(afterRootT, tb)) {
        if (loopG.getReal() == 0.0) {
            // ga == 0.0 and gb may or may not be 0.0
            // handle the root at ta first
            beforeRootT = loopT;
            beforeRootG = loopG;
            afterRootT = minTime(shiftedBy(beforeRootT, convergence), tb);
            afterRootG = g(interpolator.getInterpolatedState(afterRootT));
        } else {
            // both non-zero, the usual case, use a root finder.
            try {
                // time zero for evaluating the function f. Needs to be final
                final FieldAbsoluteDate<T> fT0 = loopT;
                final UnivariateFunction f = dt -> {
                    try {
                        return g(interpolator.getInterpolatedState(fT0.shiftedBy(dt))).getReal();
                    } catch (OrekitException oe) {
                        throw new OrekitExceptionWrapper(oe);
                    }
                };
                // tb as a double for use in f
                final T tbDouble = tb.durationFrom(fT0);
                if (forward) {
                    final Interval interval = solver.solveInterval(maxIterationCount, f, 0, tbDouble.getReal());
                    beforeRootT = fT0.shiftedBy(interval.getLeftAbscissa());
                    beforeRootG = zero.add(interval.getLeftValue());
                    afterRootT = fT0.shiftedBy(interval.getRightAbscissa());
                    afterRootG = zero.add(interval.getRightValue());
                } else {
                    final Interval interval = solver.solveInterval(maxIterationCount, f, tbDouble.getReal(), 0);
                    beforeRootT = fT0.shiftedBy(interval.getRightAbscissa());
                    beforeRootG = zero.add(interval.getRightValue());
                    afterRootT = fT0.shiftedBy(interval.getLeftAbscissa());
                    afterRootG = zero.add(interval.getLeftValue());
                }
            } catch (OrekitExceptionWrapper oew) {
                throw oew.getException();
            }
        }
        // assume tolerance is 1 ulp
        if (beforeRootT.equals(afterRootT)) {
            afterRootT = nextAfter(afterRootT);
            afterRootG = g(interpolator.getInterpolatedState(afterRootT));
        }
        // check loop is making some progress
        check((forward && afterRootT.compareTo(beforeRootT) > 0) || (!forward && afterRootT.compareTo(beforeRootT) < 0));
        // setup next iteration
        loopT = afterRootT;
        loopG = afterRootG;
    }
    // figure out the result of root finding, and return accordingly
    if (afterRootG.getReal() == 0.0 || afterRootG.getReal() > 0.0 == g0Positive) {
        // loop gave up and didn't find any crossing within this step
        return false;
    } else {
        // real crossing
        check(beforeRootT != null && !Double.isNaN(beforeRootG.getReal()));
        // variation direction, with respect to the integration direction
        increasing = !g0Positive;
        pendingEventTime = beforeRootT;
        stopTime = beforeRootG.getReal() == 0.0 ? beforeRootT : afterRootT;
        pendingEvent = true;
        afterEvent = afterRootT;
        afterG = afterRootG;
        // check increasing set correctly
        check(afterG.getReal() > 0 == increasing);
        check(increasing == gb.getReal() >= ga.getReal());
        return true;
    }
}
Also used : BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) OrekitInternalError(org.orekit.errors.OrekitInternalError) FieldOrekitStepInterpolator(org.orekit.propagation.sampling.FieldOrekitStepInterpolator) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) Interval(org.hipparchus.analysis.solvers.BracketedUnivariateSolver.Interval) Precision(org.hipparchus.util.Precision) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) Field(org.hipparchus.Field) OrekitException(org.orekit.errors.OrekitException) RealFieldElement(org.hipparchus.RealFieldElement) Action(org.orekit.propagation.events.handlers.FieldEventHandler.Action) FieldEventHandler(org.orekit.propagation.events.handlers.FieldEventHandler) BracketedUnivariateSolver(org.hipparchus.analysis.solvers.BracketedUnivariateSolver) FastMath(org.hipparchus.util.FastMath) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) OrekitException(org.orekit.errors.OrekitException) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) Interval(org.hipparchus.analysis.solvers.BracketedUnivariateSolver.Interval)

Example 4 with OrekitExceptionWrapper

use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.

the class TransformGenerator method generate.

/**
 * {@inheritDoc}
 */
public List<Transform> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
    try {
        final List<Transform> generated = new ArrayList<>();
        if (existingDate == null) {
            // no prior existing transforms, just generate a first one
            for (int i = 0; i < neighborsSize; ++i) {
                generated.add(provider.getTransform(date.shiftedBy(i * step)));
            }
        } else {
            // some transforms have already been generated
            // add the missing ones up to specified date
            AbsoluteDate t = existingDate;
            if (date.compareTo(t) > 0) {
                // forward generation
                do {
                    t = t.shiftedBy(step);
                    generated.add(generated.size(), provider.getTransform(t));
                } while (t.compareTo(date) <= 0);
            } else {
                // backward generation
                do {
                    t = t.shiftedBy(-step);
                    generated.add(0, provider.getTransform(t));
                } while (t.compareTo(date) >= 0);
            }
        }
        // return the generated transforms
        return generated;
    } catch (OrekitException oe) {
        throw new OrekitExceptionWrapper(oe);
    }
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) ArrayList(java.util.ArrayList) OrekitException(org.orekit.errors.OrekitException) AbsoluteDate(org.orekit.time.AbsoluteDate)

Example 5 with OrekitExceptionWrapper

use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.

the class JacobianPropagatorConverter method getObjectiveFunction.

/**
 * {@inheritDoc}
 */
protected MultivariateVectorFunction getObjectiveFunction() {
    return new MultivariateVectorFunction() {

        /**
         * {@inheritDoc}
         */
        public double[] value(final double[] arg) throws IllegalArgumentException, OrekitExceptionWrapper {
            try {
                final double[] value = new double[getTargetSize()];
                final NumericalPropagator prop = builder.buildPropagator(arg);
                final int stateSize = isOnlyPosition() ? 3 : 6;
                final List<SpacecraftState> sample = getSample();
                for (int i = 0; i < sample.size(); ++i) {
                    final int row = i * stateSize;
                    if (prop.getInitialState().getDate().equals(sample.get(i).getDate())) {
                        // use initial state
                        fillRows(value, row, prop.getInitialState());
                    } else {
                        // use a date detector to pick up states
                        prop.addEventDetector(new DateDetector(sample.get(i).getDate()).withHandler(new EventHandler<DateDetector>() {

                            /**
                             * {@inheritDoc}
                             */
                            @Override
                            public Action eventOccurred(final SpacecraftState state, final DateDetector detector, final boolean increasing) throws OrekitException {
                                fillRows(value, row, state);
                                return row + stateSize >= getTargetSize() ? Action.STOP : Action.CONTINUE;
                            }
                        }));
                    }
                }
                prop.propagate(sample.get(sample.size() - 1).getDate().shiftedBy(10.0));
                return value;
            } catch (OrekitException ex) {
                throw new OrekitExceptionWrapper(ex);
            }
        }
    };
}
Also used : DateDetector(org.orekit.propagation.events.DateDetector) SpacecraftState(org.orekit.propagation.SpacecraftState) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) EventHandler(org.orekit.propagation.events.handlers.EventHandler) OrekitException(org.orekit.errors.OrekitException) MultivariateVectorFunction(org.hipparchus.analysis.MultivariateVectorFunction)

Aggregations

OrekitExceptionWrapper (org.orekit.errors.OrekitExceptionWrapper)27 OrekitException (org.orekit.errors.OrekitException)25 UnivariateFunction (org.hipparchus.analysis.UnivariateFunction)10 BracketingNthOrderBrentSolver (org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver)9 SpacecraftState (org.orekit.propagation.SpacecraftState)9 AbsoluteDate (org.orekit.time.AbsoluteDate)8 Map (java.util.Map)7 UnivariateSolver (org.hipparchus.analysis.solvers.UnivariateSolver)7 MathRuntimeException (org.hipparchus.exception.MathRuntimeException)6 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)6 Frame (org.orekit.frames.Frame)5 Transform (org.orekit.frames.Transform)5 ArrayList (java.util.ArrayList)4 HashMap (java.util.HashMap)3 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)3 FiniteDifferencesDifferentiator (org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)3 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)3 UnivariateVectorFunction (org.hipparchus.analysis.UnivariateVectorFunction)2 UnivariateDifferentiableVectorFunction (org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction)2 BracketedUnivariateSolver (org.hipparchus.analysis.solvers.BracketedUnivariateSolver)2