diff --git a/pom.xml b/pom.xml index 58fd4df442..e3bfd15ce1 100644 --- a/pom.xml +++ b/pom.xml @@ -5,7 +5,7 @@ org.orekit orekit jar - 12.1.2 + 12.1.3 OREKIT http://www.orekit.org/ diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 6766bfddad..46c265fe88 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -20,6 +20,49 @@ Orekit Changes + + + Removed non thread-safe use of DecimalFormat. + + + Pass Status in UKF theoretical measurement. + + + Fixed mixed up frames in inter-satellites measurements. + + + Protected several maps against concurrent modifications. + + + Fixed NeQuick ionospheric model for perfect zenith observation. + + + Fixed introduced noises when changing covariance frame with identical frame. + + + Fixed wrong agency name length in Rinex observation writer. + + + Greatly reduced computation time of NeQuick ionospheric model. + + + Added getOneLetterCode and getTwoLettersCode to TimeSystem. + + + Protected AmbiguityCache against concurrent modifications. + + + Fixed checkstyle error in SolarRadiationPressure. + + + Fixed update of initial state's prop type after a propagation in DSST. + + >(); } + /** Constructor from measurement base. + * @param estimatedMeasurementBase estimated measurement base + * @since 12.1.3 + */ + public EstimatedMeasurement(final EstimatedMeasurementBase estimatedMeasurementBase) { + this(estimatedMeasurementBase.getObservedMeasurement(), estimatedMeasurementBase.getIteration(), + estimatedMeasurementBase.getCount(), estimatedMeasurementBase.getStates(), + estimatedMeasurementBase.getParticipants()); + this.setEstimatedValue(estimatedMeasurementBase.getEstimatedValue()); + this.setStatus(estimatedMeasurementBase.getStatus()); + } + /** Get state size. *

* Warning, the {@link #setStateDerivatives(int, double[][])} diff --git a/src/main/java/org/orekit/estimation/measurements/EstimatedMeasurementBase.java b/src/main/java/org/orekit/estimation/measurements/EstimatedMeasurementBase.java index 6ad97d7086..a39ee1cc81 100644 --- a/src/main/java/org/orekit/estimation/measurements/EstimatedMeasurementBase.java +++ b/src/main/java/org/orekit/estimation/measurements/EstimatedMeasurementBase.java @@ -66,7 +66,7 @@ public class EstimatedMeasurementBase> implemen * @param count evaluations counter * @param states states of the spacecrafts * @param participants coordinates of the participants in signal travel order - * in inertial frame + * in inertial frame of first state */ public EstimatedMeasurementBase(final T observedMeasurement, final int iteration, final int count, @@ -124,7 +124,7 @@ public SpacecraftState[] getStates() { * spacecraft for two-way range measurement). *

* @return coordinates of the measurements participants in signal travel order - * in inertial frame + * in inertial frame of first state */ public TimeStampedPVCoordinates[] getParticipants() { return participants.clone(); diff --git a/src/main/java/org/orekit/estimation/measurements/InterSatellitesRange.java b/src/main/java/org/orekit/estimation/measurements/InterSatellitesRange.java index a11c86792e..edefa2d013 100644 --- a/src/main/java/org/orekit/estimation/measurements/InterSatellitesRange.java +++ b/src/main/java/org/orekit/estimation/measurements/InterSatellitesRange.java @@ -21,6 +21,7 @@ import java.util.Map; import org.hipparchus.analysis.differentiation.Gradient; +import org.orekit.frames.Frame; import org.orekit.propagation.SpacecraftState; import org.orekit.time.AbsoluteDate; import org.orekit.time.FieldAbsoluteDate; @@ -119,10 +120,11 @@ protected EstimatedMeasurementBase theoreticalEvaluationWi final SpacecraftState[] states) { // coordinates of both satellites + final Frame frame = states[0].getFrame(); final SpacecraftState local = states[0]; - final TimeStampedPVCoordinates pvaL = local.getPVCoordinates(); + final TimeStampedPVCoordinates pvaL = local.getPVCoordinates(frame); final SpacecraftState remote = states[1]; - final TimeStampedPVCoordinates pvaR = remote.getPVCoordinates(); + final TimeStampedPVCoordinates pvaR = remote.getPVCoordinates(frame); // compute propagation times // (if state has already been set up to pre-compensate propagation delay, @@ -134,7 +136,7 @@ protected EstimatedMeasurementBase theoreticalEvaluationWi final TimeStampedPVCoordinates s1Downlink = pvaL.shiftedBy(arrivalDate.durationFrom(pvaL.getDate())); - final double tauD = signalTimeOfFlight(pvaR, s1Downlink.getPosition(), arrivalDate, local.getFrame()); + final double tauD = signalTimeOfFlight(pvaR, s1Downlink.getPosition(), arrivalDate, frame); // Transit state final double delta = getDate().durationFrom(remote.getDate()); @@ -152,15 +154,15 @@ protected EstimatedMeasurementBase theoreticalEvaluationWi final double tauU = signalTimeOfFlight(pvaL, transitState.getPosition(), transitState.getDate(), - local.getFrame()); + frame); estimated = new EstimatedMeasurementBase<>(this, iteration, evaluation, new SpacecraftState[] { local.shiftedBy(deltaMTauD), remote.shiftedBy(deltaMTauD) }, new TimeStampedPVCoordinates[] { - local.shiftedBy(delta - tauD - tauU).getPVCoordinates(), - remote.shiftedBy(delta - tauD).getPVCoordinates(), - local.shiftedBy(delta).getPVCoordinates() + local.shiftedBy(delta - tauD - tauU).getPVCoordinates(frame), + remote.shiftedBy(delta - tauD).getPVCoordinates(frame), + local.shiftedBy(delta).getPVCoordinates(frame) }); // Range value @@ -173,8 +175,8 @@ protected EstimatedMeasurementBase theoreticalEvaluationWi local.shiftedBy(deltaMTauD), remote.shiftedBy(deltaMTauD) }, new TimeStampedPVCoordinates[] { - remote.shiftedBy(delta - tauD).getPVCoordinates(), - local.shiftedBy(delta).getPVCoordinates() + remote.shiftedBy(delta - tauD).getPVCoordinates(frame), + local.shiftedBy(delta).getPVCoordinates(frame) }); // Clock offsets @@ -196,6 +198,8 @@ protected EstimatedMeasurement theoreticalEvaluation(final final int evaluation, final SpacecraftState[] states) { + final Frame frame = states[0].getFrame(); + // Range derivatives are computed with respect to spacecraft states in inertial frame // ---------------------- // @@ -222,7 +226,10 @@ protected EstimatedMeasurement theoreticalEvaluation(final final SpacecraftState local = states[0]; final TimeStampedFieldPVCoordinates pvaL = getCoordinates(local, 0, nbParams); final SpacecraftState remote = states[1]; - final TimeStampedFieldPVCoordinates pvaR = getCoordinates(remote, 6, nbParams); + final TimeStampedFieldPVCoordinates pvaR = states[1]. + getFrame(). + getTransformTo(frame, states[1].getDate()). + transformPVCoordinates(getCoordinates(remote, 6, nbParams)); // compute propagation times // (if state has already been set up to pre-compensate propagation delay, @@ -236,7 +243,7 @@ protected EstimatedMeasurement theoreticalEvaluation(final final TimeStampedFieldPVCoordinates s1Downlink = pvaL.shiftedBy(arrivalDate.durationFrom(pvaL.getDate())); final Gradient tauD = signalTimeOfFlight(pvaR, s1Downlink.getPosition(), - arrivalDate, local.getFrame()); + arrivalDate, frame); // Transit state final double delta = getDate().durationFrom(remote.getDate()); @@ -254,15 +261,15 @@ protected EstimatedMeasurement theoreticalEvaluation(final final Gradient tauU = signalTimeOfFlight(pvaL, transitStateDS.getPosition(), transitStateDS.getDate(), - local.getFrame()); + frame); estimated = new EstimatedMeasurement<>(this, iteration, evaluation, new SpacecraftState[] { local.shiftedBy(deltaMTauD.getValue()), remote.shiftedBy(deltaMTauD.getValue()) }, new TimeStampedPVCoordinates[] { - local.shiftedBy(delta - tauD.getValue() - tauU.getValue()).getPVCoordinates(), - remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(), - local.shiftedBy(delta).getPVCoordinates() + local.shiftedBy(delta - tauD.getValue() - tauU.getValue()).getPVCoordinates(frame), + remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(frame), + local.shiftedBy(delta).getPVCoordinates(frame) }); // Range value @@ -275,7 +282,7 @@ protected EstimatedMeasurement theoreticalEvaluation(final local.shiftedBy(deltaMTauD.getValue()), remote.shiftedBy(deltaMTauD.getValue()) }, new TimeStampedPVCoordinates[] { - remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(), + remote.shiftedBy(delta - tauD.getValue()).getPVCoordinates(frame), local.shiftedBy(delta).getPVCoordinates() }); diff --git a/src/main/java/org/orekit/estimation/measurements/gnss/AbstractInterSatellitesMeasurement.java b/src/main/java/org/orekit/estimation/measurements/gnss/AbstractInterSatellitesMeasurement.java index b6f433a925..d7ed375dfa 100644 --- a/src/main/java/org/orekit/estimation/measurements/gnss/AbstractInterSatellitesMeasurement.java +++ b/src/main/java/org/orekit/estimation/measurements/gnss/AbstractInterSatellitesMeasurement.java @@ -83,7 +83,10 @@ protected FieldPVCoordinatesProvider getRemotePV(final SpacecraftState final TimeStampedFieldPVCoordinates pv0 = getCoordinates(states[1], 6, freeParameters); // shift to desired date - return pv0.shiftedBy(date.durationFrom(states[1].getDate())); + final TimeStampedFieldPVCoordinates shifted = pv0.shiftedBy(date.durationFrom(states[1].getDate())); + + // transform to desired frame + return states[1].getFrame().getTransformTo(frame, states[1].getDate()).transformPVCoordinates(shifted); }; } diff --git a/src/main/java/org/orekit/estimation/measurements/gnss/AbstractOnBoardMeasurement.java b/src/main/java/org/orekit/estimation/measurements/gnss/AbstractOnBoardMeasurement.java index 1e4237744d..db3e565198 100644 --- a/src/main/java/org/orekit/estimation/measurements/gnss/AbstractOnBoardMeasurement.java +++ b/src/main/java/org/orekit/estimation/measurements/gnss/AbstractOnBoardMeasurement.java @@ -23,6 +23,7 @@ import org.orekit.estimation.measurements.ObservedMeasurement; import org.orekit.estimation.measurements.QuadraticClockModel; import org.orekit.estimation.measurements.QuadraticFieldClockModel; +import org.orekit.frames.Frame; import org.orekit.propagation.SpacecraftState; import org.orekit.time.AbsoluteDate; import org.orekit.time.ClockOffset; @@ -109,7 +110,8 @@ protected OnBoardCommonParametersWithoutDerivatives computeCommonParametersWitho final boolean clockOffsetAlreadyApplied) { // local and remote satellites - final TimeStampedPVCoordinates pvaLocal = states[0].getPVCoordinates(); + final Frame frame = states[0].getFrame(); + final TimeStampedPVCoordinates pvaLocal = states[0].getPVCoordinates(frame); final ClockOffset localClock = getSatellites(). get(0). getQuadraticClockModel(). @@ -125,7 +127,7 @@ protected OnBoardCommonParametersWithoutDerivatives computeCommonParametersWitho final double deltaT = arrivalDate.durationFrom(states[0]); final TimeStampedPVCoordinates pvaDownlink = pvaLocal.shiftedBy(deltaT); final double tauD = signalTimeOfFlight(remotePV, arrivalDate, pvaDownlink.getPosition(), - arrivalDate, states[0].getFrame()); + arrivalDate, frame); // Remote satellite at signal emission final AbsoluteDate emissionDate = arrivalDate.shiftedBy(-tauD); @@ -136,7 +138,7 @@ protected OnBoardCommonParametersWithoutDerivatives computeCommonParametersWitho localClockOffset, localClockRate, remoteClockOffset, remoteClockRate, tauD, pvaDownlink, - remotePV.getPVCoordinates(emissionDate, states[0].getFrame())); + remotePV.getPVCoordinates(emissionDate, frame)); } @@ -150,6 +152,8 @@ protected OnBoardCommonParametersWithoutDerivatives computeCommonParametersWitho protected OnBoardCommonParametersWithDerivatives computeCommonParametersWith(final SpacecraftState[] states, final boolean clockOffsetAlreadyApplied) { + final Frame frame = states[0].getFrame(); + // measurement derivatives are computed with respect to spacecraft state in inertial frame // Parameters: // - 6k..6k+2 - Position of spacecraft k (counting k from 0 to nbSat-1) in inertial frame @@ -183,7 +187,7 @@ protected OnBoardCommonParametersWithDerivatives computeCommonParametersWith(fin final TimeStampedFieldPVCoordinates pvaDownlink = pvaLocal.shiftedBy(deltaT); final Gradient tauD = signalTimeOfFlight(remotePV, arrivalDate, pvaDownlink.getPosition(), arrivalDate, - states[0].getFrame()); + frame); // Remote satellite at signal emission final FieldAbsoluteDate emissionDate = arrivalDate.shiftedBy(tauD.negate()); @@ -193,7 +197,7 @@ protected OnBoardCommonParametersWithDerivatives computeCommonParametersWith(fin localClockOffset.getOffset(), localClockOffset.getRate(), remoteClockOffset.getOffset(), remoteClockOffset.getRate(), tauD, pvaDownlink, - remotePV.getPVCoordinates(emissionDate, states[0].getFrame())); + remotePV.getPVCoordinates(emissionDate, frame)); } diff --git a/src/main/java/org/orekit/estimation/measurements/gnss/AmbiguityCache.java b/src/main/java/org/orekit/estimation/measurements/gnss/AmbiguityCache.java index 575254b05b..f89c87cedf 100644 --- a/src/main/java/org/orekit/estimation/measurements/gnss/AmbiguityCache.java +++ b/src/main/java/org/orekit/estimation/measurements/gnss/AmbiguityCache.java @@ -55,8 +55,10 @@ public AmbiguityCache() { * @return parameter driver for the emitter/receiver/wavelength triplet */ public AmbiguityDriver getAmbiguity(final String emitter, final String receiver, final double wavelength) { - return cache.computeIfAbsent(new Key(emitter, receiver, wavelength), - k -> new AmbiguityDriver(emitter, receiver, wavelength)); + final Key key = new Key(emitter, receiver, wavelength); + synchronized (cache) { + return cache.computeIfAbsent(key, k -> new AmbiguityDriver(emitter, receiver, wavelength)); + } } /** Key for the map. */ diff --git a/src/main/java/org/orekit/estimation/measurements/gnss/WindUpFactory.java b/src/main/java/org/orekit/estimation/measurements/gnss/WindUpFactory.java index 72a4dba4a7..536771fe87 100644 --- a/src/main/java/org/orekit/estimation/measurements/gnss/WindUpFactory.java +++ b/src/main/java/org/orekit/estimation/measurements/gnss/WindUpFactory.java @@ -51,15 +51,18 @@ public WindUpFactory() { public WindUp getWindUp(final SatelliteSystem system, final int prnNumber, final Dipole emitterDipole, final String receiverName) { // select satellite system - final Map> systemModifiers = - modifiers.computeIfAbsent(system, s -> new HashMap<>()); + final Map> systemModifiers; + synchronized (modifiers) { + systemModifiers = modifiers.computeIfAbsent(system, s -> new HashMap<>()); - // select satellite - final Map satelliteModifiers = + // select satellite + final Map satelliteModifiers = systemModifiers.computeIfAbsent(prnNumber, n -> new HashMap<>()); - // select receiver - return satelliteModifiers.computeIfAbsent(receiverName, r -> new WindUp(emitterDipole)); + // select receiver + return satelliteModifiers.computeIfAbsent(receiverName, r -> new WindUp(emitterDipole)); + + } } diff --git a/src/main/java/org/orekit/estimation/sequential/KalmanEstimationCommon.java b/src/main/java/org/orekit/estimation/sequential/KalmanEstimationCommon.java index a0ee0909b0..6da3867cc2 100644 --- a/src/main/java/org/orekit/estimation/sequential/KalmanEstimationCommon.java +++ b/src/main/java/org/orekit/estimation/sequential/KalmanEstimationCommon.java @@ -16,6 +16,13 @@ */ package org.orekit.estimation.sequential; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Comparator; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + import org.hipparchus.filtering.kalman.ProcessEstimate; import org.hipparchus.linear.ArrayRealVector; import org.hipparchus.linear.MatrixUtils; @@ -30,13 +37,6 @@ import org.orekit.utils.ParameterDriversList; import org.orekit.utils.ParameterDriversList.DelegatingDriver; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Comparator; -import java.util.HashMap; -import java.util.List; -import java.util.Map; - /** Class defining the process model dynamics to use with a {@link KalmanEstimator}. * @author Romain Gerbaud * @author Maxime Journot @@ -473,6 +473,11 @@ public Propagator[] getEstimatedPropagators() { return propagators; } + /** Get the normalized process noise matrix. + * + * @param stateDimension state dimension + * @return the normalized process noise matrix + */ protected RealMatrix getNormalizedProcessNoise(final int stateDimension) { final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(stateDimension, stateDimension); for (int k = 0; k < covarianceMatricesProviders.size(); ++k) { @@ -519,67 +524,116 @@ protected RealMatrix getNormalizedProcessNoise(final int stateDimension) { return KalmanEstimatorUtil.normalizeCovarianceMatrix(physicalProcessNoise, scale); } - + /** Getter for the orbitsStartColumns. + * @return the orbitsStartColumns + */ protected int[] getOrbitsStartColumns() { return orbitsStartColumns; } + /** Getter for the propagationParameterColumns. + * @return the propagationParameterColumns + */ protected Map getPropagationParameterColumns() { return propagationParameterColumns; } + /** Getter for the measurementParameterColumns. + * @return the measurementParameterColumns + */ protected Map getMeasurementParameterColumns() { return measurementParameterColumns; } + /** Getter for the estimatedPropagationParameters. + * @return the estimatedPropagationParameters + */ protected ParameterDriversList[] getEstimatedPropagationParametersArray() { return estimatedPropagationParameters; } + /** Getter for the estimatedOrbitalParameters. + * @return the estimatedOrbitalParameters + */ protected ParameterDriversList[] getEstimatedOrbitalParametersArray() { return estimatedOrbitalParameters; } + /** Getter for the covarianceIndirection. + * @return the covarianceIndirection + */ protected int[][] getCovarianceIndirection() { return covarianceIndirection; } + /** Getter for the scale. + * @return the scale + */ protected double[] getScale() { return scale; } + /** Getter for the correctedEstimate. + * @return the correctedEstimate + */ protected ProcessEstimate getCorrectedEstimate() { return correctedEstimate; } + /** Setter for the correctedEstimate. + * @param correctedEstimate the correctedEstimate + */ protected void setCorrectedEstimate(final ProcessEstimate correctedEstimate) { this.correctedEstimate = correctedEstimate; } + /** Getter for the referenceDate. + * @return the referenceDate + */ protected AbsoluteDate getReferenceDate() { return referenceDate; } + /** Increment current measurement number. */ protected void incrementCurrentMeasurementNumber() { currentMeasurementNumber += 1; } + /** Setter for the currentDate. + * @param currentDate the currentDate + */ public void setCurrentDate(final AbsoluteDate currentDate) { this.currentDate = currentDate; } + /** Set correctedSpacecraftState at index. + * + * @param correctedSpacecraftState corrected S/C state o set + * @param index index where to set in the array + */ protected void setCorrectedSpacecraftState(final SpacecraftState correctedSpacecraftState, final int index) { this.correctedSpacecraftStates[index] = correctedSpacecraftState; } + /** Set predictedSpacecraftState at index. + * + * @param predictedSpacecraftState predicted S/C state o set + * @param index index where to set in the array + */ protected void setPredictedSpacecraftState(final SpacecraftState predictedSpacecraftState, final int index) { this.predictedSpacecraftStates[index] = predictedSpacecraftState; } + /** Setter for the predictedMeasurement. + * @param predictedMeasurement the predictedMeasurement + */ protected void setPredictedMeasurement(final EstimatedMeasurement predictedMeasurement) { this.predictedMeasurement = predictedMeasurement; } + /** Setter for the correctedMeasurement. + * @param correctedMeasurement the correctedMeasurement + */ protected void setCorrectedMeasurement(final EstimatedMeasurement correctedMeasurement) { this.correctedMeasurement = correctedMeasurement; } diff --git a/src/main/java/org/orekit/estimation/sequential/UnscentedKalmanModel.java b/src/main/java/org/orekit/estimation/sequential/UnscentedKalmanModel.java index c446593b55..e30f01c453 100644 --- a/src/main/java/org/orekit/estimation/sequential/UnscentedKalmanModel.java +++ b/src/main/java/org/orekit/estimation/sequential/UnscentedKalmanModel.java @@ -336,11 +336,7 @@ private static > EstimatedMeasurement estima final EstimatedMeasurementBase estimatedMeasurementBase = observedMeasurement. estimateWithoutDerivatives(measurementNumber, measurementNumber, KalmanEstimatorUtil.filterRelevant(observedMeasurement, spacecraftStates)); - final EstimatedMeasurement estimatedMeasurement = new EstimatedMeasurement<>(estimatedMeasurementBase.getObservedMeasurement(), - estimatedMeasurementBase.getIteration(), estimatedMeasurementBase.getCount(), - estimatedMeasurementBase.getStates(), estimatedMeasurementBase.getParticipants()); - estimatedMeasurement.setEstimatedValue(estimatedMeasurementBase.getEstimatedValue()); - return estimatedMeasurement; + return new EstimatedMeasurement<>(estimatedMeasurementBase); } /** Update parameter drivers with a normalised state, adjusting state according to the driver limits. diff --git a/src/main/java/org/orekit/files/ilrs/CPF.java b/src/main/java/org/orekit/files/ilrs/CPF.java index aaf0199108..99581fe4df 100644 --- a/src/main/java/org/orekit/files/ilrs/CPF.java +++ b/src/main/java/org/orekit/files/ilrs/CPF.java @@ -56,13 +56,13 @@ public class CPF implements EphemerisFile { private CartesianDerivativesFilter filter; /** CPF file header. */ - private CPFHeader header; + private final CPFHeader header; /** Map containing satellite information. */ - private Map ephemeris; + private final Map ephemeris; /** List of comments contained in the file. */ - private List comments; + private final List comments; /** * Constructor. @@ -114,7 +114,11 @@ public List getComments() { * @since 11.0.1 */ public void addSatelliteCoordinates(final String id, final List coord) { - ephemeris.computeIfAbsent(id, i -> new CPFEphemeris(i)).coordinates.addAll(coord); + final CPFEphemeris e; + synchronized (this) { + e = ephemeris.computeIfAbsent(id, CPFEphemeris::new); + } + e.coordinates.addAll(coord); } /** @@ -124,7 +128,11 @@ public void addSatelliteCoordinates(final String id, final List c * @since 11.0.1 */ public void addSatelliteCoordinate(final String id, final CPFCoordinate coord) { - ephemeris.computeIfAbsent(id, i -> new CPFEphemeris(i)).coordinates.add(coord); + final CPFEphemeris e; + synchronized (this) { + e = ephemeris.computeIfAbsent(id, CPFEphemeris::new); + } + e.coordinates.add(coord); } /** diff --git a/src/main/java/org/orekit/files/rinex/clock/RinexClock.java b/src/main/java/org/orekit/files/rinex/clock/RinexClock.java index cb36aca64e..764494c5ab 100644 --- a/src/main/java/org/orekit/files/rinex/clock/RinexClock.java +++ b/src/main/java/org/orekit/files/rinex/clock/RinexClock.java @@ -403,9 +403,11 @@ public Map> getSystemObservationTypes() { */ public void addSystemObservationType(final SatelliteSystem satSystem, final ObservationType observationType) { - systemObservationTypes. - computeIfAbsent(satSystem, s -> new ArrayList<>()). - add(observationType); + final List list; + synchronized (systemObservationTypes) { + list = systemObservationTypes.computeIfAbsent(satSystem, s -> new ArrayList<>()); + } + list.add(observationType); } /** Getter for the file time system. @@ -663,7 +665,11 @@ public Map> getClockData() { */ public void addClockData(final String id, final ClockDataLine clockDataLine) { - clockData.computeIfAbsent(id, i -> new ArrayList<>()).add(clockDataLine); + final List list; + synchronized (clockData) { + list = clockData.computeIfAbsent(id, i -> new ArrayList<>()); + } + list.add(clockDataLine); final AbsoluteDate epoch = clockDataLine.getEpoch(); if (epoch.isBefore(earliestEpoch)) { earliestEpoch = epoch; diff --git a/src/main/java/org/orekit/files/rinex/observation/RinexObservationHeader.java b/src/main/java/org/orekit/files/rinex/observation/RinexObservationHeader.java index 869f43a397..0c848501fd 100644 --- a/src/main/java/org/orekit/files/rinex/observation/RinexObservationHeader.java +++ b/src/main/java/org/orekit/files/rinex/observation/RinexObservationHeader.java @@ -706,8 +706,10 @@ public List getPhaseShiftCorrections() { * @param scaleFactorCorrection scale factor correction */ public void addScaleFactorCorrection(final SatelliteSystem satelliteSystem, final ScaleFactorCorrection scaleFactorCorrection) { - final List sfc = scaleFactorCorrections.computeIfAbsent(satelliteSystem, - k -> new ArrayList<>()); + final List sfc; + synchronized (scaleFactorCorrections) { + sfc = scaleFactorCorrections.computeIfAbsent(satelliteSystem, k -> new ArrayList<>()); + } sfc.add(scaleFactorCorrection); } @@ -759,7 +761,10 @@ public int getNbSat() { * @since 12.0 */ public void setNbObsPerSatellite(final SatInSystem sat, final ObservationType type, final int nbObs) { - final Map satNbObs = nbObsPerSat.computeIfAbsent(sat, k -> new HashMap<>()); + final Map satNbObs; + synchronized (nbObsPerSat) { + satNbObs = nbObsPerSat.computeIfAbsent(sat, k -> new HashMap<>()); + } satNbObs.put(type, nbObs); } diff --git a/src/main/java/org/orekit/files/rinex/observation/RinexObservationWriter.java b/src/main/java/org/orekit/files/rinex/observation/RinexObservationWriter.java index 0cadd0e273..e62867cbdc 100644 --- a/src/main/java/org/orekit/files/rinex/observation/RinexObservationWriter.java +++ b/src/main/java/org/orekit/files/rinex/observation/RinexObservationWriter.java @@ -282,7 +282,7 @@ public void writeHeader(final RinexObservationHeader header) // OBSERVER / AGENCY outputField(header.getObserverName(), 20, true); - outputField(header.getAgencyName(), 40, true); + outputField(header.getAgencyName(), 60, true); finishHeaderLine(RinexLabels.OBSERVER_AGENCY); // REC # / TYPE / VERS @@ -883,6 +883,9 @@ private void outputField(final String format, final double value, final int next */ private void outputField(final String field, final int next, final boolean leftJustified) throws IOException { final int padding = next - (field == null ? 0 : field.length()) - column; + if (padding < 0) { + throw new OrekitException(OrekitMessages.FIELD_TOO_LONG, field, next - column); + } if (leftJustified && field != null) { output.append(field); } diff --git a/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java b/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java index 58b55a68cf..255b8c6d30 100644 --- a/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java +++ b/src/main/java/org/orekit/forces/radiation/SolarRadiationPressure.java @@ -16,6 +16,9 @@ */ package org.orekit.forces.radiation; +import java.lang.reflect.Array; +import java.util.List; + import org.hipparchus.CalculusFieldElement; import org.hipparchus.geometry.euclidean.threed.FieldVector3D; import org.hipparchus.geometry.euclidean.threed.Vector3D; @@ -33,9 +36,6 @@ import org.orekit.utils.OccultationEngine; import org.orekit.utils.ParameterDriver; -import java.lang.reflect.Array; -import java.util.List; - /** Solar radiation pressure force model. *

* Since Orekit 11.0, it is possible to take into account @@ -333,7 +333,7 @@ private > T getLightingRatio(final FieldSpacec @SuppressWarnings("unchecked") final OccultationEngine.FieldOccultationAngles[] angles = - (OccultationEngine.FieldOccultationAngles[]) Array.newInstance(OccultationEngine.FieldOccultationAngles.class, n); + (OccultationEngine.FieldOccultationAngles[]) Array.newInstance(OccultationEngine.FieldOccultationAngles.class, n); for (int i = 0; i < n; ++i) { angles[i] = occultingBodies.get(i).angles(state); } diff --git a/src/main/java/org/orekit/frames/FixedTransformProvider.java b/src/main/java/org/orekit/frames/FixedTransformProvider.java index 7b21012657..65fd8a7de5 100644 --- a/src/main/java/org/orekit/frames/FixedTransformProvider.java +++ b/src/main/java/org/orekit/frames/FixedTransformProvider.java @@ -21,8 +21,8 @@ import java.util.HashMap; import java.util.Map; -import org.hipparchus.Field; import org.hipparchus.CalculusFieldElement; +import org.hipparchus.Field; import org.orekit.time.AbsoluteDate; import org.orekit.time.FieldAbsoluteDate; @@ -55,16 +55,13 @@ public Transform getTransform(final AbsoluteDate date) { } /** {@inheritDoc} */ + @SuppressWarnings("unchecked") @Override public > FieldTransform getTransform(final FieldAbsoluteDate date) { - - @SuppressWarnings("unchecked") - final FieldTransform ft = - (FieldTransform) cached.computeIfAbsent(date.getField(), - f -> new FieldTransform<>((Field) f, transform)); - - return ft; - + synchronized (cached) { + return (FieldTransform) cached.computeIfAbsent(date.getField(), + f -> new FieldTransform<>((Field) f, transform)); + } } /** Replace the instance with a data transfer object for serialization. diff --git a/src/main/java/org/orekit/frames/MODProvider.java b/src/main/java/org/orekit/frames/MODProvider.java index 1bb9655169..ef3201bcba 100644 --- a/src/main/java/org/orekit/frames/MODProvider.java +++ b/src/main/java/org/orekit/frames/MODProvider.java @@ -20,8 +20,8 @@ import java.util.HashMap; import java.util.Map; -import org.hipparchus.Field; import org.hipparchus.CalculusFieldElement; +import org.hipparchus.Field; import org.hipparchus.geometry.euclidean.threed.FieldRotation; import org.hipparchus.geometry.euclidean.threed.Rotation; import org.hipparchus.geometry.euclidean.threed.RotationConvention; @@ -94,16 +94,18 @@ public Transform getTransform(final AbsoluteDate date) { } /** {@inheritDoc} */ + @SuppressWarnings("unchecked") @Override public > FieldTransform getTransform(final FieldAbsoluteDate date) { // compute the precession angles phiA, omegaA, chiA final T[] angles = precessionFunction.value(date); - @SuppressWarnings("unchecked") - final FieldRotation fR4 = - (FieldRotation) fieldR4.computeIfAbsent(date.getField(), - f -> new FieldRotation<>((Field) f, r4)); + final FieldRotation fR4; + synchronized (fieldR4) { + fR4 = (FieldRotation) fieldR4.computeIfAbsent(date.getField(), + f -> new FieldRotation<>((Field) f, r4)); + } // complete precession final FieldRotation precession = fR4.compose(new FieldRotation<>(RotationOrder.ZXZ, RotationConvention.FRAME_TRANSFORM, diff --git a/src/main/java/org/orekit/gnss/TimeSystem.java b/src/main/java/org/orekit/gnss/TimeSystem.java index 3ede2749ac..89690ae14a 100644 --- a/src/main/java/org/orekit/gnss/TimeSystem.java +++ b/src/main/java/org/orekit/gnss/TimeSystem.java @@ -63,12 +63,12 @@ public enum TimeSystem { */ SBAS(null, "SB", "S", TimeScales::getUTC), - /** GMT (should only by used in RUN BY / DATE entries). + /** GMT (should only be used in RUN BY / DATE entries). * @since 12.0 */ GMT("GMT", null, null, TimeScales::getUTC), - /** Unknown (should only by used in RUN BY / DATE entries). */ + /** Unknown (should only be used in RUN BY / DATE entries). */ UNKNOWN("LCL", null, null, TimeScales::getGPS); /** Parsing key map. */ @@ -138,6 +138,22 @@ public String getKey() { return key; } + /** Get the two letters code. + * @return two letters code (may be null for non-GNSS time systems) + * @since 12.2 + */ + public String getTwoLettersCode() { + return twoLettersCode; + } + + /** Get the one letter code. + * @return one letter code (may be null for non-GNSS time systems) + * @since 12.2 + */ + public String getOneLetterCode() { + return oneLetterCode; + } + /** Parse a string to get the time system. *

* The string must be the time system. diff --git a/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java b/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java index 061270c780..0908387e4f 100644 --- a/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java +++ b/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java @@ -21,7 +21,6 @@ import org.hipparchus.util.FastMath; import org.hipparchus.util.FieldSinCos; import org.hipparchus.util.MathArrays; -import org.hipparchus.util.SinCos; import org.orekit.time.DateComponents; import org.orekit.time.DateTimeComponents; import org.orekit.time.TimeComponents; @@ -76,21 +75,21 @@ class FieldNeQuickParameters > { /** * Build a new instance. - * @param field field of the elements * @param dateTime current date time components - * @param f2 F2 coefficients used by the F2 layer - * @param fm3 Fm3 coefficients used by the F2 layer + * @param flattenF2 F2 coefficients used by the F2 layer (flatten array) + * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array) * @param latitude latitude of a point along the integration path, in radians * @param longitude longitude of a point along the integration path, in radians * @param alpha effective ionisation level coefficients * @param modipGrip modip grid */ - FieldNeQuickParameters(final Field field, final DateTimeComponents dateTime, final double[][][] f2, - final double[][][] fm3, final T latitude, final T longitude, + FieldNeQuickParameters(final DateTimeComponents dateTime, + final double[] flattenF2, final double[] flattenFm3, + final T latitude, final T longitude, final double[] alpha, final double[][] modipGrip) { // Zero - final T zero = field.getZero(); + final T zero = latitude.getField().getZero(); // MODIP in degrees final T modip = computeMODIP(latitude, longitude, modipGrip); @@ -106,25 +105,8 @@ class FieldNeQuickParameters > { // Effective solar zenith angle in radians final T xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude); - // Coefficients for F2 layer parameters - // Compute the array of interpolated coefficients for foF2 (Eq. 44) - final T[][] af2 = MathArrays.buildArray(field, 76, 13); - for (int j = 0; j < 76; j++) { - for (int k = 0; k < 13; k++ ) { - af2[j][k] = azr.multiply(0.01).negate().add(1.0).multiply(f2[0][j][k]).add(azr.multiply(0.01).multiply(f2[1][j][k])); - } - } - - // Compute the array of interpolated coefficients for M(3000)F2 (Eq. 46) - final T[][] am3 = MathArrays.buildArray(field, 49, 9); - for (int j = 0; j < 49; j++) { - for (int k = 0; k < 9; k++ ) { - am3[j][k] = azr.multiply(0.01).negate().add(1.0).multiply(fm3[0][j][k]).add(azr.multiply(0.01).multiply(fm3[1][j][k])); - } - } - // E layer maximum density height in km (Eq. 78) - this.hmE = field.getZero().newInstance(120.0); + this.hmE = zero.newInstance(120.0); // E layer critical frequency in MHz final T foE = computefoE(date.getMonth(), az, xeff, latitude); // E layer maximum density in 10^11 m-3 (Eq. 36) @@ -133,19 +115,21 @@ class FieldNeQuickParameters > { // Time argument (Eq. 49) final double t = FastMath.toRadians(15 * hours) - FastMath.PI; // Compute Fourier time series for foF2 and M(3000)F2 - final T[] cf2 = computeCF2(field, af2, t); - final T[] cm3 = computeCm3(field, am3, t); + final T[] scT = sinCos(zero.newInstance(t), 6); + final T[] cf2 = computeCF2(flattenF2, azr, scT); + final T[] cm3 = computeCm3(flattenFm3, azr, scT); // F2 layer critical frequency in MHz - final T foF2 = computefoF2(field, modip, cf2, latitude, longitude); + final T[] scL = sinCos(longitude, 8); + final T foF2 = computefoF2(modip, cf2, latitude, scL); // Maximum Usable Frequency factor - final T mF2 = computeMF2(field, modip, cm3, latitude, longitude); + final T mF2 = computeMF2(modip, cm3, latitude, scL); // F2 layer maximum density in 10^11 m-3 this.nmF2 = foF2.multiply(foF2).multiply(0.124); // F2 layer maximum density height in km - this.hmF2 = computehmF2(field, foE, foF2, mF2); + this.hmF2 = computehmF2(foE, foF2, mF2); // F1 layer critical frequency in MHz - final T foF1 = computefoF1(field, foE, foF2); + final T foF1 = computefoF1(foE, foF2); // F1 layer maximum density in 10^11 m-3 final T nmF1; if (foF1.getReal() <= 0.0 && foE.getReal() > 2.0) { @@ -166,10 +150,10 @@ class FieldNeQuickParameters > { this.beBot = zero.newInstance(5.0); // Layer amplitude coefficients - this.amplitudes = computeLayerAmplitudes(field, nmE, nmF1, foF1); + this.amplitudes = computeLayerAmplitudes(nmE, nmF1, foF1); // Topside thickness parameter - this.h0 = computeH0(field, date.getMonth(), azr); + this.h0 = computeH0(date.getMonth(), azr); } /** @@ -323,9 +307,8 @@ private T computeMODIP(final T lat, final T lon, final double[][] stModip) { final T y = b.subtract(bF); // MODIP (Ref Eq. 16) - final T modip = interpolate(z1, z2, z3, z4, y); + return interpolate(z1, z2, z3, z4, y); - return modip; } /** @@ -425,21 +408,20 @@ private T computefoE(final int month, final T az, final T seasp = ee.subtract(1.0).divide(ee.add(1.0)).multiply(seas); // Critical frequency (Eq. 35) final T coef = seasp.multiply(0.019).negate().add(1.112); - final T foE = FastMath.sqrt(coef .multiply(coef).multiply(sqAz).multiply(FastMath.cos(xeff).pow(0.6)).add(0.49)); - return foE; + return FastMath.sqrt(coef .multiply(coef).multiply(sqAz).multiply(FastMath.cos(xeff).pow(0.6)).add(0.49)); + } /** * Computes the F2 layer height of maximum electron density. - * @param field field of the elements * @param foE E layer layer critical frequency in MHz * @param foF2 F2 layer layer critical frequency in MHz * @param mF2 maximum usable frequency factor * @return hmF2 in km */ - private T computehmF2(final Field field, final T foE, final T foF2, final T mF2) { + private T computehmF2(final T foE, final T foF2, final T mF2) { // Zero - final T zero = field.getZero(); + final T zero = foE.getField().getZero(); // Ratio final T fo = foF2.divide(foE); final T ratio = join(fo, zero.newInstance(1.75), zero.newInstance(20.0), fo.subtract(1.75)); @@ -453,163 +435,192 @@ private T computehmF2(final Field field, final T foE, final T foF2, final T m // hmF2 Eq. 80 final T mF2Sq = mF2.square(); final T temp = FastMath.sqrt(mF2Sq.multiply(0.0196).add(1.0).divide(mF2Sq.multiply(1.2967).subtract(1.0))); - final T height = mF2.multiply(1490.0).multiply(temp).divide(mF2.add(deltaM)).subtract(176.0); - return height; + return mF2.multiply(1490.0).multiply(temp).divide(mF2.add(deltaM)).subtract(176.0); + + } + + /** Compute sines and cosines. + * @param a argument + * @param n number of terms + * @return sin(a), cos(a), sin(2a), cos(2a) … sin(n a), cos(n a) array + * @since 12.1.3 + */ + private T[] sinCos(final T a, final int n) { + + final FieldSinCos sc0 = FastMath.sinCos(a); + FieldSinCos sci = sc0; + final T[] sc = MathArrays.buildArray(a.getField(), 2 * n); + int isc = 0; + sc[isc++] = sci.sin(); + sc[isc++] = sci.cos(); + for (int i = 1; i < n; i++) { + sci = FieldSinCos.sum(sc0, sci); + sc[isc++] = sci.sin(); + sc[isc++] = sci.cos(); + } + + return sc; + } /** * Computes cf2 coefficients. - * @param field field of the elements - * @param af2 interpolated coefficients for foF2 - * @param t time argument + * @param flattenF2 F2 coefficients used by the F2 layer (flatten array) + * @param azr effective sunspot number (Eq. 19) + * @param scT sines/cosines array of time argument * @return the cf2 coefficients array */ - private T[] computeCF2(final Field field, final T[][] af2, final double t) { - // Eq. 50 - final T[] cf2 = MathArrays.buildArray(field, 76); + private T[] computeCF2(final double[] flattenF2, final T azr, final T[] scT) { + + // interpolation coefficients for effective spot number + final T azr01 = azr.multiply(0.01); + final T omazr01 = azr01.negate().add(1); + + // Eq. 44 and Eq. 50 merged into one loop + final T[] cf2 = MathArrays.buildArray(azr.getField(), 76); + int index = 0; for (int i = 0; i < cf2.length; i++) { - T sum = field.getZero(); - for (int k = 0; k < 6; k++) { - final SinCos sc = FastMath.sinCos((k + 1) * t); - sum = sum.add(af2[i][2 * k + 1].multiply(sc.sin()).add(af2[i][2 * (k + 1)].multiply(sc.cos()))); - } - cf2[i] = af2[i][0].add(sum); + // CHECKSTYLE: stop Indentation check + cf2[i] = omazr01.multiply(flattenF2[index ]).add(azr01.multiply(flattenF2[index + 1])). + add(omazr01.multiply(flattenF2[index + 2]).add(azr01.multiply(flattenF2[index + 3])).multiply(scT[ 0])). + add(omazr01.multiply(flattenF2[index + 4]).add(azr01.multiply(flattenF2[index + 5])).multiply(scT[ 1])). + add(omazr01.multiply(flattenF2[index + 6]).add(azr01.multiply(flattenF2[index + 7])).multiply(scT[ 2])). + add(omazr01.multiply(flattenF2[index + 8]).add(azr01.multiply(flattenF2[index + 9])).multiply(scT[ 3])). + add(omazr01.multiply(flattenF2[index + 10]).add(azr01.multiply(flattenF2[index + 11])).multiply(scT[ 4])). + add(omazr01.multiply(flattenF2[index + 12]).add(azr01.multiply(flattenF2[index + 13])).multiply(scT[ 5])). + add(omazr01.multiply(flattenF2[index + 14]).add(azr01.multiply(flattenF2[index + 15])).multiply(scT[ 6])). + add(omazr01.multiply(flattenF2[index + 16]).add(azr01.multiply(flattenF2[index + 17])).multiply(scT[ 7])). + add(omazr01.multiply(flattenF2[index + 18]).add(azr01.multiply(flattenF2[index + 19])).multiply(scT[ 8])). + add(omazr01.multiply(flattenF2[index + 20]).add(azr01.multiply(flattenF2[index + 21])).multiply(scT[ 9])). + add(omazr01.multiply(flattenF2[index + 22]).add(azr01.multiply(flattenF2[index + 23])).multiply(scT[10])). + add(omazr01.multiply(flattenF2[index + 24]).add(azr01.multiply(flattenF2[index + 25])).multiply(scT[11])); + index += 26; + // CHECKSTYLE: resume Indentation check } return cf2; } /** * Computes Cm3 coefficients. - * @param field field of the elements - * @param am3 interpolated coefficients for foF2 - * @param t time argument + * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array) + * @param azr effective sunspot number (Eq. 19) + * @param scT sines/cosines array of time argument * @return the Cm3 coefficients array */ - private T[] computeCm3(final Field field, final T[][] am3, final double t) { - // Eq. 51 - final T[] cm3 = MathArrays.buildArray(field, 49); + private T[] computeCm3(final double[] flattenFm3, final T azr, final T[] scT) { + + // interpolation coefficients for effective spot number + final T azr01 = azr.multiply(0.01); + final T omazr01 = azr01.negate().add(1); + + // Eq. 44 and Eq. 51 merged into one loop + final T[] cm3 = MathArrays.buildArray(azr.getField(), 49); + int index = 0; for (int i = 0; i < cm3.length; i++) { - T sum = field.getZero(); - for (int k = 0; k < 4; k++) { - final SinCos sc = FastMath.sinCos((k + 1) * t); - sum = sum.add(am3[i][2 * k + 1].multiply(sc.sin()).add(am3[i][2 * (k + 1)].multiply(sc.cos()))); - } - cm3[i] = am3[i][0].add(sum); + cm3[i] = omazr01.multiply(flattenFm3[index ]).add(azr01.multiply(flattenFm3[index + 1])). + add(omazr01.multiply(flattenFm3[index + 2]).add(azr01.multiply(flattenFm3[index + 3])).multiply(scT[ 0])). + add(omazr01.multiply(flattenFm3[index + 4]).add(azr01.multiply(flattenFm3[index + 5])).multiply(scT[ 1])). + add(omazr01.multiply(flattenFm3[index + 6]).add(azr01.multiply(flattenFm3[index + 7])).multiply(scT[ 2])). + add(omazr01.multiply(flattenFm3[index + 8]).add(azr01.multiply(flattenFm3[index + 9])).multiply(scT[ 3])). + add(omazr01.multiply(flattenFm3[index + 10]).add(azr01.multiply(flattenFm3[index + 11])).multiply(scT[ 4])). + add(omazr01.multiply(flattenFm3[index + 12]).add(azr01.multiply(flattenFm3[index + 13])).multiply(scT[ 5])). + add(omazr01.multiply(flattenFm3[index + 14]).add(azr01.multiply(flattenFm3[index + 15])).multiply(scT[ 6])). + add(omazr01.multiply(flattenFm3[index + 16]).add(azr01.multiply(flattenFm3[index + 17])).multiply(scT[ 7])); + index += 18; } return cm3; } /** * This method computes the F2 layer critical frequency. - * @param field field of the elements * @param modip modified DIP latitude, in degrees * @param cf2 Fourier time series for foF2 * @param latitude latitude in radians - * @param longitude longitude in radians + * @param scL sines/cosines array of longitude argument * @return the F2 layer critical frequency, MHz */ - private T computefoF2(final Field field, final T modip, final T[] cf2, - final T latitude, final T longitude) { - - // One - final T one = field.getOne(); + private T computefoF2(final T modip, final T[] cf2, + final T latitude, final T[] scL) { // Legendre grades (Eq. 63) final int[] q = new int[] { 12, 12, 9, 5, 2, 1, 1, 1, 1 }; - // Array for geographic terms - final T[] g = MathArrays.buildArray(field, cf2.length); - g[0] = one; + T frequency = cf2[0]; // MODIP coefficients Eq. 57 final T sinMODIP = FastMath.sin(FastMath.toRadians(modip)); - final T[] m = MathArrays.buildArray(field, 12); - m[0] = one; + final T[] m = MathArrays.buildArray(latitude.getField(), 12); + m[0] = latitude.getField().getOne(); for (int i = 1; i < q[0]; i++) { m[i] = sinMODIP.multiply(m[i - 1]); - g[i] = m[i]; - } - - // Latitude coefficients (Eq. 58) - final T cosLat = FastMath.cos(latitude); - final T[] p = MathArrays.buildArray(field, 8); - p[0] = cosLat; - for (int n = 2; n < 9; n++) { - p[n - 1] = cosLat.multiply(p[n - 2]); + frequency = frequency.add(m[i].multiply(cf2[i])); } // latitude and longitude terms int index = 12; + final T cosLat1 = FastMath.cos(latitude); + T cosLatI = cosLat1; for (int i = 1; i < q.length; i++) { - final FieldSinCos sc = FastMath.sinCos(longitude.multiply(i)); + final T c = cosLatI.multiply(scL[2 * i - 1]); + final T s = cosLatI.multiply(scL[2 * i - 2]); for (int j = 0; j < q[i]; j++) { - g[index++] = m[j].multiply(p[i - 1]).multiply(sc.cos()); - g[index++] = m[j].multiply(p[i - 1]).multiply(sc.sin()); + frequency = frequency.add(m[j].multiply(c).multiply(cf2[index++])); + frequency = frequency.add(m[j].multiply(s).multiply(cf2[index++])); } + cosLatI = cosLatI.multiply(cosLat1); } - // Compute foF2 by linear combination - final T frequency = one.linearCombination(g, cf2); return frequency; + } /** * This method computes the Maximum Usable Frequency factor. - * @param field field of the elements * @param modip modified DIP latitude, in degrees * @param cm3 Fourier time series for M(3000)F2 * @param latitude latitude in radians - * @param longitude longitude in radians + * @param scL sines/cosines array of longitude argument * @return the Maximum Usable Frequency factor */ - private T computeMF2(final Field field, final T modip, final T[] cm3, - final T latitude, final T longitude) { + private T computeMF2(final T modip, final T[] cm3, + final T latitude, final T[] scL) { - // One - final T one = field.getOne(); // Legendre grades (Eq. 71) final int[] r = new int[] { 7, 8, 6, 3, 2, 1, 1 }; - // Array for geographic terms - final T[] g = MathArrays.buildArray(field, cm3.length); - g[0] = one; + T m3000 = cm3[0]; // MODIP coefficients Eq. 57 final T sinMODIP = FastMath.sin(FastMath.toRadians(modip)); - final T[] m = MathArrays.buildArray(field, 12); - m[0] = one; + final T[] m = MathArrays.buildArray(latitude.getField(), 12); + m[0] = latitude.getField().getOne(); for (int i = 1; i < 12; i++) { m[i] = sinMODIP.multiply(m[i - 1]); if (i < 7) { - g[i] = m[i]; + m3000 = m3000.add(m[i].multiply(cm3[i])); } } - // Latitude coefficients (Eq. 58) - final T cosLat = FastMath.cos(latitude); - final T[] p = MathArrays.buildArray(field, 8); - p[0] = cosLat; - for (int n = 2; n < 9; n++) { - p[n - 1] = cosLat.multiply(p[n - 2]); - } - // latitude and longitude terms int index = 7; + final T cosLat1 = FastMath.cos(latitude); + T cosLatI = cosLat1; for (int i = 1; i < r.length; i++) { - final FieldSinCos sc = FastMath.sinCos(longitude.multiply(i)); + final T c = cosLatI.multiply(scL[2 * i - 1]); + final T s = cosLatI.multiply(scL[2 * i - 2]); for (int j = 0; j < r[i]; j++) { - g[index++] = m[j].multiply(p[i - 1]).multiply(sc.cos()); - g[index++] = m[j].multiply(p[i - 1]).multiply(sc.sin()); + m3000 = m3000.add(m[j].multiply(c).multiply(cm3[index++])); + m3000 = m3000.add(m[j].multiply(s).multiply(cm3[index++])); } + cosLatI = cosLatI.multiply(cosLat1); // Eq. 58 } - // Compute m3000 by linear combination - final T m3000 = one.linearCombination(g, cm3); return m3000; + } /** @@ -618,13 +629,12 @@ private T computeMF2(final Field field, final T modip, final T[] cm3, * This computation performs the algorithm exposed in Annex F * of the reference document. *

- * @param field field of the elements * @param foE the E layer critical frequency, MHz * @return the F1 layer critical frequency, MHz * @param foF2 the F2 layer critical frequency, MHz */ - private T computefoF1(final Field field, final T foE, final T foF2) { - final T zero = field.getZero(); + private T computefoF1(final T foE, final T foF2) { + final T zero = foE.getField().getZero(); final T temp = join(foE.multiply(1.4), zero, zero.newInstance(1000.0), foE.subtract(2.0)); final T temp2 = join(zero, temp, zero.newInstance(1000.0), foE.subtract(temp)); final T value = join(temp2, temp2.multiply(0.85), zero.newInstance(60.0), foF2.multiply(0.85).subtract(temp2)); @@ -645,18 +655,17 @@ private T computefoF1(final Field field, final T foE, final T foF2) { *
  • double[2] = A3 → E layer amplitude * *

    - * @param field field of the elements * @param nmE E layer maximum density in 10^11 m-3 * @param nmF1 F1 layer maximum density in 10^11 m-3 * @param foF1 F1 layer critical frequency in MHz * @return a three components array containing the layer amplitudes */ - private T[] computeLayerAmplitudes(final Field field, final T nmE, final T nmF1, final T foF1) { + private T[] computeLayerAmplitudes(final T nmE, final T nmF1, final T foF1) { // Zero - final T zero = field.getZero(); + final T zero = nmE.getField().getZero(); // Initialize array - final T[] amplitude = MathArrays.buildArray(field, 3); + final T[] amplitude = MathArrays.buildArray(nmE.getField(), 3); // F2 layer amplitude (Eq. 90) final T a1 = nmF2.multiply(4.0); @@ -671,7 +680,7 @@ private T[] computeLayerAmplitudes(final Field field, final T nmE, final T nm T a3a = nmE.multiply(4.0); for (int i = 0; i < 5; i++) { a2a = nmF1.subtract(epst(a1, hmF2, b2Bot, hmF1)).subtract(epst(a3a, hmE, beTop, hmF1)).multiply(4.0); - a2a = join(a2a, nmF1.multiply(0.8), field.getOne(), a2a.subtract(nmF1.multiply(0.8))); + a2a = join(a2a, nmF1.multiply(0.8), nmE.getField().getOne(), a2a.subtract(nmF1.multiply(0.8))); a3a = nmE.subtract(epst(a2a, hmF1, b1Bot, hmE)).subtract(epst(a1, hmF2, b2Bot, hmE)).multiply(4.0); } amplitude[1] = a2a; @@ -684,15 +693,14 @@ private T[] computeLayerAmplitudes(final Field field, final T nmE, final T nm /** * This method computes the topside thickness parameter H0. * - * @param field field of the elements * @param month current month * @param azr effective sunspot number * @return H0 in km */ - private T computeH0(final Field field, final int month, final T azr) { + private T computeH0(final int month, final T azr) { // One - final T one = field.getOne(); + final T one = azr.getField().getOne(); // Auxiliary parameter ka (Eq. 99 and 100) final T ka; @@ -717,8 +725,8 @@ private T computeH0(final Field field, final int month, final T azr) { final T v = x.multiply(0.041163).subtract(0.183981).multiply(x).add(1.424472); // Topside thickness parameter (Eq. 106) - final T h = hA.divide(v); - return h; + return hA.divide(v); + } /** @@ -769,9 +777,8 @@ private T interpolate(final T z1, final T z2, final T a1 = g2.multiply(9.0).subtract(g4); final T a2 = g3.subtract(g1); final T a3 = g4.subtract(g2); - final T zx = delta.multiply(a3).add(a2).multiply(delta).add(a1).multiply(delta).add(a0).multiply(0.0625); + return delta.multiply(a3).add(a2).multiply(delta).add(a1).multiply(delta).add(a0).multiply(0.0625); - return zx; } /** @@ -809,8 +816,8 @@ private T epst(final T x, final T y, final T z, final T w) { final T ex = clipExp(w.subtract(y).divide(z)); final T opex = ex.add(1.0); - final T epst = x.multiply(ex).divide(opex.multiply(opex)); - return epst; + return x.multiply(ex).divide(opex.multiply(opex)); + } } diff --git a/src/main/java/org/orekit/models/earth/ionosphere/NeQuickModel.java b/src/main/java/org/orekit/models/earth/ionosphere/NeQuickModel.java index 1af899a1c9..1824457b5a 100644 --- a/src/main/java/org/orekit/models/earth/ionosphere/NeQuickModel.java +++ b/src/main/java/org/orekit/models/earth/ionosphere/NeQuickModel.java @@ -87,11 +87,11 @@ public class NeQuickModel implements IonosphericModel { /** Month used for loading CCIR coefficients. */ private int month; - /** F2 coefficients used by the F2 layer. */ - private double[][][] f2; + /** F2 coefficients used by the F2 layer (flatten array for cache efficiency). */ + private double[] flattenF2; - /** Fm3 coefficients used by the F2 layer. */ - private double[][][] fm3; + /** Fm3 coefficients used by the F2 layer(flatten array for cache efficiency). */ + private double[] flattenFm3; /** UTC time scale. */ private final TimeScale utc; @@ -118,9 +118,9 @@ public NeQuickModel(final double[] alpha) { public NeQuickModel(final double[] alpha, final TimeScale utc) { // F2 layer values - this.month = 0; - this.f2 = null; - this.fm3 = null; + this.month = 0; + this.flattenF2 = null; + this.flattenFm3 = null; // Read modip grid final MODIPLoader parser = new MODIPLoader(); parser.loadMODIPGrid(); @@ -309,7 +309,7 @@ private double stecIntegration(final Segment seg, final DateTimeComponents dateT // Compute electron density double density = 0.0; for (int i = 0; i < heightS.length; i++) { - final NeQuickParameters parameters = new NeQuickParameters(dateTime, f2, fm3, + final NeQuickParameters parameters = new NeQuickParameters(dateTime, flattenF2, flattenFm3, latitudeS[i], longitudeS[i], alpha, stModip); density += electronDensity(heightS[i], parameters); @@ -327,8 +327,8 @@ private double stecIntegration(final Segment seg, final DateTimeComponents dateT * @return result of the integration */ private > T stecIntegration(final Field field, - final FieldSegment seg, - final DateTimeComponents dateTime) { + final FieldSegment seg, + final DateTimeComponents dateTime) { // Integration points final T[] heightS = seg.getHeights(); final T[] latitudeS = seg.getLatitudes(); @@ -337,7 +337,7 @@ private > T stecIntegration(final Field fie // Compute electron density T density = field.getZero(); for (int i = 0; i < heightS.length; i++) { - final FieldNeQuickParameters parameters = new FieldNeQuickParameters<>(field, dateTime, f2, fm3, + final FieldNeQuickParameters parameters = new FieldNeQuickParameters<>(dateTime, flattenF2, flattenFm3, latitudeS[i], longitudeS[i], alpha, stModip); density = density.add(electronDensity(field, heightS[i], parameters)); @@ -609,7 +609,7 @@ private void loadsIfNeeded(final DateComponents date) { final int currentMonth = date.getMonth(); // Check if date have changed or if f2 and fm3 arrays are null - if (currentMonth != month || f2 == null || fm3 == null) { + if (currentMonth != month || flattenF2 == null || flattenFm3 == null) { this.month = currentMonth; // Read file @@ -617,11 +617,32 @@ private void loadsIfNeeded(final DateComponents date) { loader.loadCCIRCoefficients(date); // Update arrays - this.f2 = loader.getF2(); - this.fm3 = loader.getFm3(); + this.flattenF2 = flatten(loader.getF2()); + this.flattenFm3 = flatten(loader.getFm3()); } } + /** Flatten a 3-dimensions array. + *

    + * This method convert 3-dimensions arrays into 1-dimension arrays + * optimized to avoid cache misses when looping over all elements. + *

    + * @param original original array a[i][j][k] + * @return flatten array, for embedded loops on j (outer), k (intermediate), i (inner) + */ + private double[] flatten(final double[][][] original) { + final double[] flatten = new double[original.length * original[0].length * original[0][0].length]; + int index = 0; + for (int j = 0; j < original[0].length; j++) { + for (int k = 0; k < original[0][0].length; k++) { + for (final double[][] doubles : original) { + flatten[index++] = doubles[j][k]; + } + } + } + return flatten; + } + /** * A clipped exponential function. *

    @@ -746,7 +767,7 @@ public void loadData(final InputStream input, final String name) line = line.trim(); // Read grid data - if (line.length() > 0) { + if (!line.isEmpty()) { final String[] modip_line = SEPARATOR.split(line); for (int column = 0; column < modip_line.length; column++) { array[lineNumber - 1][column] = Double.parseDouble(modip_line[column]); @@ -891,9 +912,9 @@ public void loadData(final InputStream input, final String name) line = line.trim(); // Read grid data - if (line.length() > 0) { + if (!line.isEmpty()) { final String[] ccir_line = SEPARATOR.split(line); - for (int i = 0; i < ccir_line.length; i++) { + for (final String field : ccir_line) { if (index < NUMBER_F2_COEFFICIENTS) { // Parse F2 coefficients @@ -905,7 +926,7 @@ public void loadData(final InputStream input, final String name) currentColumnF2 = 0; currentRowF2++; } - f2Temp[currentRowF2][currentColumnF2][currentDepthF2++] = Double.parseDouble(ccir_line[i]); + f2Temp[currentRowF2][currentColumnF2][currentDepthF2++] = Double.parseDouble(field); index++; } else { // Parse Fm3 coefficients @@ -917,7 +938,7 @@ public void loadData(final InputStream input, final String name) currentColumnFm3 = 0; currentRowFm3++; } - fm3Temp[currentRowFm3][currentColumnFm3][currentDepthFm3++] = Double.parseDouble(ccir_line[i]); + fm3Temp[currentRowFm3][currentColumnFm3][currentDepthFm3++] = Double.parseDouble(field); index++; } @@ -995,20 +1016,18 @@ private static class Ray { final double cosD = scLatRec.sin() * scLatSat.sin() + scLatRec.cos() * scLatSat.cos() * scLon21.cos(); final double sinD = FastMath.sqrt(1.0 - cosD * cosD); - final double z = FastMath.atan2(sinD, cosD - (r1 / r2)); + final double z = FastMath.atan2(sinD, cosD - (r1 / r2)); + final SinCos scZ = FastMath.sinCos(z); // Ray-perigee computation in meters (Eq. 156) - this.rp = r1 * FastMath.sin(z); + this.rp = r1 * scZ.sin(); // Ray-perigee latitude and longitude if (FastMath.abs(FastMath.abs(lat1) - 0.5 * FastMath.PI) < THRESHOLD) { + // receiver is almost at North or South pole // Ray-perigee latitude (Eq. 157) - if (lat1 < 0) { - this.latP = -z; - } else { - this.latP = z; - } + this.latP = FastMath.copySign(z, lat1); // Ray-perigee longitude (Eq. 164) if (z < 0) { @@ -1017,20 +1036,24 @@ private static class Ray { this.lonP = lon2 + FastMath.PI; } + } else if (FastMath.abs(scZ.sin()) < THRESHOLD) { + // satellite is almost on receiver zenith + + this.latP = recP.getLatitude(); + this.lonP = recP.getLongitude(); + } else { // Ray-perigee latitude (Eq. 158 to 163) - final double deltaP = 0.5 * FastMath.PI - z; - final SinCos scDeltaP = FastMath.sinCos(deltaP); final double sinAz = scLon21.sin() * scLatSat.cos() / sinD; final double cosAz = (scLatSat.sin() - cosD * scLatRec.sin()) / (sinD * scLatRec.cos()); - final double sinLatP = scLatRec.sin() * scDeltaP.cos() - scLatRec.cos() * scDeltaP.sin() * cosAz; + final double sinLatP = scLatRec.sin() * scZ.sin() - scLatRec.cos() * scZ.cos() * cosAz; final double cosLatP = FastMath.sqrt(1.0 - sinLatP * sinLatP); this.latP = FastMath.atan2(sinLatP, cosLatP); // Ray-perigee longitude (Eq. 165 to 167) - final double sinLonP = -sinAz * scDeltaP.sin() / cosLatP; - final double cosLonP = (scDeltaP.cos() - scLatRec.sin() * sinLatP) / (scLatRec.cos() * cosLatP); + final double sinLonP = -sinAz * scZ.cos() / cosLatP; + final double cosLonP = (scZ.sin() - scLatRec.sin() * sinLatP) / (scLatRec.cos() * cosLatP); this.lonP = FastMath.atan2(sinLonP, cosLonP) + lon1; } @@ -1038,19 +1061,15 @@ private static class Ray { // Sine and cosine of ray-perigee latitude this.scLatP = FastMath.sinCos(latP); - final SinCos scLon = FastMath.sinCos(lon2 - lonP); - // Sine and cosine of azimuth of satellite as seen from ray-perigee - final double psi = greatCircleAngle(scLatSat, scLon); - final SinCos scPsi = FastMath.sinCos(psi); - if (FastMath.abs(FastMath.abs(latP) - 0.5 * FastMath.PI) < THRESHOLD) { + if (FastMath.abs(FastMath.abs(latP) - 0.5 * FastMath.PI) < THRESHOLD || + FastMath.abs(scZ.sin()) < THRESHOLD) { // Eq. 172 and 173 this.sinAzP = 0.0; - if (latP < 0.0) { - this.cosAzP = 1; - } else { - this.cosAzP = -1; - } + this.cosAzP = -FastMath.copySign(1, latP); } else { + final SinCos scLon = FastMath.sinCos(lon2 - lonP); + // Sine and cosine of azimuth of satellite as seen from ray-perigee + final SinCos scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon)); // Eq. 174 and 175 this.sinAzP = scLatSat.cos() * scLon.sin() / scPsi.sin(); this.cosAzP = (scLatSat.sin() - scLatP.sin() * scPsi.cos()) / (scLatP.cos() * scPsi.sin()); @@ -1191,25 +1210,23 @@ private static class FieldRay > { final T lon2 = satP.getLongitude(); final FieldSinCos scLatSat = FastMath.sinCos(lat2); final FieldSinCos scLatRec = FastMath.sinCos(lat1); + final FieldSinCos scLon21 = FastMath.sinCos(lon2.subtract(lon1)); // Zenith angle computation (Eq. 153 to 155) final T cosD = scLatRec.sin().multiply(scLatSat.sin()). - add(scLatRec.cos().multiply(scLatSat.cos()).multiply(FastMath.cos(lon2.subtract(lon1)))); + add(scLatRec.cos().multiply(scLatSat.cos()).multiply(scLon21.cos())); final T sinD = FastMath.sqrt(cosD.multiply(cosD).negate().add(1.0)); final T z = FastMath.atan2(sinD, cosD.subtract(r1.divide(r2))); + final FieldSinCos scZ = FastMath.sinCos(z); // Ray-perigee computation in meters (Eq. 156) - this.rp = r1.multiply(FastMath.sin(z)); + this.rp = r1.multiply(scZ.sin()); // Ray-perigee latitude and longitude if (FastMath.abs(FastMath.abs(lat1).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD) { // Ray-perigee latitude (Eq. 157) - if (lat1.getReal() < 0) { - this.latP = z.negate(); - } else { - this.latP = z; - } + this.latP = FastMath.copySign(z, lat1); // Ray-perigee longitude (Eq. 164) if (z.getReal() < 0) { @@ -1218,20 +1235,24 @@ private static class FieldRay > { this.lonP = lon2.add(pi); } + } else if (FastMath.abs(scZ.sin().getReal()) < THRESHOLD) { + // satellite is almost on receiver zenith + + this.latP = recP.getLatitude(); + this.lonP = recP.getLongitude(); + } else { // Ray-perigee latitude (Eq. 158 to 163) - final T deltaP = z.negate().add(pi.multiply(0.5)); - final FieldSinCos scDeltaP = FastMath.sinCos(deltaP); final T sinAz = FastMath.sin(lon2.subtract(lon1)).multiply(scLatSat.cos()).divide(sinD); final T cosAz = scLatSat.sin().subtract(cosD.multiply(scLatRec.sin())).divide(sinD.multiply(scLatRec.cos())); - final T sinLatP = scLatRec.sin().multiply(scDeltaP.cos()).subtract(scLatRec.cos().multiply(scDeltaP.sin()).multiply(cosAz)); + final T sinLatP = scLatRec.sin().multiply(scZ.sin()).subtract(scLatRec.cos().multiply(scZ.cos()).multiply(cosAz)); final T cosLatP = FastMath.sqrt(sinLatP.multiply(sinLatP).negate().add(1.0)); this.latP = FastMath.atan2(sinLatP, cosLatP); // Ray-perigee longitude (Eq. 165 to 167) - final T sinLonP = sinAz.negate().multiply(scDeltaP.sin()).divide(cosLatP); - final T cosLonP = scDeltaP.cos().subtract(scLatRec.sin().multiply(sinLatP)).divide(scLatRec.cos().multiply(cosLatP)); + final T sinLonP = sinAz.negate().multiply(scZ.cos()).divide(cosLatP); + final T cosLonP = scZ.sin().subtract(scLatRec.sin().multiply(sinLatP)).divide(scLatRec.cos().multiply(cosLatP)); this.lonP = FastMath.atan2(sinLonP, cosLonP).add(lon1); } @@ -1239,19 +1260,15 @@ private static class FieldRay > { // Sine and cosine of ray-perigee latitude this.scLatP = FastMath.sinCos(latP); - final FieldSinCos scLon = FastMath.sinCos(lon2.subtract(lonP)); - // Sine and cosie of azimuth of satellite as seen from ray-perigee - final T psi = greatCircleAngle(scLatSat, scLon); - final FieldSinCos scPsi = FastMath.sinCos(psi); - if (FastMath.abs(FastMath.abs(latP).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD) { + if (FastMath.abs(FastMath.abs(latP).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD || + FastMath.abs(scZ.sin().getReal()) < THRESHOLD) { // Eq. 172 and 173 this.sinAzP = field.getZero(); - if (latP.getReal() < 0.0) { - this.cosAzP = field.getOne(); - } else { - this.cosAzP = field.getOne().negate(); - } + this.cosAzP = FastMath.copySign(field.getOne(), latP).negate(); } else { + final FieldSinCos scLon = FastMath.sinCos(lon2.subtract(lonP)); + // Sine and cosine of azimuth of satellite as seen from ray-perigee + final FieldSinCos scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon)); // Eq. 174 and 175 this.sinAzP = scLatSat.cos().multiply(scLon.sin()).divide(scPsi.sin()); this.cosAzP = scLatSat.sin().subtract(scLatP.sin().multiply(scPsi.cos())).divide(scLatP.cos().multiply(scPsi.sin())); @@ -1342,6 +1359,9 @@ private T greatCircleAngle(final FieldSinCos scLat, final FieldSinCos scLo /** Performs the computation of the coordinates along the integration path. */ private static class Segment { + /** Threshold for zenith segment. */ + private static final double THRESHOLD = 1.0e-3; + /** Latitudes [rad]. */ private final double[] latitudes; @@ -1368,7 +1388,7 @@ private static class Segment { this.deltaN = (s2 - s1) / n; // Segments - final double[] s = getSegments(n, s1, s2); + final double[] s = getSegments(n, s1); // Useful parameters final double rp = ray.getRadius(); @@ -1383,20 +1403,26 @@ private static class Segment { // Heights (Eq. 178) heights[i] = FastMath.sqrt(s[i] * s[i] + rp * rp) - RE; - // Great circle parameters (Eq. 179 to 181) - final double tanDs = s[i] / rp; - final double cosDs = 1.0 / FastMath.sqrt(1.0 + tanDs * tanDs); - final double sinDs = tanDs * cosDs; - - // Latitude (Eq. 182 to 183) - final double sinLatS = scLatP.sin() * cosDs + scLatP.cos() * sinDs * ray.getCosineAz(); - final double cosLatS = FastMath.sqrt(1.0 - sinLatS * sinLatS); - latitudes[i] = FastMath.atan2(sinLatS, cosLatS); - - // Longitude (Eq. 184 to 187) - final double sinLonS = sinDs * ray.getSineAz() * scLatP.cos(); - final double cosLonS = cosDs - scLatP.sin() * sinLatS; - longitudes[i] = FastMath.atan2(sinLonS, cosLonS) + ray.getLongitude(); + if (rp < THRESHOLD) { + // zenith segment + latitudes[i] = ray.getLatitude(); + longitudes[i] = ray.getLongitude(); + } else { + // Great circle parameters (Eq. 179 to 181) + final double tanDs = s[i] / rp; + final double cosDs = 1.0 / FastMath.sqrt(1.0 + tanDs * tanDs); + final double sinDs = tanDs * cosDs; + + // Latitude (Eq. 182 to 183) + final double sinLatS = scLatP.sin() * cosDs + scLatP.cos() * sinDs * ray.getCosineAz(); + final double cosLatS = FastMath.sqrt(1.0 - sinLatS * sinLatS); + latitudes[i] = FastMath.atan2(sinLatS, cosLatS); + + // Longitude (Eq. 184 to 187) + final double sinLonS = sinDs * ray.getSineAz() * scLatP.cos(); + final double cosLonS = cosDs - scLatP.sin() * sinLatS; + longitudes[i] = FastMath.atan2(sinLonS, cosLonS) + ray.getLongitude(); + } } } @@ -1404,10 +1430,9 @@ private static class Segment { * Computes the distance of a point from the ray-perigee. * @param n number of points used for the integration * @param s1 lower boundary - * @param s2 upper boundary * @return the distance of a point from the ray-perigee in km */ - private double[] getSegments(final int n, final double s1, final double s2) { + private double[] getSegments(final int n, final double s1) { // Eq. 196 final double g = 0.5773502691896 * deltaN; // Eq. 197 @@ -1460,6 +1485,9 @@ public double getInterval() { /** Performs the computation of the coordinates along the integration path. */ private static class FieldSegment > { + /** Threshold for zenith segment. */ + private static final double THRESHOLD = 1.0e-3; + /** Latitudes [rad]. */ private final T[] latitudes; @@ -1487,7 +1515,7 @@ private static class FieldSegment > { this.deltaN = s2.subtract(s1).divide(n); // Segments - final T[] s = getSegments(field, n, s1, s2); + final T[] s = getSegments(field, n, s1); // Useful parameters final T rp = ray.getRadius(); @@ -1502,20 +1530,26 @@ private static class FieldSegment > { // Heights (Eq. 178) heights[i] = FastMath.sqrt(s[i].multiply(s[i]).add(rp.multiply(rp))).subtract(RE); - // Great circle parameters (Eq. 179 to 181) - final T tanDs = s[i].divide(rp); - final T cosDs = FastMath.sqrt(tanDs.multiply(tanDs).add(1.0)).reciprocal(); - final T sinDs = tanDs.multiply(cosDs); - - // Latitude (Eq. 182 to 183) - final T sinLatS = scLatP.sin().multiply(cosDs).add(scLatP.cos().multiply(sinDs).multiply(ray.getCosineAz())); - final T cosLatS = FastMath.sqrt(sinLatS.multiply(sinLatS).negate().add(1.0)); - latitudes[i] = FastMath.atan2(sinLatS, cosLatS); - - // Longitude (Eq. 184 to 187) - final T sinLonS = sinDs.multiply(ray.getSineAz()).multiply(scLatP.cos()); - final T cosLonS = cosDs.subtract(scLatP.sin().multiply(sinLatS)); - longitudes[i] = FastMath.atan2(sinLonS, cosLonS).add(ray.getLongitude()); + if (rp.getReal() < THRESHOLD) { + // zenith segment + latitudes[i] = ray.getLatitude(); + longitudes[i] = ray.getLongitude(); + } else { + // Great circle parameters (Eq. 179 to 181) + final T tanDs = s[i].divide(rp); + final T cosDs = FastMath.sqrt(tanDs.multiply(tanDs).add(1.0)).reciprocal(); + final T sinDs = tanDs.multiply(cosDs); + + // Latitude (Eq. 182 to 183) + final T sinLatS = scLatP.sin().multiply(cosDs).add(scLatP.cos().multiply(sinDs).multiply(ray.getCosineAz())); + final T cosLatS = FastMath.sqrt(sinLatS.multiply(sinLatS).negate().add(1.0)); + latitudes[i] = FastMath.atan2(sinLatS, cosLatS); + + // Longitude (Eq. 184 to 187) + final T sinLonS = sinDs.multiply(ray.getSineAz()).multiply(scLatP.cos()); + final T cosLonS = cosDs.subtract(scLatP.sin().multiply(sinLatS)); + longitudes[i] = FastMath.atan2(sinLonS, cosLonS).add(ray.getLongitude()); + } } } @@ -1524,10 +1558,9 @@ private static class FieldSegment > { * @param field field of the elements * @param n number of points used for the integration * @param s1 lower boundary - * @param s2 upper boundary * @return the distance of a point from the ray-perigee in km */ - private T[] getSegments(final Field field, final int n, final T s1, final T s2) { + private T[] getSegments(final Field field, final int n, final T s1) { // Eq. 196 final T g = deltaN.multiply(0.5773502691896); // Eq. 197 diff --git a/src/main/java/org/orekit/models/earth/ionosphere/NeQuickParameters.java b/src/main/java/org/orekit/models/earth/ionosphere/NeQuickParameters.java index 1ad4b57e68..97042983e8 100644 --- a/src/main/java/org/orekit/models/earth/ionosphere/NeQuickParameters.java +++ b/src/main/java/org/orekit/models/earth/ionosphere/NeQuickParameters.java @@ -17,7 +17,6 @@ package org.orekit.models.earth.ionosphere; import org.hipparchus.util.FastMath; -import org.hipparchus.util.MathArrays; import org.hipparchus.util.SinCos; import org.orekit.time.DateComponents; import org.orekit.time.DateTimeComponents; @@ -74,15 +73,16 @@ class NeQuickParameters { /** * Build a new instance. * @param dateTime current date time components - * @param f2 F2 coefficients used by the F2 layer - * @param fm3 Fm3 coefficients used by the F2 layer + * @param flattenF2 F2 coefficients used by the F2 layer (flatten array) + * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array) * @param latitude latitude of a point along the integration path, in radians * @param longitude longitude of a point along the integration path, in radians * @param alpha effective ionisation level coefficients * @param modipGrip modip grid */ - NeQuickParameters(final DateTimeComponents dateTime, final double[][][] f2, - final double[][][] fm3, final double latitude, final double longitude, + NeQuickParameters(final DateTimeComponents dateTime, + final double[] flattenF2, final double[] flattenFm3, + final double latitude, final double longitude, final double[] alpha, final double[][] modipGrip) { // MODIP in degrees @@ -99,23 +99,6 @@ class NeQuickParameters { // Effective solar zenith angle in radians final double xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude); - // Coefficients for F2 layer parameters - // Compute the array of interpolated coefficients for foF2 (Eq. 44) - final double[][] af2 = new double[76][13]; - for (int j = 0; j < 76; j++) { - for (int k = 0; k < 13; k++ ) { - af2[j][k] = f2[0][j][k] * (1.0 - (azr * 0.01)) + f2[1][j][k] * (azr * 0.01); - } - } - - // Compute the array of interpolated coefficients for M(3000)F2 (Eq. 46) - final double[][] am3 = new double[49][9]; - for (int j = 0; j < 49; j++) { - for (int k = 0; k < 9; k++ ) { - am3[j][k] = fm3[0][j][k] * (1.0 - (azr * 0.01)) + fm3[1][j][k] * (azr * 0.01); - } - } - // E layer maximum density height in km (Eq. 78) this.hmE = 120.0; // E layer critical frequency in MHz @@ -126,12 +109,14 @@ class NeQuickParameters { // Time argument (Eq. 49) final double t = FastMath.toRadians(15 * hours) - FastMath.PI; // Compute Fourier time series for foF2 and M(3000)F2 - final double[] cf2 = computeCF2(af2, t); - final double[] cm3 = computeCm3(am3, t); + final double[] scT = sinCos(t, 6); + final double[] cf2 = computeCF2(flattenF2, azr, scT); + final double[] cm3 = computeCm3(flattenFm3, azr, scT); // F2 layer critical frequency in MHz - final double foF2 = computefoF2(modip, cf2, latitude, longitude); + final double[] scL = sinCos(longitude, 8); + final double foF2 = computefoF2(modip, cf2, latitude, scL); // Maximum Usable Frequency factor - final double mF2 = computeMF2(modip, cm3, latitude, longitude); + final double mF2 = computeMF2(modip, cm3, latitude, scL); // F2 layer maximum density in 10^11 m-3 this.nmF2 = 0.124 * foF2 * foF2; // F2 layer maximum density height in km @@ -309,9 +294,8 @@ private double computeMODIP(final double lat, final double lon, final double[][] final double y = b - bF; // MODIP (Ref Eq. 16) - final double modip = interpolate(z1, z2, z3, z4, y); + return interpolate(z1, z2, z3, z4, y); - return modip; } /** @@ -404,8 +388,8 @@ private double computefoE(final int month, final double az, final double seasp = seas * ((ee - 1.0) / (ee + 1.0)); // Critical frequency (Eq. 35) final double coef = 1.112 - 0.019 * seasp; - final double foE = FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49); - return foE; + return FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49); + } /** @@ -429,46 +413,96 @@ private double computehmF2(final double foE, final double foF2, final double mF2 // hmF2 Eq. 80 final double mF2Sq = mF2 * mF2; final double temp = FastMath.sqrt((0.0196 * mF2Sq + 1) / (1.2967 * mF2Sq - 1.0)); - final double height = ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0; - return height; + return ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0; + + } + + /** Compute sines and cosines. + * @param a argument + * @param n number of terms + * @return sin(a), cos(a), sin(2a), cos(2a) … sin(n a), cos(n a) array + * @since 12.1.3 + */ + private double[] sinCos(final double a, final int n) { + + final SinCos sc0 = FastMath.sinCos(a); + SinCos sci = sc0; + final double[] sc = new double[2 * n]; + int isc = 0; + sc[isc++] = sci.sin(); + sc[isc++] = sci.cos(); + for (int i = 1; i < n; i++) { + sci = SinCos.sum(sc0, sci); + sc[isc++] = sci.sin(); + sc[isc++] = sci.cos(); + } + + return sc; + } /** * Computes cf2 coefficients. - * @param af2 interpolated coefficients for foF2 - * @param t time argument + * @param flattenF2 F2 coefficients used by the F2 layer (flatten array) + * @param azr effective sunspot number (Eq. 19) + * @param scT sines/cosines array of time argument * @return the cf2 coefficients array */ - private double[] computeCF2(final double[][] af2, final double t) { - // Eq. 50 + private double[] computeCF2(final double[] flattenF2, final double azr, final double[] scT) { + + // interpolation coefficients for effective spot number + final double azr01 = azr * 0.01; + final double omazr01 = 1 - azr01; + + // Eq. 44 and Eq. 50 merged into one loop final double[] cf2 = new double[76]; + int index = 0; for (int i = 0; i < cf2.length; i++) { - double sum = 0.0; - for (int k = 0; k < 6; k++) { - final SinCos sc = FastMath.sinCos((k + 1) * t); - sum += af2[i][2 * k + 1] * sc.sin() + af2[i][2 * (k + 1)] * sc.cos(); - } - cf2[i] = af2[i][0] + sum; + cf2[i] = omazr01 * flattenF2[index ] + azr01 * flattenF2[index + 1] + + (omazr01 * flattenF2[index + 2] + azr01 * flattenF2[index + 3]) * scT[ 0] + + (omazr01 * flattenF2[index + 4] + azr01 * flattenF2[index + 5]) * scT[ 1] + + (omazr01 * flattenF2[index + 6] + azr01 * flattenF2[index + 7]) * scT[ 2] + + (omazr01 * flattenF2[index + 8] + azr01 * flattenF2[index + 9]) * scT[ 3] + + (omazr01 * flattenF2[index + 10] + azr01 * flattenF2[index + 11]) * scT[ 4] + + (omazr01 * flattenF2[index + 12] + azr01 * flattenF2[index + 13]) * scT[ 5] + + (omazr01 * flattenF2[index + 14] + azr01 * flattenF2[index + 15]) * scT[ 6] + + (omazr01 * flattenF2[index + 16] + azr01 * flattenF2[index + 17]) * scT[ 7] + + (omazr01 * flattenF2[index + 18] + azr01 * flattenF2[index + 19]) * scT[ 8] + + (omazr01 * flattenF2[index + 20] + azr01 * flattenF2[index + 21]) * scT[ 9] + + (omazr01 * flattenF2[index + 22] + azr01 * flattenF2[index + 23]) * scT[10] + + (omazr01 * flattenF2[index + 24] + azr01 * flattenF2[index + 25]) * scT[11]; + index += 26; } return cf2; } /** * Computes Cm3 coefficients. - * @param am3 interpolated coefficients for foF2 - * @param t time argument + * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array) + * @param azr effective sunspot number (Eq. 19) + * @param scT sines/cosines array of time argument * @return the Cm3 coefficients array */ - private double[] computeCm3(final double[][] am3, final double t) { - // Eq. 51 + private double[] computeCm3(final double[] flattenFm3, final double azr, final double[] scT) { + + // interpolation coefficients for effective spot number + final double azr01 = azr * 0.01; + final double omazr01 = 1 - azr01; + + // Eq. 44 and Eq. 51 merged into one loop final double[] cm3 = new double[49]; + int index = 0; for (int i = 0; i < cm3.length; i++) { - double sum = 0.0; - for (int k = 0; k < 4; k++) { - final SinCos sc = FastMath.sinCos((k + 1) * t); - sum += am3[i][2 * k + 1] * sc.sin() + am3[i][2 * (k + 1)] * sc.cos(); - } - cm3[i] = am3[i][0] + sum; + cm3[i] = omazr01 * flattenFm3[index ] + azr01 * flattenFm3[index + 1] + + (omazr01 * flattenFm3[index + 2] + azr01 * flattenFm3[index + 3]) * scT[ 0] + + (omazr01 * flattenFm3[index + 4] + azr01 * flattenFm3[index + 5]) * scT[ 1] + + (omazr01 * flattenFm3[index + 6] + azr01 * flattenFm3[index + 7]) * scT[ 2] + + (omazr01 * flattenFm3[index + 8] + azr01 * flattenFm3[index + 9]) * scT[ 3] + + (omazr01 * flattenFm3[index + 10] + azr01 * flattenFm3[index + 11]) * scT[ 4] + + (omazr01 * flattenFm3[index + 12] + azr01 * flattenFm3[index + 13]) * scT[ 5] + + (omazr01 * flattenFm3[index + 14] + azr01 * flattenFm3[index + 15]) * scT[ 6] + + (omazr01 * flattenFm3[index + 16] + azr01 * flattenFm3[index + 17]) * scT[ 7]; + index += 18; } return cm3; } @@ -478,20 +512,18 @@ private double[] computeCm3(final double[][] am3, final double t) { * @param modip modified DIP latitude, in degrees * @param cf2 Fourier time series for foF2 * @param latitude latitude in radians - * @param longitude longitude in radians + * @param scL sines/cosines array of longitude argument * @return the F2 layer critical frequency, MHz */ private double computefoF2(final double modip, final double[] cf2, - final double latitude, final double longitude) { + final double latitude, final double[] scL) { // Legendre grades (Eq. 63) final int[] q = new int[] { 12, 12, 9, 5, 2, 1, 1, 1, 1 }; - // Array for geographic terms - final double[] g = new double[cf2.length]; - g[0] = 1.0; + double frequency = cf2[0]; // MODIP coefficients Eq. 57 final double sinMODIP = FastMath.sin(FastMath.toRadians(modip)); @@ -499,30 +531,25 @@ private double computefoF2(final double modip, final double[] cf2, m[0] = 1.0; for (int i = 1; i < q[0]; i++) { m[i] = sinMODIP * m[i - 1]; - g[i] = m[i]; - } - - // Latitude coefficients (Eq. 58) - final double cosLat = FastMath.cos(latitude); - final double[] p = new double[8]; - p[0] = cosLat; - for (int n = 2; n < 9; n++) { - p[n - 1] = cosLat * p[n - 2]; + frequency += m[i] * cf2[i]; } // latitude and longitude terms int index = 12; + final double cosLat1 = FastMath.cos(latitude); + double cosLatI = cosLat1; for (int i = 1; i < q.length; i++) { - final SinCos sc = FastMath.sinCos(i * longitude); + final double c = cosLatI * scL[2 * i - 1]; + final double s = cosLatI * scL[2 * i - 2]; for (int j = 0; j < q[i]; j++) { - g[index++] = m[j] * p[i - 1] * sc.cos(); - g[index++] = m[j] * p[i - 1] * sc.sin(); + frequency += m[j] * c * cf2[index++]; + frequency += m[j] * s * cf2[index++]; } + cosLatI *= cosLat1; // Eq. 58 } - // Compute foF2 by linear combination - final double frequency = MathArrays.linearCombination(cf2, g); return frequency; + } /** @@ -530,20 +557,18 @@ private double computefoF2(final double modip, final double[] cf2, * @param modip modified DIP latitude, in degrees * @param cm3 Fourier time series for M(3000)F2 * @param latitude latitude in radians - * @param longitude longitude in radians + * @param scL sines/cosines array of longitude argument * @return the Maximum Usable Frequency factor */ private double computeMF2(final double modip, final double[] cm3, - final double latitude, final double longitude) { + final double latitude, final double[] scL) { // Legendre grades (Eq. 71) final int[] r = new int[] { 7, 8, 6, 3, 2, 1, 1 }; - // Array for geographic terms - final double[] g = new double[cm3.length]; - g[0] = 1.0; + double m3000 = cm3[0]; // MODIP coefficients Eq. 57 final double sinMODIP = FastMath.sin(FastMath.toRadians(modip)); @@ -552,31 +577,26 @@ private double computeMF2(final double modip, final double[] cm3, for (int i = 1; i < 12; i++) { m[i] = sinMODIP * m[i - 1]; if (i < 7) { - g[i] = m[i]; + m3000 += m[i] * cm3[i]; } } - // Latitude coefficients (Eq. 58) - final double cosLat = FastMath.cos(latitude); - final double[] p = new double[8]; - p[0] = cosLat; - for (int n = 2; n < 9; n++) { - p[n - 1] = cosLat * p[n - 2]; - } - // latitude and longitude terms int index = 7; + final double cosLat1 = FastMath.cos(latitude); + double cosLatI = cosLat1; for (int i = 1; i < r.length; i++) { - final SinCos sc = FastMath.sinCos(i * longitude); + final double c = cosLatI * scL[2 * i - 1]; + final double s = cosLatI * scL[2 * i - 2]; for (int j = 0; j < r[i]; j++) { - g[index++] = m[j] * p[i - 1] * sc.cos(); - g[index++] = m[j] * p[i - 1] * sc.sin(); + m3000 += m[j] * c * cm3[index++]; + m3000 += m[j] * s * cm3[index++]; } + cosLatI *= cosLat1; // Eq. 58 } - // Compute m3000 by linear combination - final double m3000 = MathArrays.linearCombination(g, cm3); return m3000; + } /** @@ -674,8 +694,8 @@ private double computeH0(final int month, final double azr) { final double v = (0.041163 * x - 0.183981) * x + 1.424472; // Topside thickness parameter (Eq. 106) - final double h = hA / v; - return h; + return hA / v; + } /** @@ -725,9 +745,8 @@ private double interpolate(final double z1, final double z2, final double a1 = 9.0 * g2 - g4; final double a2 = g3 - g1; final double a3 = g4 - g2; - final double zx = 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3))); + return 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3))); - return zx; } /** @@ -765,8 +784,8 @@ private double epst(final double x, final double y, final double z, final double w) { final double ex = clipExp((w - y) / z); final double opex = 1.0 + ex; - final double epst = x * ex / (opex * opex); - return epst; + return x * ex / (opex * opex); + } } diff --git a/src/main/java/org/orekit/propagation/FieldStateCovariance.java b/src/main/java/org/orekit/propagation/FieldStateCovariance.java index 7001bf8b65..a675e90adf 100644 --- a/src/main/java/org/orekit/propagation/FieldStateCovariance.java +++ b/src/main/java/org/orekit/propagation/FieldStateCovariance.java @@ -214,7 +214,7 @@ public FieldStateCovariance changeCovarianceType(final FieldOrbit orbit, f final PositionAngleType outAngleType) { // Handle case where the covariance is already expressed in the output type - if (outOrbitType == orbitType && (outAngleType == angleType || outOrbitType == OrbitType.CARTESIAN)) { + if (outOrbitType == orbitType && (outOrbitType == OrbitType.CARTESIAN || outAngleType == angleType)) { if (lof == null) { return new FieldStateCovariance<>(orbitalCovariance, epoch, frame, orbitType, angleType); } @@ -261,6 +261,11 @@ public FieldStateCovariance changeCovarianceFrame(final FieldOrbit orbit, // Verify current covariance frame if (lof != null) { + // Check specific case where output covariance will be the same + if (lofOut == lof) { + return new FieldStateCovariance<>(orbitalCovariance, epoch, lof); + } + // Change the covariance local orbital frame return changeFrameAndCreate(orbit, epoch, lof, lofOut, orbitalCovariance); @@ -297,6 +302,11 @@ public FieldStateCovariance changeCovarianceFrame(final FieldOrbit orbit, } else { + // Check specific case where output covariance will be the same + if (frame == frameOut) { + return new FieldStateCovariance<>(orbitalCovariance, epoch, frame, orbitType, angleType); + } + // Change covariance frame return changeFrameAndCreate(orbit, epoch, frame, frameOut, orbitalCovariance, orbitType, angleType); @@ -435,6 +445,12 @@ private FieldStateCovariance changeTypeAndCreate(final FieldOrbit orbit, final PositionAngleType outAngleType, final FieldMatrix inputCov) { + // Check if type change is really necessary, if not then return input covariance + if (StateCovariance.inputAndOutputOrbitTypesAreCartesian(inOrbitType, outOrbitType) || + StateCovariance.inputAndOutputAreIdentical(inOrbitType, inAngleType, outOrbitType, outAngleType)) { + return new FieldStateCovariance<>(inputCov, date, covFrame, inOrbitType, inAngleType); + } + // Notations: // I: Input orbit type // O: Output orbit type diff --git a/src/main/java/org/orekit/propagation/StateCovariance.java b/src/main/java/org/orekit/propagation/StateCovariance.java index 6b24e7693d..af955dc7d4 100644 --- a/src/main/java/org/orekit/propagation/StateCovariance.java +++ b/src/main/java/org/orekit/propagation/StateCovariance.java @@ -143,6 +143,31 @@ public static void checkFrameAndOrbitTypeConsistency(final Frame covarianceFrame } } + /** + * Checks if input/output orbit and angle types are identical. + * + * @param inOrbitType input orbit type + * @param inAngleType input angle type + * @param outOrbitType output orbit type + * @param outAngleType output angle type + * @return flag defining if input/output orbit and angle types are identical + */ + public static boolean inputAndOutputAreIdentical(final OrbitType inOrbitType, final PositionAngleType inAngleType, + final OrbitType outOrbitType, final PositionAngleType outAngleType) { + return inOrbitType == outOrbitType && inAngleType == outAngleType; + } + + /** + * Checks if input and output orbit types are both {@code OrbitType.CARTESIAN}. + * + * @param inOrbitType input orbit type + * @param outOrbitType output orbit type + * @return flag defining if input and output orbit types are both {@code OrbitType.CARTESIAN} + */ + public static boolean inputAndOutputOrbitTypesAreCartesian(final OrbitType inOrbitType, final OrbitType outOrbitType) { + return inOrbitType == OrbitType.CARTESIAN && outOrbitType == OrbitType.CARTESIAN; + } + /** {@inheritDoc}. */ @Override public AbsoluteDate getDate() { @@ -215,7 +240,7 @@ public StateCovariance changeCovarianceType(final Orbit orbit, final OrbitType o final PositionAngleType outAngleType) { // Handle case where the covariance is already expressed in the output type - if (outOrbitType == orbitType && (outAngleType == angleType || outOrbitType == OrbitType.CARTESIAN)) { + if (outOrbitType == orbitType && (outOrbitType == OrbitType.CARTESIAN || outAngleType == angleType)) { if (lof == null) { return new StateCovariance(orbitalCovariance, epoch, frame, orbitType, angleType); } @@ -261,6 +286,11 @@ public StateCovariance changeCovarianceFrame(final Orbit orbit, final LOF lofOut // Verify current covariance frame if (lof != null) { + // Check specific case where output covariance will be the same + if (lofOut == lof) { + return new StateCovariance(orbitalCovariance, epoch, lof); + } + // Change the covariance local orbital frame return changeFrameAndCreate(orbit, epoch, lof, lofOut, orbitalCovariance); @@ -296,6 +326,11 @@ public StateCovariance changeCovarianceFrame(final Orbit orbit, final Frame fram } else { + // Check specific case where output covariance will be the same + if (frame == frameOut) { + return new StateCovariance(orbitalCovariance, epoch, frame, orbitType, angleType); + } + // Change covariance frame return changeFrameAndCreate(orbit, epoch, frame, frameOut, orbitalCovariance, orbitType, angleType); @@ -393,12 +428,19 @@ public StateCovariance shiftedBy(final Orbit orbit, final double dt) { * @param inputCov input covariance * @return the covariance expressed in the target orbit type with the target position angle */ - private static StateCovariance changeTypeAndCreate(final Orbit orbit, final AbsoluteDate date, + private static StateCovariance changeTypeAndCreate(final Orbit orbit, + final AbsoluteDate date, final Frame covFrame, final OrbitType inOrbitType, final PositionAngleType inAngleType, final OrbitType outOrbitType, final PositionAngleType outAngleType, final RealMatrix inputCov) { + // Check if type change is really necessary, if not then return input covariance + if (inputAndOutputOrbitTypesAreCartesian(inOrbitType, outOrbitType) || + inputAndOutputAreIdentical(inOrbitType, inAngleType, outOrbitType, outAngleType)) { + return new StateCovariance(inputCov, date, covFrame, inOrbitType, inAngleType); + } + // Notations: // I: Input orbit type // O: Output orbit type diff --git a/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java b/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java index d8b176c9f6..2845bc56ae 100644 --- a/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java +++ b/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java @@ -445,6 +445,20 @@ public SpacecraftState propagate(final AbsoluteDate tStart, final AbsoluteDate t } + /** Reset initial state with a given propagation type. + * + *

    By default this method returns the same as {@link #resetInitialState(SpacecraftState)} + *

    Its purpose is mostly to be derived in DSSTPropagator + * + * @param state new initial state to consider + * @param stateType type of the new state (mean or osculating) + * @since 12.1.3 + */ + public void resetInitialState(final SpacecraftState state, final PropagationType stateType) { + // Default behavior, do not take propagation type into account + resetInitialState(state); + } + /** Set up State Transition Matrix and Jacobian matrix handling. * @since 11.1 */ @@ -505,7 +519,7 @@ private SpacecraftState integrateDynamics(final AbsoluteDate tEnd, final boolean finalState = updateAdditionalStatesAndDerivatives(finalState, mathFinalState); if (resetAtEnd || forceResetAtEnd) { - resetInitialState(finalState); + resetInitialState(finalState, propagationType); setStartDate(finalState.getDate()); } diff --git a/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java b/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java index b5ca913ae1..b5864364cc 100644 --- a/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java +++ b/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java @@ -450,6 +450,20 @@ public FieldSpacecraftState propagate(final FieldAbsoluteDate tStart, fina } + /** Reset initial state with a given propagation type. + * + *

    By default this method returns the same as method resetInitialState(FieldSpacecraftState) + *

    Its purpose is mostly to be derived in FieldDSSTPropagator + * + * @param state new initial state to consider + * @param stateType type of the new state (mean or osculating) + * @since 12.1.3 + */ + public void resetInitialState(final FieldSpacecraftState state, final PropagationType stateType) { + // Default behavior, do not take propagation type into account + resetInitialState(state); + } + /** Propagation with or without event detection. * @param tEnd target date to which orbit should be propagated * @return state at end of propagation @@ -500,7 +514,7 @@ private FieldSpacecraftState integrateDynamics(final FieldAbsoluteDate tEn finalState = updateAdditionalStatesAndDerivatives(finalState, mathFinalState); if (resetAtEnd) { - resetInitialState(finalState); + resetInitialState(finalState, propagationType); setStartDate(finalState.getDate()); } diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTHarvester.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTHarvester.java index 0c082d4ac6..fe4f24fa2c 100644 --- a/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTHarvester.java +++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTHarvester.java @@ -274,8 +274,11 @@ public void initializeFieldShortPeriodTerms(final SpacecraftState reference, type, dsParameters); // create a copy of the list to protect against inadvertent modification - fieldShortPeriodTerms.computeIfAbsent(forceModel, x -> new ArrayList<>()) - .addAll(terms); + final List> list; + synchronized (fieldShortPeriodTerms) { + list = fieldShortPeriodTerms.computeIfAbsent(forceModel, x -> new ArrayList<>()); + } + list.addAll(terms); } @@ -324,8 +327,10 @@ public void setReferenceState(final SpacecraftState reference) { final Gradient zero = dsState.getDate().getField().getZero(); final Gradient[] shortPeriod = new Gradient[6]; Arrays.fill(shortPeriod, zero); - final List> terms = fieldShortPeriodTerms - .computeIfAbsent(forceModel, x -> new ArrayList<>(0)); + final List> terms; + synchronized (fieldShortPeriodTerms) { + terms = fieldShortPeriodTerms.computeIfAbsent(forceModel, x -> new ArrayList<>(0)); + } for (final FieldShortPeriodTerms spt : terms) { final Gradient[] spVariation = spt.value(dsState.getOrbit()); for (int i = 0; i < spVariation .length; i++) { diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagator.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagator.java index a031e60225..e769d8104b 100644 --- a/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagator.java +++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagator.java @@ -16,6 +16,14 @@ */ package org.orekit.propagation.semianalytical.dsst; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collection; +import java.util.Collections; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + import org.hipparchus.linear.RealMatrix; import org.hipparchus.ode.ODEIntegrator; import org.hipparchus.ode.ODEStateAndDerivative; @@ -28,7 +36,6 @@ import org.orekit.attitudes.AttitudeProvider; import org.orekit.data.DataContext; import org.orekit.errors.OrekitException; -import org.orekit.errors.OrekitInternalError; import org.orekit.errors.OrekitMessages; import org.orekit.frames.Frame; import org.orekit.orbits.EquinoctialOrbit; @@ -59,14 +66,6 @@ import org.orekit.utils.TimeSpanMap; import org.orekit.utils.TimeSpanMap.Span; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Collection; -import java.util.Collections; -import java.util.HashSet; -import java.util.List; -import java.util.Set; - /** * This class propagates {@link org.orekit.orbits.Orbit orbits} using the DSST theory. *

    @@ -278,17 +277,7 @@ public void setInitialState(final SpacecraftState initialState) { */ public void setInitialState(final SpacecraftState initialState, final PropagationType stateType) { - switch (stateType) { - case MEAN: - initialIsOsculating = false; - break; - case OSCULATING: - initialIsOsculating = true; - break; - default: - throw new OrekitInternalError(null); - } - resetInitialState(initialState); + resetInitialState(initialState, stateType); } /** Reset the initial state. @@ -305,6 +294,20 @@ public void resetInitialState(final SpacecraftState state) { super.setStartDate(state.getDate()); } + /** {@inheritDoc}. + * + *

    Change parameter {@link #initialIsOsculating()} accordingly + * @since 12.1.3 + */ + @Override + public void resetInitialState(final SpacecraftState state, final PropagationType stateType) { + // Reset initial state + resetInitialState(state); + + // Change state of initial osculating, if needed + initialIsOsculating = stateType == PropagationType.OSCULATING; + } + /** Set the selected short periodic coefficients that must be stored as additional states. * @param selectedCoefficients short periodic coefficients that must be stored as additional states * (null means no coefficients are selected, empty set means all coefficients are selected) diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagator.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagator.java index 30eac9acc0..ad51495d81 100644 --- a/src/main/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagator.java +++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagator.java @@ -16,6 +16,14 @@ */ package org.orekit.propagation.semianalytical.dsst; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collection; +import java.util.Collections; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + import org.hipparchus.CalculusFieldElement; import org.hipparchus.Field; import org.hipparchus.ode.FieldODEIntegrator; @@ -58,14 +66,6 @@ import org.orekit.utils.ParameterObserver; import org.orekit.utils.TimeSpanMap; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Collection; -import java.util.Collections; -import java.util.HashSet; -import java.util.List; -import java.util.Set; - /** * This class propagates {@link org.orekit.orbits.FieldOrbit orbits} using the DSST theory. *

    @@ -309,17 +309,7 @@ public void setInitialState(final FieldSpacecraftState initialState) { */ public void setInitialState(final FieldSpacecraftState initialState, final PropagationType stateType) { - switch (stateType) { - case MEAN: - initialIsOsculating = false; - break; - case OSCULATING: - initialIsOsculating = true; - break; - default: - throw new OrekitInternalError(null); - } - resetInitialState(initialState); + resetInitialState(initialState, stateType); } /** Reset the initial state. @@ -336,6 +326,20 @@ public void resetInitialState(final FieldSpacecraftState state) { super.setStartDate(state.getDate()); } + /** {@inheritDoc}. + * + *

    Change parameter {@link #initialIsOsculating()} accordingly + * @since 12.1.3 + */ + @Override + public void resetInitialState(final FieldSpacecraftState state, final PropagationType stateType) { + // Reset initial state + resetInitialState(state); + + // Change state of initial osculating, if needed + initialIsOsculating = stateType == PropagationType.OSCULATING; + } + /** Set the selected short periodic coefficients that must be stored as additional states. * @param selectedCoefficients short periodic coefficients that must be stored as additional states * (null means no coefficients are selected, empty set means all coefficients are selected) diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/FieldHansenTesseralLinear.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/FieldHansenTesseralLinear.java index b1365fd255..1995bdf2df 100644 --- a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/FieldHansenTesseralLinear.java +++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/FieldHansenTesseralLinear.java @@ -303,7 +303,7 @@ private FieldGradient getValueGradient(final T e2, final T chi, final T chi2) final T coef = zero.subtract(mnm1 + 1.5); final T derivative = coef.multiply(chi2).multiply(value). add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(0)).divide(chi)); - return new FieldGradient<>(value, derivative); + return new FieldGradient(value, derivative); } /** Generate the serie expansion in e². diff --git a/src/main/java/org/orekit/time/AbstractTimeScales.java b/src/main/java/org/orekit/time/AbstractTimeScales.java index b80a70bddd..7518fcaec6 100644 --- a/src/main/java/org/orekit/time/AbstractTimeScales.java +++ b/src/main/java/org/orekit/time/AbstractTimeScales.java @@ -78,16 +78,18 @@ protected abstract EOPHistory getEopHistory(IERSConventions conventions, @Override public UT1Scale getUT1(final IERSConventions conventions, final boolean simpleEOP) { - return ut1Map.computeIfAbsent( - new Pair<>(conventions, simpleEOP), - k -> getUT1(getEopHistory(conventions, simpleEOP))); + final Pair key = new Pair<>(conventions, simpleEOP); + synchronized (this) { + return ut1Map.computeIfAbsent(key, k -> getUT1(getEopHistory(conventions, simpleEOP))); + } } @Override public GMSTScale getGMST(final IERSConventions conventions, final boolean simpleEOP) { - return gmstMap.computeIfAbsent( - new Pair<>(conventions, simpleEOP), - k -> new GMSTScale(getUT1(conventions, simpleEOP))); + final Pair key = new Pair<>(conventions, simpleEOP); + synchronized (this) { + return gmstMap.computeIfAbsent(key, k -> new GMSTScale(getUT1(conventions, simpleEOP))); + } } @Override diff --git a/src/main/java/org/orekit/time/DateComponents.java b/src/main/java/org/orekit/time/DateComponents.java index f836e6e2bd..d0b648a415 100644 --- a/src/main/java/org/orekit/time/DateComponents.java +++ b/src/main/java/org/orekit/time/DateComponents.java @@ -17,7 +17,6 @@ package org.orekit.time; import java.io.Serializable; -import java.text.DecimalFormat; import java.text.DecimalFormatSymbols; import java.util.Locale; import java.util.regex.Matcher; @@ -124,12 +123,6 @@ public class DateComponents implements Serializable, Comparable /** Formatting symbols used in {@link #toString()}. */ private static final DecimalFormatSymbols US_SYMBOLS = new DecimalFormatSymbols(Locale.US); - /** Format for years. */ - private static final DecimalFormat FOUR_DIGITS = new DecimalFormat("0000", US_SYMBOLS); - - /** Format for months and days. */ - private static final DecimalFormat TWO_DIGITS = new DecimalFormat("00", US_SYMBOLS); - /** Offset between J2000 epoch and modified julian day epoch. */ private static final int MJD_TO_J2000 = 51544; @@ -473,11 +466,7 @@ public int getDayOfYear() { * @return string representation of the date. */ public String toString() { - return new StringBuilder(). - append(FOUR_DIGITS.format(year)).append('-'). - append(TWO_DIGITS.format(month)).append('-'). - append(TWO_DIGITS.format(day)). - toString(); + return String.format(Locale.US, "%04d-%02d-%02d", year, month, day); } /** {@inheritDoc} */ diff --git a/src/main/java/org/orekit/time/PreloadedTimeScales.java b/src/main/java/org/orekit/time/PreloadedTimeScales.java index 48631ef9ed..ca32861d0b 100644 --- a/src/main/java/org/orekit/time/PreloadedTimeScales.java +++ b/src/main/java/org/orekit/time/PreloadedTimeScales.java @@ -110,10 +110,12 @@ public UTCScale getUTC() { @Override protected EOPHistory getEopHistory(final IERSConventions conventions, final boolean simpleEOP) { + final Collection data; + synchronized (this) { + data = eopMap.computeIfAbsent(conventions, c -> eopSupplier.apply(c, this)); + } return new EOPHistory(conventions, EOPHistory.DEFAULT_INTERPOLATION_DEGREE, - eopMap.computeIfAbsent(conventions, c -> eopSupplier.apply(c, this)), - simpleEOP, - this); + data, simpleEOP, this); } @Override diff --git a/src/main/java/org/orekit/utils/units/UnitsCache.java b/src/main/java/org/orekit/utils/units/UnitsCache.java index 1972a14872..5e61d68dcd 100644 --- a/src/main/java/org/orekit/utils/units/UnitsCache.java +++ b/src/main/java/org/orekit/utils/units/UnitsCache.java @@ -50,7 +50,9 @@ public Unit getUnits(final String specification) { return Unit.NONE; } - return cache.computeIfAbsent(specification, s -> Unit.parse(specification)); + synchronized (cache) { + return cache.computeIfAbsent(specification, s -> Unit.parse(specification)); + } } diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_ca.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_ca.utf8 index 56517df3ec..355f9fd3ed 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_ca.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_ca.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_da.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_da.utf8 index 03dfe749fe..1f2a48860c 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_da.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_da.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_de.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_de.utf8 index 3732af6ee1..3f92c19873 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_de.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_de.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_el.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_el.utf8 index 5c3a6ccb88..e1ae7eab9e 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_el.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_el.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_en.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_en.utf8 index d4a40f432c..08cb7c8b5a 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_en.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_en.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = number of planes {0} is inconsistent with number of # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = Infinite value appears during computation of atmospheric density in NRLMSISE00 model + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = field "{0}" is too long, maximum length is {1} characters diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_es.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_es.utf8 index 0c6363d924..3021e0fa79 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_es.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_es.utf8 @@ -888,3 +888,5 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_fr.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_fr.utf8 index 5c8ccff1ef..34d7d7bce8 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_fr.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_fr.utf8 @@ -888,3 +888,5 @@ WALKER_INCONSISTENT_PLANES = le nombre {0} de plans est incohérent avec le nomb # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = Valeur infinie lors du calcul de la densité atmosphérique dans le modèle NRLMSISE00 +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = le champ « {0} » est trop long, la longueur maximale est de {1} caractères diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_gl.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_gl.utf8 index 01cef4f554..10996b5300 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_gl.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_gl.utf8 @@ -888,3 +888,5 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_it.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_it.utf8 index 839812ba7a..8a71ca5bf5 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_it.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_it.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_no.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_no.utf8 index b5836a1f62..996b7416c9 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_no.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_no.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/main/resources/assets/org/orekit/localization/OrekitMessages_ro.utf8 b/src/main/resources/assets/org/orekit/localization/OrekitMessages_ro.utf8 index aaf4d32ab9..bb6be4cd49 100644 --- a/src/main/resources/assets/org/orekit/localization/OrekitMessages_ro.utf8 +++ b/src/main/resources/assets/org/orekit/localization/OrekitMessages_ro.utf8 @@ -887,3 +887,6 @@ WALKER_INCONSISTENT_PLANES = numărul de planuri {0} este neconsecvent cu număr # Infinite value appears during computation of atmospheric density in NRLMSISE00 model INFINITE_NRLMSISE00_DENSITY = valori infinite au fost detectate în timpul calculului densității atmosferice pentru modelul NRLMSISE00 + +# field "{0}" is too long, maximum length is {1} characters +FIELD_TOO_LONG = diff --git a/src/site/markdown/downloads.md.vm b/src/site/markdown/downloads.md.vm index 80738794d1..197349a5b0 100644 --- a/src/site/markdown/downloads.md.vm +++ b/src/site/markdown/downloads.md.vm @@ -45,7 +45,7 @@ with groupID org.orekit and artifactId orekit so maven internal mechanism will download automatically all artifacts and dependencies as required. -#set ( $versions = {"12.1.2": "2024-07-13", "12.1.1": "2024-06-25", "12.1": "2024-06-24", "12.0.2": "2024-03-15", "12.0.1": "2023-12-31", "12.0": "2023-11-08", "11.3.3": "2023-06-30", "11.3.2": "2023-02-17", "11.3.1": "2022-12-24", "11.3": "2022-10-25", "11.2.1": "2022-08-01", "11.2": "2022-06-20", "11.1.2": "2022-04-27", "11.1.1": "2022-03-17", "11.1": "2022-02-14", "11.0.2": "2021-11-24", "11.0.1": "2021-10-22", "11.0": "2021-09-20", "10.3.1": "2021-06-16", "10.3": "2020-12-21", "10.2": "2020-07-14", "10.1": "2020-02-19", "10.0": "2019-06-24", "9.3.1": "2019-03-16", "9.3": "2019-01-25", "9.2": "2018-05-26","9.1": "2017-11-26","9.0.1": "2017-11-03","9.0": "2017-07-26","8.0.1": "2017-11-03","8.0": "2016-06-30","7.2.1": "2017-11-03","7.2": "2016-04-05","7.1": "2016-02-07","7.0": "2015-01-11","6.1": "2013-12-13","6.0": "2013-04-23","5.0.3": "2011-07-13","5.0.2": "2011-07-11","5.0.1": "2011-04-18"} ) +#set ( $versions = {"12.1.3": "2024-09-04", "12.1.2": "2024-07-13", "12.1.1": "2024-06-25", "12.1": "2024-06-24", "12.0.2": "2024-03-15", "12.0.1": "2023-12-31", "12.0": "2023-11-08", "11.3.3": "2023-06-30", "11.3.2": "2023-02-17", "11.3.1": "2022-12-24", "11.3": "2022-10-25", "11.2.1": "2022-08-01", "11.2": "2022-06-20", "11.1.2": "2022-04-27", "11.1.1": "2022-03-17", "11.1": "2022-02-14", "11.0.2": "2021-11-24", "11.0.1": "2021-10-22", "11.0": "2021-09-20", "10.3.1": "2021-06-16", "10.3": "2020-12-21", "10.2": "2020-07-14", "10.1": "2020-02-19", "10.0": "2019-06-24", "9.3.1": "2019-03-16", "9.3": "2019-01-25", "9.2": "2018-05-26","9.1": "2017-11-26","9.0.1": "2017-11-03","9.0": "2017-07-26","8.0.1": "2017-11-03","8.0": "2016-06-30","7.2.1": "2017-11-03","7.2": "2016-04-05","7.1": "2016-02-07","7.0": "2015-01-11","6.1": "2013-12-13","6.0": "2013-04-23","5.0.3": "2011-07-13","5.0.2": "2011-07-11","5.0.1": "2011-04-18"} ) #foreach( $version in $versions.entrySet() ) | package | link | diff --git a/src/site/markdown/faq.md b/src/site/markdown/faq.md index e5993cd9bf..c485983e63 100644 --- a/src/site/markdown/faq.md +++ b/src/site/markdown/faq.md @@ -158,6 +158,7 @@ Math to Hipparchus Orekit 12.1 | Hipparchus 3.1 Orekit 12.1.1 | Hipparchus 3.1 Orekit 12.1.2 | Hipparchus 3.1 + Orekit 12.1.3 | Hipparchus 3.1 ### Maven failed to compile Orekit and complained about a missing artifact. diff --git a/src/test/java/org/orekit/errors/OrekitMessagesTest.java b/src/test/java/org/orekit/errors/OrekitMessagesTest.java index ffa5eac8b5..45a4ccebe1 100644 --- a/src/test/java/org/orekit/errors/OrekitMessagesTest.java +++ b/src/test/java/org/orekit/errors/OrekitMessagesTest.java @@ -32,7 +32,7 @@ public class OrekitMessagesTest { @Test public void testMessageNumber() { - Assertions.assertEquals(296, OrekitMessages.values().length); + Assertions.assertEquals(297, OrekitMessages.values().length); } @Test diff --git a/src/test/java/org/orekit/estimation/measurements/EstimatedMeasurementTest.java b/src/test/java/org/orekit/estimation/measurements/EstimatedMeasurementTest.java new file mode 100644 index 0000000000..ad9e31546e --- /dev/null +++ b/src/test/java/org/orekit/estimation/measurements/EstimatedMeasurementTest.java @@ -0,0 +1,36 @@ +package org.orekit.estimation.measurements; + +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.EnumSource; +import org.mockito.Mockito; +import org.orekit.propagation.SpacecraftState; +import org.orekit.utils.TimeStampedPVCoordinates; + +class EstimatedMeasurementTest { + + @ParameterizedTest + @EnumSource(EstimatedMeasurementBase.Status.class) + void testConstructor(final EstimatedMeasurementBase.Status expectedtStatus) { + // GIVEN + final Position mockedPosition = Mockito.mock(Position.class); + final int expectedIteration = 2; + final int expectedCount = 3; + final SpacecraftState[] states = new SpacecraftState[0]; + final double[] expectedEstimatedValue = new double[] {1., 2., 3.}; + final TimeStampedPVCoordinates[] timeStampedPVCoordinates = new TimeStampedPVCoordinates[0]; + final EstimatedMeasurementBase estimatedMeasurementBase = new EstimatedMeasurementBase<>(mockedPosition, + expectedIteration, expectedCount, states, timeStampedPVCoordinates); + estimatedMeasurementBase.setEstimatedValue(expectedEstimatedValue); + estimatedMeasurementBase.setStatus(expectedtStatus); + // WHEN + final EstimatedMeasurement estimatedMeasurement = new EstimatedMeasurement<>(estimatedMeasurementBase); + // THEN + Assertions.assertEquals(mockedPosition, estimatedMeasurement.getObservedMeasurement()); + Assertions.assertEquals(expectedCount, estimatedMeasurement.getCount()); + Assertions.assertEquals(expectedIteration, estimatedMeasurement.getIteration()); + Assertions.assertEquals(expectedtStatus, estimatedMeasurement.getStatus()); + Assertions.assertArrayEquals(expectedEstimatedValue, estimatedMeasurement.getEstimatedValue()); + } + +} diff --git a/src/test/java/org/orekit/estimation/sequential/SequentialNumericalOrbitDeterminationTest.java b/src/test/java/org/orekit/estimation/sequential/SequentialNumericalOrbitDeterminationTest.java index 346c50d281..6fcdfcf8b4 100644 --- a/src/test/java/org/orekit/estimation/sequential/SequentialNumericalOrbitDeterminationTest.java +++ b/src/test/java/org/orekit/estimation/sequential/SequentialNumericalOrbitDeterminationTest.java @@ -422,76 +422,6 @@ public void testW3BExtended() throws URISyntaxException, IOException { predictionDistanceAccuracy, predictionVelocityAccuracy, false, false); } - @Test - public void testW3BUnscented() throws URISyntaxException, IOException { - // Batch LS result: -0.2154; - final double dragCoef = -0.0214; - - // Batch LS results: 8.002e-6 - final double leakAccelerationNorm0 = 5.954e-6; - - // Batch LS results: 3.058e-10 - final double leakAccelerationNorm1 = 1.619e-10; - - // Batch LS results - // final double[] CastleAzElBias = { 0.062701342, -0.003613508 }; - // final double CastleRangeBias = 11274.4677; - final double[] castleAzElBias = { 0.062344, -0.004106}; - final double castleRangeBias = 11333.0998; - - // Batch LS results - // final double[] FucAzElBias = { -0.053526137, 0.075483886 }; - // final double FucRangeBias = 13467.8256; - final double[] fucAzElBias = { -0.053870, 0.075641 }; - final double fucRangeBias = 13461.7291; - - // Batch LS results - // final double[] KumAzElBias = { -0.023574208, -0.054520756 }; - // final double KumRangeBias = 13512.57594; - final double[] kumAzElBias = { -0.023393, -0.055078 }; - final double kumRangeBias = 13515.6244; - - // Batch LS results - // final double[] PreAzElBias = { 0.030201539, 0.009747877 }; - // final double PreRangeBias = 13594.11889; - final double[] preAzElBias = { 0.030250, 0.010083 }; - final double preRangeBias = 13534.0334; - - // Batch LS results - // final double[] UraAzElBias = { 0.167814449, -0.12305252 }; - // final double UraRangeBias = 13450.26738; - final double[] uraAzElBias = { 0.167700, -0.122408 }; - final double uraRangeBias = 13417.6695; - - //statistics for the range residual (min, max, mean, std) - final double[] refStatRange = { -144.9733, 14.7486, -3.8990, 11.9050 }; - - //statistics for the azimuth residual (min, max, mean, std) - final double[] refStatAzi = { -0.041873, 0.018087, -0.004536, 0.008995 }; - - //statistics for the elevation residual (min, max, mean, std) - final double[] refStatEle = { -0.025583, 0.043560, 0.001857, 0.010625 }; - - // Expected covariance - final double dragVariance = 0.018813; - final double leakXVariance = 2.117E-13; - final double leakYVariance = 5.540E-13; - final double leakZVariance = 1.73244E-11; - - // Prediction position/velocity accuracies - // FIXME: debug - Comparison with batch LS is bad - final double predictionDistanceAccuracy = 285.31; - final double predictionVelocityAccuracy = 0.101; - - testW3B(dragCoef, leakAccelerationNorm0, leakAccelerationNorm1, - castleAzElBias, castleRangeBias, fucAzElBias, fucRangeBias, kumAzElBias, kumRangeBias, - preAzElBias, preRangeBias, uraAzElBias, uraRangeBias, - refStatRange, refStatAzi, refStatEle, dragVariance, - leakXVariance, leakYVariance, leakZVariance, - predictionDistanceAccuracy, predictionVelocityAccuracy, false, true); - } - - // Orbit determination for range, azimuth elevation measurements protected void testW3B(final double dragCoef, final double leakAccelerationNorm0, final double leakAccelerationNorm1, final double[] castleAzElBias, final double castleRangeBias, diff --git a/src/test/java/org/orekit/files/rinex/navigation/NavigationFileParserTest.java b/src/test/java/org/orekit/files/rinex/navigation/NavigationFileParserTest.java index 20827135aa..bbf2c55a2b 100644 --- a/src/test/java/org/orekit/files/rinex/navigation/NavigationFileParserTest.java +++ b/src/test/java/org/orekit/files/rinex/navigation/NavigationFileParserTest.java @@ -21,7 +21,6 @@ import java.lang.reflect.Field; import java.lang.reflect.InvocationTargetException; import java.lang.reflect.Method; -import java.net.URISyntaxException; import java.util.Arrays; import java.util.List; @@ -92,7 +91,7 @@ public void testWorkAroundWrongFormatNumber() throws IOException { } @Test - public void testGpsRinex301Truncated() throws URISyntaxException, IOException { + public void testGpsRinex301Truncated() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_GPS_Rinex301.n"; @@ -107,7 +106,7 @@ public void testGpsRinex301Truncated() throws URISyntaxException, IOException { } @Test - public void testGpsRinex301() throws URISyntaxException, IOException { + public void testGpsRinex301() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_GPS_Rinex301.n"; @@ -195,7 +194,7 @@ public void testGpsRinex301() throws URISyntaxException, IOException { } @Test - public void testGpsRinex400() throws URISyntaxException, IOException { + public void testGpsRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_GPS_Rinex400.n"; @@ -253,7 +252,7 @@ public void testGpsRinex400() throws URISyntaxException, IOException { } @Test - public void testSBASRinex301() throws URISyntaxException, IOException { + public void testSBASRinex301() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_SBAS_Rinex301.n"; @@ -320,7 +319,7 @@ public void testSBASRinex301() throws URISyntaxException, IOException { } @Test - public void testBeidouRinex302() throws URISyntaxException, IOException { + public void testBeidouRinex302() throws IOException { final String ex = "/gnss/navigation/Example_Beidou_Rinex302.n"; final RinexNavigation file = new RinexNavigationParser(). @@ -390,7 +389,7 @@ public void testBeidouRinex302() throws URISyntaxException, IOException { } @Test - public void testBeidouRinex400() throws URISyntaxException, IOException { + public void testBeidouRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Beidou_Rinex400.n"; @@ -479,7 +478,7 @@ public void testBeidouRinex400() throws URISyntaxException, IOException { } @Test - public void testGalileoRinex302() throws URISyntaxException, IOException { + public void testGalileoRinex302() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Galileo_Rinex302.n"; @@ -565,7 +564,7 @@ public void testGalileoRinex302() throws URISyntaxException, IOException { } @Test - public void testGalileoRinex400() throws URISyntaxException, IOException { + public void testGalileoRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Galileo_Rinex400.n"; @@ -603,7 +602,7 @@ public void testGalileoRinex400() throws URISyntaxException, IOException { } @Test - public void testQZSSRinex302() throws URISyntaxException, IOException { + public void testQZSSRinex302() throws IOException { final String ex = "/gnss/navigation/Example_QZSS_Rinex302.n"; final RinexNavigation file = new RinexNavigationParser(). @@ -692,7 +691,7 @@ public void testQZSSRinex302() throws URISyntaxException, IOException { } @Test - public void testQZSSRinex400() throws URISyntaxException, IOException { + public void testQZSSRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_QZSS_Rinex400.n"; @@ -745,7 +744,7 @@ public void testQZSSRinex400() throws URISyntaxException, IOException { } @Test - public void testGLONASSRinex303() throws URISyntaxException, IOException { + public void testGLONASSRinex303() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Glonass_Rinex303.n"; @@ -806,7 +805,7 @@ public void testGLONASSRinex303() throws URISyntaxException, IOException { } @Test - public void testIRNSSRinex303() throws URISyntaxException, IOException { + public void testIRNSSRinex303() throws IOException { final String ex = "/gnss/navigation/Example_IRNSS_Rinex303.n"; final RinexNavigation file = new RinexNavigationParser(). @@ -876,7 +875,7 @@ public void testIRNSSRinex303() throws URISyntaxException, IOException { } @Test - public void testIRNSSRinex400() throws URISyntaxException, IOException { + public void testIRNSSRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_IRNSS_Rinex400.n"; @@ -912,7 +911,7 @@ public void testIRNSSRinex400() throws URISyntaxException, IOException { } @Test - public void testMixedRinex304() throws URISyntaxException, IOException { + public void testMixedRinex304() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Mixed_Rinex304.n"; @@ -973,7 +972,7 @@ public void testMixedRinex304() throws URISyntaxException, IOException { } @Test - public void testMixedRinex305() throws URISyntaxException, IOException { + public void testMixedRinex305() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Mixed_Rinex305.n"; @@ -1033,7 +1032,7 @@ public void testMixedRinex305() throws URISyntaxException, IOException { } @Test - public void testQZSSRinex304() throws URISyntaxException, IOException { + public void testQZSSRinex304() throws IOException { final String ex = "/gnss/navigation/Example_QZSS_Rinex304.n"; final RinexNavigation file = new RinexNavigationParser(). @@ -1104,7 +1103,7 @@ public void testQZSSRinex304() throws URISyntaxException, IOException { } @Test - public void testGpsRinex304() throws URISyntaxException, IOException { + public void testGpsRinex304() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_GPS_Rinex304.n"; @@ -1194,7 +1193,7 @@ public void testGpsRinex304() throws URISyntaxException, IOException { } @Test - public void testGalileoRinex304() throws URISyntaxException, IOException { + public void testGalileoRinex304() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Galileo_Rinex304.n"; @@ -1273,7 +1272,7 @@ public void testGalileoRinex304() throws URISyntaxException, IOException { } @Test - public void testSBASRinex304() throws URISyntaxException, IOException { + public void testSBASRinex304() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_SBAS_Rinex304.n"; @@ -1340,7 +1339,7 @@ public void testSBASRinex304() throws URISyntaxException, IOException { } @Test - public void testSBASRinex400() throws URISyntaxException, IOException { + public void testSBASRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_SBAS_Rinex400.n"; @@ -1377,7 +1376,7 @@ public void testSBASRinex400() throws URISyntaxException, IOException { } @Test - public void testIRNSSRinex304() throws URISyntaxException, IOException { + public void testIRNSSRinex304() throws IOException { final String ex = "/gnss/navigation/Example_IRNSS_Rinex304.n"; final RinexNavigation file = new RinexNavigationParser(). @@ -1447,7 +1446,7 @@ public void testIRNSSRinex304() throws URISyntaxException, IOException { } @Test - public void testBeidouRinex304() throws URISyntaxException, IOException { + public void testBeidouRinex304() throws IOException { final String ex = "/gnss/navigation/Example_Beidou_Rinex304.n"; final RinexNavigation file = new RinexNavigationParser(). @@ -1517,7 +1516,7 @@ public void testBeidouRinex304() throws URISyntaxException, IOException { } @Test - public void testStoRinex400() throws URISyntaxException, IOException { + public void testStoRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Sto_Rinex400.n"; @@ -1575,7 +1574,13 @@ public void testStoRinex400() throws URISyntaxException, IOException { Assertions.assertEquals(TimeSystem.BEIDOU, list.get( 2).getDefinedTimeSystem()); Assertions.assertEquals(TimeSystem.GPS, list.get( 2).getReferenceTimeSystem()); Assertions.assertEquals(TimeSystem.BEIDOU, list.get( 3).getDefinedTimeSystem()); + Assertions.assertEquals("BDT", list.get( 3).getDefinedTimeSystem().getKey()); + Assertions.assertEquals("BD", list.get( 3).getDefinedTimeSystem().getTwoLettersCode()); + Assertions.assertEquals("C", list.get( 3).getDefinedTimeSystem().getOneLetterCode()); Assertions.assertEquals(TimeSystem.UTC, list.get( 3).getReferenceTimeSystem()); + Assertions.assertEquals("UTC", list.get( 3).getReferenceTimeSystem().getKey()); + Assertions.assertEquals("UT", list.get( 3).getReferenceTimeSystem().getTwoLettersCode()); + Assertions.assertNull( list.get( 3).getReferenceTimeSystem().getOneLetterCode()); Assertions.assertEquals(UtcId.NTSC, list.get( 3).getUtcId()); Assertions.assertEquals(TimeSystem.GALILEO, list.get( 4).getDefinedTimeSystem()); Assertions.assertEquals(TimeSystem.GPS, list.get( 4).getReferenceTimeSystem()); @@ -1612,7 +1617,7 @@ public void testStoRinex400() throws URISyntaxException, IOException { } @Test - public void testEopRinex400() throws URISyntaxException, IOException { + public void testEopRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Eop_Rinex400.n"; @@ -1704,7 +1709,7 @@ public void testEopRinex400() throws URISyntaxException, IOException { } @Test - public void testIonRinex400() throws URISyntaxException, IOException { + public void testIonRinex400() throws IOException { // Parse file final String ex = "/gnss/navigation/Example_Ion_Rinex400.n"; @@ -1804,7 +1809,7 @@ public void testIonRinex400() throws URISyntaxException, IOException { } @Test - public void testGPSRinex2() throws URISyntaxException, IOException { + public void testGPSRinex2() throws IOException { // Parse file final String ex = "/gnss/navigation/brdc0130.22n"; @@ -1862,7 +1867,7 @@ public void testGPSRinex2() throws URISyntaxException, IOException { } @Test - public void testGlonassRinex2() throws URISyntaxException, IOException { + public void testGlonassRinex2() throws IOException { // Parse file final String ex = "/gnss/navigation/brdc0130.22g"; @@ -1940,7 +1945,7 @@ public void testUnknownRinexVersion() throws IOException { Assertions.fail("an exception should have been thrown"); } catch (OrekitException oe) { Assertions.assertEquals(OrekitMessages.UNSUPPORTED_FILE_FORMAT_VERSION, oe.getSpecifier()); - Assertions.assertEquals(9.99, ((Double) oe.getParts()[0]).doubleValue(), 1.0e-10); + Assertions.assertEquals(9.99, (Double) oe.getParts()[0], 1.0e-10); Assertions.assertEquals(ex, oe.getParts()[1]); } } @@ -1999,7 +2004,7 @@ public void testDefensiveProgrammingExceptions() { m.invoke(sbasParserField.get(null), "", parseInfo); Assertions.fail("an exception should have been thrown"); } catch (InvocationTargetException e) { - Assertions.assertTrue(e.getCause() instanceof OrekitInternalError); + Assertions.assertInstanceOf(OrekitInternalError.class, e.getCause()); } } diff --git a/src/test/java/org/orekit/files/rinex/observation/RinexObservationWriterTest.java b/src/test/java/org/orekit/files/rinex/observation/RinexObservationWriterTest.java index 63a2de1fc5..6a66695e85 100644 --- a/src/test/java/org/orekit/files/rinex/observation/RinexObservationWriterTest.java +++ b/src/test/java/org/orekit/files/rinex/observation/RinexObservationWriterTest.java @@ -65,6 +65,23 @@ public void testWriteHeaderTwice() throws IOException { } } + @Test + public void testTooLongAgency() throws IOException { + final RinexObservation robs = load("rinex/bbbb0000.00o"); + robs.getHeader().setAgencyName("much too long agency name exceeding 40 characters"); + final CharArrayWriter caw = new CharArrayWriter(); + try (RinexObservationWriter writer = new RinexObservationWriter(caw, "dummy")) { + writer.setReceiverClockModel(robs.extractClockModel(2)); + writer.prepareComments(robs.getComments()); + writer.writeHeader(robs.getHeader()); + Assertions.fail("an exception should have been thrown"); + } catch (OrekitException oe) { + Assertions.assertEquals(OrekitMessages.FIELD_TOO_LONG, oe.getSpecifier()); + Assertions.assertEquals("much too long agency name exceeding 40 characters", oe.getParts()[0]); + Assertions.assertEquals(40, (Integer) oe.getParts()[1]); + } + } + @Test public void testNoWriteHeader() throws IOException { final RinexObservation robs = load("rinex/aiub0000.00o"); diff --git a/src/test/java/org/orekit/models/earth/ionosphere/NeQuickModelTest.java b/src/test/java/org/orekit/models/earth/ionosphere/NeQuickModelTest.java index c4e45e7b8b..d242d42f98 100644 --- a/src/test/java/org/orekit/models/earth/ionosphere/NeQuickModelTest.java +++ b/src/test/java/org/orekit/models/earth/ionosphere/NeQuickModelTest.java @@ -266,8 +266,8 @@ private > void doTestFieldAntiMeridian(final F // Date final FieldAbsoluteDate date = - new FieldAbsoluteDate(field, - new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC())); + new FieldAbsoluteDate<>(field, + new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC())); // Geodetic points final FieldGeodeticPoint recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(-31.80)), @@ -281,4 +281,50 @@ private > void doTestFieldAntiMeridian(final F } + @Test + public void testZenith() { + + // Model + final NeQuickModel model = new NeQuickModel(medium); + + // Date + final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC()); + + // Geodetic points + final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(51.678), FastMath.toRadians(-9.448), 0.0); + final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(51.678), FastMath.toRadians(-9.448), 600000.0); + double stec = model.stec(date, recP, satP); + Assertions.assertEquals(26.346, stec, 0.001); + + } + + @Test + public void testFieldZenith() { + doTestFieldZenith(Binary64Field.getInstance()); + } + + private > void doTestFieldZenith(final Field field) { + + final T zero = field.getZero(); + + // Model + final NeQuickModel model = new NeQuickModel(medium); + + // Date + final FieldAbsoluteDate date = + new FieldAbsoluteDate<>(field, + new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC())); + + // Geodetic points + final FieldGeodeticPoint recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(51.678)), + FastMath.toRadians(zero.newInstance(-9.448)), + zero.newInstance(0.0)); + final FieldGeodeticPoint satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(51.678)), + FastMath.toRadians(zero.newInstance(-9.448)), + zero.newInstance(600000.0)); + T stec = model.stec(date, recP, satP); + Assertions.assertEquals(26.346, stec.getReal(), 0.001); + + } + } diff --git a/src/test/java/org/orekit/propagation/FieldStateCovarianceTest.java b/src/test/java/org/orekit/propagation/FieldStateCovarianceTest.java index 900f527c2d..bf65944b9d 100644 --- a/src/test/java/org/orekit/propagation/FieldStateCovarianceTest.java +++ b/src/test/java/org/orekit/propagation/FieldStateCovarianceTest.java @@ -18,6 +18,7 @@ import org.hipparchus.CalculusFieldElement; import org.hipparchus.Field; +import org.hipparchus.FieldElement; import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.hipparchus.linear.BlockFieldMatrix; import org.hipparchus.linear.FieldMatrix; @@ -37,9 +38,11 @@ import org.orekit.frames.LOF; import org.orekit.frames.LOFType; import org.orekit.orbits.FieldCartesianOrbit; +import org.orekit.orbits.FieldKeplerianOrbit; import org.orekit.orbits.FieldOrbit; import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngleType; +import org.orekit.time.AbsoluteDate; import org.orekit.time.FieldAbsoluteDate; import org.orekit.time.TimeScalesFactory; import org.orekit.utils.Constants; @@ -47,6 +50,11 @@ import org.orekit.utils.IERSConventions; import org.orekit.utils.PVCoordinates; +import java.util.Arrays; + +import static org.junit.jupiter.api.Assertions.assertDoesNotThrow; +import static org.junit.jupiter.api.Assertions.assertEquals; + class FieldStateCovarianceTest { private final Field field = Binary64Field.getInstance(); @@ -1068,4 +1076,103 @@ void testConversionToNonFieldEquivalentFrame() { } -} \ No newline at end of file + @Test + public void testIssue1485CartesianOrbitTypeAndNullAngleType() { + // GIVEN + final Field field = Binary64Field.getInstance(); + + // WHEN & THEN + doTestIssue1485CartesianOrbitTypeAndNullAngleType(field); + } + + > void doTestIssue1485CartesianOrbitTypeAndNullAngleType(final Field field) { + // GIVEN + final T one = field.getOne(); + final FieldAbsoluteDate date = new FieldAbsoluteDate<>(field, AbsoluteDate.ARBITRARY_EPOCH); + final FieldMatrix expectedMatrix = MatrixUtils.createFieldIdentityMatrix(field, 6); + final PositionAngleType anomalyType = PositionAngleType.MEAN; + final Frame inFrame = FramesFactory.getEME2000(); + final Frame outFrame = FramesFactory.getGCRF(); + final double mu = Constants.EGM96_EARTH_MU; + + final FieldKeplerianOrbit orbit = + new FieldKeplerianOrbit<>(one.multiply(1e7), one.multiply(0.01), one.multiply(1.), one.multiply(2.), + one.multiply(3.), one.multiply(4.), anomalyType, inFrame, date, one.multiply(mu)); + + final FieldStateCovariance originalCovariance = + new FieldStateCovariance<>(expectedMatrix, date, inFrame, OrbitType.CARTESIAN, null); + + // WHEN & THEN + assertDoesNotThrow(() -> originalCovariance.changeCovarianceFrame(orbit, outFrame)); + + // THEN + final FieldStateCovariance convertedCovariance = originalCovariance.changeCovarianceFrame(orbit, outFrame); + assertEquals(getNorm1(originalCovariance.getMatrix()), getNorm1(convertedCovariance.getMatrix())); + } + + @Test + void testIssue1485SameLOF() { + // GIVEN + final Field field = Binary64Field.getInstance(); + final LOF lof = LOFType.TNW; + + // WHEN & THEN + doTestIssue1485SameFrameDefinition(field, lof); + } + + @Test + void testIssue1485SameFrame() { + // GIVEN + final Field field = Binary64Field.getInstance(); + final Frame eme2000 = FramesFactory.getEME2000(); + + // WHEN & THEN + doTestIssue1485SameFrameDefinition(field, eme2000); + } + + > void doTestIssue1485SameFrameDefinition(Field field, T frameOrLOF) { + // GIVEN + final I one = field.getOne(); + final FieldAbsoluteDate date = new FieldAbsoluteDate<>(field, AbsoluteDate.ARBITRARY_EPOCH); + final FieldMatrix expectedMatrix = MatrixUtils.createFieldIdentityMatrix(field, 6); + final PositionAngleType anomalyType = PositionAngleType.MEAN; + final Frame frame = FramesFactory.getEME2000(); + final double mu = Constants.EGM96_EARTH_MU; + final FieldKeplerianOrbit orbit = + new FieldKeplerianOrbit<>(one.multiply(1e7), one.multiply(0.01), one.multiply(1.), one.multiply(2.), + one.multiply(3.), one.multiply(4.), anomalyType, frame, date, one.multiply(mu)); + // WHEN + final FieldStateCovariance originalCovariance; + final FieldStateCovariance covariance; + if (frameOrLOF instanceof LOF) { + originalCovariance = new FieldStateCovariance<>(expectedMatrix, date, (LOF) frameOrLOF); + covariance = originalCovariance.changeCovarianceFrame(orbit, (LOF) frameOrLOF); + } else { + originalCovariance = + new FieldStateCovariance<>(expectedMatrix, date, (Frame) frameOrLOF, orbit.getType(), anomalyType); + covariance = originalCovariance.changeCovarianceFrame(orbit, (Frame) frameOrLOF); + } + + // THEN + final FieldMatrix actualMatrix = covariance.getMatrix(); + Assertions.assertEquals(0, getNorm1(actualMatrix.subtract(expectedMatrix))); + } + + /** + * Compute L1 norm of given field matrix + * + * @param fieldMatrix field matrix + * @param type of field element + * @return real L1 norm + */ + private > double getNorm1(final FieldMatrix fieldMatrix) { + + return Arrays.stream( + Arrays.stream(fieldMatrix.getData()) + .flatMap(Arrays::stream) + .mapToDouble(FieldElement::getReal) + .toArray()) + .sum(); + } + +} diff --git a/src/test/java/org/orekit/propagation/StateCovarianceTest.java b/src/test/java/org/orekit/propagation/StateCovarianceTest.java index e3fc690f90..ad70dcca90 100644 --- a/src/test/java/org/orekit/propagation/StateCovarianceTest.java +++ b/src/test/java/org/orekit/propagation/StateCovarianceTest.java @@ -32,8 +32,10 @@ import org.orekit.forces.gravity.potential.ICGEMFormatReader; import org.orekit.frames.Frame; import org.orekit.frames.FramesFactory; +import org.orekit.frames.LOF; import org.orekit.frames.LOFType; import org.orekit.orbits.CartesianOrbit; +import org.orekit.orbits.KeplerianOrbit; import org.orekit.orbits.Orbit; import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngleType; @@ -44,6 +46,9 @@ import org.orekit.utils.IERSConventions; import org.orekit.utils.PVCoordinates; +import static org.junit.jupiter.api.Assertions.assertDoesNotThrow; +import static org.junit.jupiter.api.Assertions.assertEquals; + public class StateCovarianceTest { final private double DEFAULT_VALLADO_THRESHOLD = 1e-6; @@ -82,11 +87,11 @@ public void testConstructor() { final AbsoluteDate initialDate = new AbsoluteDate(); final Frame initialInertialFrame = FramesFactory.getGCRF(); final StateCovariance stateCovariance = new StateCovariance(getValladoInitialCovarianceMatrix(), initialDate, initialInertialFrame, OrbitType.CARTESIAN, PositionAngleType.MEAN); - Assertions.assertEquals(OrbitType.CARTESIAN, stateCovariance.getOrbitType()); - Assertions.assertEquals(PositionAngleType.MEAN, stateCovariance.getPositionAngleType()); - Assertions.assertEquals(initialInertialFrame, stateCovariance.getFrame()); + assertEquals(OrbitType.CARTESIAN, stateCovariance.getOrbitType()); + assertEquals(PositionAngleType.MEAN, stateCovariance.getPositionAngleType()); + assertEquals(initialInertialFrame, stateCovariance.getFrame()); Assertions.assertNull(stateCovariance.getLOF()); - Assertions.assertEquals(initialDate, stateCovariance.getDate()); + assertEquals(initialDate, stateCovariance.getDate()); } public void setUp() { @@ -120,11 +125,11 @@ public static void compareCovariance(final RealMatrix reference, final RealMatri for (int row = 0; row < reference.getRowDimension(); row++) { for (int column = 0; column < reference.getColumnDimension(); column++) { if (reference.getEntry(row, column) == 0) { - Assertions.assertEquals(reference.getEntry(row, column), computed.getEntry(row, column), + assertEquals(reference.getEntry(row, column), computed.getEntry(row, column), threshold); } else { - Assertions.assertEquals(reference.getEntry(row, column), computed.getEntry(row, column), + assertEquals(reference.getEntry(row, column), computed.getEntry(row, column), FastMath.abs(threshold * reference.getEntry(row, column))); } } @@ -961,10 +966,77 @@ void testIssue1052() { Assertions.assertDoesNotThrow(() -> propagator.propagate(initialDate.shiftedBy(1))); // Assert that propagated covariance is in the same LOF as the initial covariance - Assertions.assertEquals(LOFType.QSW, propagatedCovariance.getLOF()); + assertEquals(LOFType.QSW, propagatedCovariance.getLOF()); } + @Test + void testIssue1485CartesianOrbitTypeAndNullAngleType() { + // GIVEN + final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH; + final RealMatrix expectedMatrix = MatrixUtils.createRealIdentityMatrix(6); + final PositionAngleType anomalyType = PositionAngleType.MEAN; + final Frame inFrame = FramesFactory.getEME2000(); + final Frame outFrame = FramesFactory.getGCRF(); + final double mu = Constants.EGM96_EARTH_MU; + + final KeplerianOrbit orbit = new KeplerianOrbit(1e7, 0.01, 1., 2., 3., 4., anomalyType, inFrame, date, mu); + + final StateCovariance originalCovariance = + new StateCovariance(expectedMatrix, date, inFrame, OrbitType.CARTESIAN, null); + + // WHEN & THEN + assertDoesNotThrow(()-> originalCovariance.changeCovarianceFrame(orbit, outFrame)); + + // THEN + final StateCovariance convertedCovariance = originalCovariance.changeCovarianceFrame(orbit, outFrame); + assertEquals(originalCovariance.getMatrix().getNorm1(), convertedCovariance.getMatrix().getNorm1(), 1e-15); + } + + @Test + void testIssue1485SameLOF() { + // GIVEN + final LOF lof = LOFType.TNW; + + // WHEN & THEN + doTestIssue1485SameFrameDefinition(lof); + } + + @Test + void testIssue1485SameFrame() { + // GIVEN + final Frame eme2000 = FramesFactory.getEME2000(); + + // WHEN & THEN + doTestIssue1485SameFrameDefinition(eme2000); + } + + void doTestIssue1485SameFrameDefinition(T frameOrLOF) { + // GIVEN + final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH; + final RealMatrix expectedMatrix = MatrixUtils.createRealIdentityMatrix(6); + final PositionAngleType anomalyType = PositionAngleType.MEAN; + final Frame frame = FramesFactory.getEME2000(); + final double mu = Constants.EGM96_EARTH_MU; + + final KeplerianOrbit orbit = new KeplerianOrbit(1e7, 0.01, 1., 2., 3., 4., anomalyType, frame, date, mu); + + // WHEN + final StateCovariance originalCovariance; + final StateCovariance covariance; + if (frameOrLOF instanceof LOF) { + originalCovariance = new StateCovariance(expectedMatrix, date, (LOF) frameOrLOF); + covariance = originalCovariance.changeCovarianceFrame(orbit, (LOF) frameOrLOF); + } else { + originalCovariance = new StateCovariance(expectedMatrix, date, (Frame) frameOrLOF, orbit.getType(), anomalyType); + covariance = originalCovariance.changeCovarianceFrame(orbit, (Frame) frameOrLOF); + } + + // THEN + final RealMatrix actualMatrix = covariance.getMatrix(); + assertEquals(0., actualMatrix.subtract(expectedMatrix).getNorm1()); + } + private NumericalPropagator buildDefaultPropagator(final SpacecraftState state, final OrbitType orbitType, final PositionAngleType positionAngleType) { @@ -991,3 +1063,4 @@ private ODEIntegrator buildDefaultODEIntegrator(final Orbit orbit, final OrbitTy } } + diff --git a/src/test/java/org/orekit/propagation/integration/AbstractIntegratedPropagatorTest.java b/src/test/java/org/orekit/propagation/integration/AbstractIntegratedPropagatorTest.java index b07b7593c4..a8862fa827 100644 --- a/src/test/java/org/orekit/propagation/integration/AbstractIntegratedPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/integration/AbstractIntegratedPropagatorTest.java @@ -89,6 +89,44 @@ public void testIssue1254() { Assertions.assertEquals(0., ephemeris.getMinDate().durationFrom(minDate), 0.); Assertions.assertEquals(0., ephemeris.getMaxDate().durationFrom(maxDate), 0.); } + + /** Test issue 1461. + *

    Test for the new generic method AbstractIntegratedPropagator.reset(SpacecraftState, PropagationType) + */ + @Test + public void testIssue1461() { + // GIVEN + // GEO orbit + final Orbit startOrbit = new EquinoctialOrbit(42165765.0, 0.0, 0.0, 0.0, 0.0, 0.0, PositionAngleType.TRUE, + FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, + Constants.IERS2010_EARTH_MU); + + // Init numerical propagator + final SpacecraftState state = new SpacecraftState(startOrbit); + final NumericalPropagator propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(300.)); + propagator.setInitialState(state); + + // WHEN + propagator.resetInitialState(state, PropagationType.OSCULATING); + final SpacecraftState oscState = propagator.getInitialState(); + + propagator.resetInitialState(state, PropagationType.MEAN); + final SpacecraftState meanState = propagator.getInitialState(); + + // THEN + + // Check that all three states are identical + final double dpOsc = oscState.getPosition().distance(state.getPosition()); + final double dvOsc = oscState.getPVCoordinates().getVelocity().distance(state.getPVCoordinates().getVelocity()); + + final double dpMean = meanState.getPosition().distance(state.getPosition()); + final double dvMean = meanState.getPVCoordinates().getVelocity().distance(state.getPVCoordinates().getVelocity()); + + Assertions.assertEquals(0., dpOsc, 0.); + Assertions.assertEquals(0., dvOsc, 0.); + Assertions.assertEquals(0., dpMean, 0.); + Assertions.assertEquals(0., dvMean, 0.); + } private static class TestAbstractIntegratedPropagator extends AbstractIntegratedPropagator { diff --git a/src/test/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagatorTest.java b/src/test/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagatorTest.java index 1551d929a8..499db2ffc6 100644 --- a/src/test/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagatorTest.java @@ -16,18 +16,30 @@ */ package org.orekit.propagation.integration; +import org.hipparchus.CalculusFieldElement; +import org.hipparchus.Field; import org.hipparchus.complex.Complex; import org.hipparchus.complex.ComplexField; import org.hipparchus.ode.FieldODEIntegrator; +import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator; +import org.hipparchus.util.Binary64Field; import org.junit.jupiter.api.Assertions; import org.junit.jupiter.api.Test; import org.mockito.Mockito; import org.orekit.attitudes.AttitudeProvider; import org.orekit.frames.Frame; +import org.orekit.frames.FramesFactory; +import org.orekit.orbits.EquinoctialOrbit; +import org.orekit.orbits.FieldEquinoctialOrbit; +import org.orekit.orbits.FieldOrbit; import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngleType; +import org.orekit.propagation.FieldSpacecraftState; import org.orekit.propagation.PropagationType; +import org.orekit.propagation.numerical.FieldNumericalPropagator; +import org.orekit.time.AbsoluteDate; import org.orekit.time.FieldAbsoluteDate; +import org.orekit.utils.Constants; class FieldAbstractIntegratedPropagatorTest { @@ -50,6 +62,53 @@ void testGetResetAtEnd(final boolean expectedResetAtEnd) { final boolean actualResetAtEnd = testAbstractIntegratedPropagator.getResetAtEnd(); Assertions.assertEquals(expectedResetAtEnd, actualResetAtEnd); } + + /** Test issue 1461. + *

    Test for the new generic method AbstractIntegratedPropagator.reset(SpacecraftState, PropagationType) + */ + @Test + void testIssue1461() { + doTestIssue1461(Binary64Field.getInstance()); + } + + /** Method for running test for issue 1461. */ + private > void doTestIssue1461(Field field) { + // GIVEN + + final T zero = field.getZero(); + // GEO orbit + final FieldOrbit startOrbit = new FieldEquinoctialOrbit<>(field, + new EquinoctialOrbit(42165765.0, 0.0, 0.0, 0.0, 0.0, 0.0, PositionAngleType.TRUE, + FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, + Constants.IERS2010_EARTH_MU)); + + // Init numerical propagator + final FieldSpacecraftState state = new FieldSpacecraftState<>(startOrbit); + final FieldNumericalPropagator propagator = new FieldNumericalPropagator<>(field, + new ClassicalRungeKuttaFieldIntegrator(field, zero.newInstance(300.))); + propagator.setInitialState(state); + + // WHEN + propagator.resetInitialState(state, PropagationType.OSCULATING); + final FieldSpacecraftState oscState = propagator.getInitialState(); + + propagator.resetInitialState(state, PropagationType.MEAN); + final FieldSpacecraftState meanState = propagator.getInitialState(); + + // THEN + + // Check that all three states are identical + final double dpOsc = oscState.getPosition().distance(state.getPosition()).getReal(); + final double dvOsc = oscState.getPVCoordinates().getVelocity().distance(state.getPVCoordinates().getVelocity()).getReal(); + + final double dpMean = meanState.getPosition().distance(state.getPosition()).getReal(); + final double dvMean = meanState.getPVCoordinates().getVelocity().distance(state.getPVCoordinates().getVelocity()).getReal(); + + Assertions.assertEquals(0., dpOsc, 0.); + Assertions.assertEquals(0., dvOsc, 0.); + Assertions.assertEquals(0., dpMean, 0.); + Assertions.assertEquals(0., dvMean, 0.); + } private static class TestFieldAbstractIntegratedPropagator extends FieldAbstractIntegratedPropagator { diff --git a/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java b/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java index 8e4cc05df9..b9d9800220 100644 --- a/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java @@ -16,6 +16,10 @@ */ package org.orekit.propagation.semianalytical.dsst; +import static org.hamcrest.Matchers.contains; +import static org.hamcrest.Matchers.is; +import static org.hamcrest.Matchers.nullValue; + import java.io.IOException; import java.text.ParseException; import java.util.ArrayList; @@ -50,9 +54,9 @@ import org.junit.jupiter.api.BeforeEach; import org.junit.jupiter.api.Test; import org.orekit.OrekitMatchers; -import org.orekit.TestUtils; import org.orekit.Utils; import org.orekit.attitudes.AttitudeProvider; +import org.orekit.attitudes.FrameAlignedProvider; import org.orekit.attitudes.LofOffset; import org.orekit.bodies.CelestialBody; import org.orekit.bodies.CelestialBodyFactory; @@ -117,10 +121,6 @@ import org.orekit.utils.TimeSpanMap; import org.orekit.utils.TimeStampedPVCoordinates; -import static org.hamcrest.Matchers.contains; -import static org.hamcrest.Matchers.is; -import static org.hamcrest.Matchers.nullValue; - public class DSSTPropagatorTest { private DSSTPropagator dsstProp; @@ -1197,6 +1197,61 @@ public void testIssue672() { Assertions.assertNotNull(state); } + /** Test issue 1461. + *

    Test for the dedicated method DSSTPropagator.reset(SpacecraftState, PropagationType) + *

    This should change the status of attribute "initialIsOsculating" depending on input PropagationType + */ + @Test + public void testIssue1461() { + + // GIVEN + // ----- + + final double step = 60.; + final int nStep = 10; + final Frame gcrf = FramesFactory.getGCRF(); + final AttitudeProvider attitude = new FrameAlignedProvider(gcrf); + + // LEO orbit as input + final Orbit initialOrbit = new KeplerianOrbit(7000e3, 1.e-3, FastMath.toRadians(98.), 0., 0., 0., + PositionAngleType.ECCENTRIC, gcrf, + AbsoluteDate.ARBITRARY_EPOCH, Constants.EGM96_EARTH_MU); + final AbsoluteDate t0 = initialOrbit.getDate(); + final SpacecraftState initialState = new SpacecraftState(initialOrbit); + final ClassicalRungeKuttaIntegrator integrator = new ClassicalRungeKuttaIntegrator(step); + + // DSST is configured to propagate in MEAN elements + final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.MEAN); + + // The initial state is OSCULATING + propagator.setInitialState(initialState, PropagationType.OSCULATING); + + // Using a J2-only perturbed force model, the SMA of the mean / averaged orbit should remained constant during propagation + propagator.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0))); + propagator.setAttitudeProvider(attitude); + + // Initial mean SMA + final double initialMeanSma = DSSTPropagator.computeMeanState(initialState, + attitude, + propagator.getAllForceModels()).getA(); + // WHEN + // ---- + + // Propagating step by step and getting mean SMAs + final List propagatedMeanSma = new ArrayList<>(); + for (int i = 1; i <= nStep; i++) { + propagatedMeanSma.add(propagator.propagate(t0.shiftedBy(i * step)).getA()); + } + + // THEN + // ---- + + // Mean SMA should remain constant and equal to initial mean sma + for (Double meanSma : propagatedMeanSma) { + Assertions.assertEquals(initialMeanSma, meanSma, 0.); + } + } + /** * Check that the DSST can include the derivatives of force model parameters in the * Jacobian. Uses a very basic force model for a quick test. Checks both mean and @@ -1488,6 +1543,7 @@ public void updateShortPeriodTerms(double[] parameters, SpacecraftState... meanStates) { } + @SuppressWarnings("unchecked") @Override public > void updateShortPeriodTerms( T[] parameters, FieldSpacecraftState... meanStates) { diff --git a/src/test/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagatorTest.java b/src/test/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagatorTest.java index d0334582a5..f6c6b4eda7 100644 --- a/src/test/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/semianalytical/dsst/FieldDSSTPropagatorTest.java @@ -49,6 +49,7 @@ import org.orekit.OrekitMatchers; import org.orekit.Utils; import org.orekit.attitudes.AttitudeProvider; +import org.orekit.attitudes.FrameAlignedProvider; import org.orekit.attitudes.LofOffset; import org.orekit.bodies.CelestialBody; import org.orekit.bodies.CelestialBodyFactory; @@ -73,6 +74,7 @@ import org.orekit.orbits.FieldEquinoctialOrbit; import org.orekit.orbits.FieldKeplerianOrbit; import org.orekit.orbits.FieldOrbit; +import org.orekit.orbits.KeplerianOrbit; import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngleType; import org.orekit.propagation.FieldBoundedPropagator; @@ -1279,6 +1281,69 @@ private > void doTestIssue672(final Field f Assertions.assertNotNull(state); } + /** Test issue 1461. + *

    Test for the dedicated method DSSTPropagator.reset(SpacecraftState, PropagationType) + *

    This should change the status of attribute "initialIsOsculating" depending on input PropagationType + */ + @Test + public void testIssue1461() { + doTestIssue1461(Binary64Field.getInstance()); + } + + /** Method for running test for issue 1461. */ + private > void doTestIssue1461(final Field field) { + + // GIVEN + // ----- + + final T zero = field.getZero(); + + final double step = 60.; + final int nStep = 10; + final Frame gcrf = FramesFactory.getGCRF(); + final AttitudeProvider attitude = new FrameAlignedProvider(gcrf); + + // LEO orbit as input + final FieldOrbit initialOrbit = new FieldKeplerianOrbit<>(field, + new KeplerianOrbit(7000e3, 1.e-3, FastMath.toRadians(98.), 0., 0., 0., + PositionAngleType.ECCENTRIC, gcrf, + AbsoluteDate.ARBITRARY_EPOCH, Constants.EGM96_EARTH_MU)); + final FieldAbsoluteDate t0 = initialOrbit.getDate(); + final FieldSpacecraftState initialState = new FieldSpacecraftState<>(initialOrbit); + final ClassicalRungeKuttaFieldIntegrator integrator = new ClassicalRungeKuttaFieldIntegrator<>(field, zero.newInstance(step)); + + // DSST is configured to propagate in MEAN elements + final FieldDSSTPropagator propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.MEAN); + + // The initial state is OSCULATING + propagator.setInitialState(initialState, PropagationType.OSCULATING); + + // Using a J2-only perturbed force model, the SMA of the mean / averaged orbit should remained constant during propagation + propagator.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0))); + propagator.setAttitudeProvider(attitude); + + // Initial mean SMA + final double initialMeanSma = FieldDSSTPropagator.computeMeanState(initialState, + attitude, + propagator.getAllForceModels()).getA().getReal(); + // WHEN + // ---- + + // Propagating step by step and getting mean SMAs + final List propagatedMeanSma = new ArrayList<>(); + for (int i = 1; i <= nStep; i++) { + propagatedMeanSma.add(propagator.propagate(t0.shiftedBy(i * step)).getA().getReal()); + } + + // THEN + // ---- + + // Mean SMA should remain constant and equal to initial mean sma + for (Double meanSma : propagatedMeanSma) { + Assertions.assertEquals(initialMeanSma, meanSma, 0.); + } + } + private > FieldSpacecraftState getGEOState(final Field field) { final T zero = field.getZero(); // No shadow at this date diff --git a/src/test/java/org/orekit/time/AbsoluteDateTest.java b/src/test/java/org/orekit/time/AbsoluteDateTest.java index c26ee666a0..d8792fc152 100644 --- a/src/test/java/org/orekit/time/AbsoluteDateTest.java +++ b/src/test/java/org/orekit/time/AbsoluteDateTest.java @@ -1328,10 +1328,10 @@ public void testToString() { // test proleptic checkToString(new AbsoluteDate(123, 4, 5, 6, 7, 8.9, utc), "0123-04-05T06:07:08.900"); - // there is not way to produce valid RFC3339 for these cases + // there is no way to produce valid RFC3339 for these cases // I would rather print something useful than throw an exception // so these cases don't check for a correct answer, just an informative one - checkToString(new AbsoluteDate(-123, 4, 5, 6, 7, 8.9, utc), "-0123-04-05T06:07:08.900"); + checkToString(new AbsoluteDate(-123, 4, 5, 6, 7, 8.9, utc), "-123-04-05T06:07:08.900"); checkToString(new AbsoluteDate(-1230, 4, 5, 6, 7, 8.9, utc), "-1230-04-05T06:07:08.900"); // test far future checkToString(new AbsoluteDate(12300, 4, 5, 6, 7, 8.9, utc), "12300-04-05T06:07:08.900");