Hello orekit community,
I met some issues when I tried to do orbit determination with ground-based optical measurements and inter-satellite range measurements. I called the objects of orekit from Matlab. The main code is as follows:
clear
clc
import org.orekit.data.DataProvidersManager
import org.orekit.data.DirectoryCrawler
import org.orekit.frames.FramesFactory
import org.orekit.time.TimeScalesFactory
import org.orekit.time.AbsoluteDate;
import org.hipparchus.geometry.euclidean.threed.Vector3D
import org.hipparchus.util.FastMath
import org.orekit.utils.IERSConventions
import org.orekit.utils.TimeStampedPVCoordinates
import org.orekit.utils.PVCoordinates;
import org.orekit.models.earth.ReferenceEllipsoid
import org.orekit.bodies.GeodeticPoint
import org.orekit.estimation.measurements.Range
import org.orekit.estimation.measurements.AngularRaDec
import org.orekit.estimation.measurements.ObservableSatellite
import org.orekit.estimation.measurements.AngularAzEl
import org.orekit.estimation.measurements.Position
import org.orekit.frames.TopocentricFrame
import org.orekit.estimation.measurements.GroundStation
import org.orekit.estimation.measurements.InterSatellitesRange
import org.orekit.estimation.iod.IodLaplace
import org.orekit.estimation.iod.IodGibbs
import org.orekit.utils.Constants
import org.orekit.propagation.analytical.KeplerianPropagator
import org.orekit.attitudes.NadirPointing
import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder
import org.orekit.propagation.conversion.NumericalPropagatorBuilder
import org.orekit.orbits.PositionAngle
import org.orekit.orbits.CartesianOrbit
import org.orekit.orbits.OrbitType;
import org.orekit.bodies.CelestialBodyFactory
import org.orekit.forces.gravity.potential.GravityFieldFactory
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel
import org.orekit.forces.gravity.ThirdBodyAttraction
import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient
import org.orekit.forces.radiation.SolarRadiationPressure
import org.orekit.forces.gravity.Relativity
import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation
import org.orekit.models.earth.atmosphere.NRLMSISE00
import org.orekit.models.earth.atmosphere.HarrisPriester
import org.orekit.forces.drag.IsotropicDrag
import org.orekit.forces.drag.DragForce
import org.hipparchus.linear.QRDecomposer
import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer
import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer
import org.orekit.estimation.leastsquares.BatchLSEstimator
import org.orekit.estimation.sequential.KalmanEstimatorBuilder
import org.orekit.propagation.numerical.NumericalPropagator
import java.io.File;
import org.orekit.estimation.measurements.generation.Generator
%% Configure Orekit.
DM = org.orekit.data.DataContext.getDefault().getDataProvidersManager();
dir_orekit_data = 'D:\Codes_Tao\orekit\orekit_11.3.2\orekit-data';
crawler = org.orekit.data.DirectoryCrawler(File(dir_orekit_data));
DM.clearProviders()
DM.addProvider(crawler)
%% Defines frames
gcrf = FramesFactory.getGCRF();
icrf = FramesFactory.getICRF();
eci = FramesFactory.getEME2000();
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
ecef = itrf;
wgs84ellipsoid = ReferenceEllipsoid.getWgs84(itrf);
utc = TimeScalesFactory.getUTC();
mu = 3.986004415e+14; % gravitation coefficient
%% observations
load SimuObsDataForIOD_LEO.mat
i = 9;
MJD1 = SimuObsData(i).MJD;
RA1 = SimuObsData(i).RA;
DE1 = SimuObsData(i).DE;
StateTrueKep = SimuObsData(i).SatStateTrue;
[r, v] = orb2eci(398600.4418, StateTrueKep);
StateTrue = [r; v];
i = 10;
MJD2 = SimuObsData(i).MJD;
RA2 = SimuObsData(i).RA;
DE2 = SimuObsData(i).DE;
MJD = [MJD1; MJD2];
RA = [RA1; RA2];
DE = [DE1; DE2];
%% Ground station (radar)
r2d = 180/pi;
d2r = pi/180;
StationLBH = [117.5745, 40.3959, 960];
sigma_azi = 0.1 * d2r; % rad
sigma_elv = 0.1 * d2r; % rad
sigma_rho = 0.03 * 1000; % m
sigma_ra = 1/3600 * d2r; % rad
sigma_de = 1/3600 * d2r; % rad
weight_rho = 1.0;
weight_azi = 1.0;
weight_elv = 1.0;
weight_ra = 1.0;
weight_de = 1.0;
% noise
RA_noise = sigma_ra * randn(length(MJD), 1);
DE_noise = sigma_de * randn(length(MJD), 1);
LargeRaErrFlag = 1;
if LargeRaErrFlag
RA_noise = RA_noise ./ cosd(DE);
end
RA = RA + RA_noise;
DE = DE + DE_noise;
longitude = FastMath.toRadians(StationLBH(1));
latitude = FastMath.toRadians(StationLBH(2));
altitude = StationLBH(3);
geodeticP = GeodeticPoint(latitude, longitude, altitude);
topocentricFrame = TopocentricFrame(wgs84ellipsoid, geodeticP, 'station');
groundStation = GroundStation(topocentricFrame);
fac = Facility(StationLBH(2), StationLBH(1), StationLBH(3)/1000);
fac.sigma_azi = sigma_azi;
fac.sigma_elv = sigma_elv;
fac.sigma_range = sigma_rho;
%% IOD
mjd_od = MJD1(1); % OD eopch
[x0, P0] = IOD_LS(MJD1, fac.positionECEF, RA1, DE1, sigma_ra*r2d, sigma_de*r2d, LargeRaErrFlag);
fprintf('initial state: [%.3f %.3f %.3f %.3f %.3f %.3f] km & km/sec\n', x0/1000)
fprintf('initial cov (diag): [%.3f %.3f %.3f %.3f %.3f %.3f] km^2 & km^2/s^2\n', diag(P0)/1000)
[year, month, day, hr, min, sec] = invjday(mjd_od + 2400000.5);
initialDate = AbsoluteDate(year, month, day, hr, min, sec, utc);
posisat = Vector3D(x0(1), x0(2), x0(3));
velosat = Vector3D(x0(4), x0(5), x0(6));
pvsat = PVCoordinates(posisat, velosat);
initialOrbit = CartesianOrbit(pvsat, eci, initialDate, mu);
%% Estimator and propagator parameters
% Estimator parameters
estimator_position_scale = 1.0; % m
estimator_convergence_thres = 1e-2;
estimator_max_iterations = 50;
estimator_max_evaluations = 50;
% Orbit propagator parameters
prop_min_step = 0.001; % s
prop_max_step = 300.0; % s
prop_position_error = 10.0; % m
%% Setting up a numerical propagator.
nadirPointing = NadirPointing(gcrf, wgs84ellipsoid);
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error);
propagatorBuilder = NumericalPropagatorBuilder(initialOrbit,...
integratorBuilder, PositionAngle.MEAN, estimator_position_scale, nadirPointing);
propagatorBuilder.setAttitudeProvider(nadirPointing);
% Setting force model
propagatorBuilder.setMass(1000);
% Earth gravity field with degree 64 and order 64
gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(21, 21);
gravityAttractionModel = HolmesFeatherstoneAttractionModel(ecef, gravityProvider);
propagatorBuilder.addForceModel(gravityAttractionModel);
% Moon and Sun perturbations
moon = CelestialBodyFactory.getMoon();
sun = CelestialBodyFactory.getSun();
moon_3dbodyattraction = ThirdBodyAttraction(moon);
propagatorBuilder.addForceModel(moon_3dbodyattraction);
sun_3dbodyattraction = ThirdBodyAttraction(sun);
propagatorBuilder.addForceModel(sun_3dbodyattraction);
% Solar radiation pressure
isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(10, 1.5);
solarRadiationPressure = SolarRadiationPressure(sun, wgs84ellipsoid.getEquatorialRadius(), isotropicRadiationSingleCoeff);
propagatorBuilder.addForceModel(solarRadiationPressure);
% Atmospheric drag
AVERAGE = javaMethod('valueOf','org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation$StrengthLevel','AVERAGE');
msafe = MarshallSolarActivityFutureEstimation(...
MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,...
AVERAGE);
atmosphere = NRLMSISE00(msafe, sun, wgs84ellipsoid);
isotropicDrag = IsotropicDrag(forceParam.areaDrag, forceParam.cd);
dragForce = DragForce(atmosphere, isotropicDrag);
propagatorBuilder.addForceModel(dragForce)
%% inter satellite range measurement
load('SimuObsDataInterSatelliteRange.mat')
initialStateLocal = SimuObsDataInterSatelliteRange.initialState2;
forceParam = SimuObsDataInterSatelliteRange.forceParam;
builderLocalSat = NumericalPropagatorBuilder(initialStateLocal.getOrbit,...
integratorBuilder, PositionAngle.MEAN, estimator_position_scale, nadirPointing);
builderLocalSat.setMass(700);
builderLocalSat.addForceModel(gravityAttractionModel);
builderLocalSat.addForceModel(moon_3dbodyattraction);
builderLocalSat.addForceModel(sun_3dbodyattraction);
generator = Generator();
satLocal = generator.addPropagator(builderLocalSat.buildPropagator(builderLocalSat.getSelectedNormalizedParameters()));
%% Setting up the estimator
matrixDecomposer = QRDecomposer(1e-11);
optimizer = GaussNewtonOptimizer(matrixDecomposer, false);
propagatorBuilderArray = javaArray('org.orekit.propagation.conversion.NumericalPropagatorBuilder', 1);
propagatorBuilderArray(1) = propagatorBuilder;
estimator = BatchLSEstimator(optimizer, propagatorBuilderArray);
estimator.setParametersConvergenceThreshold(estimator_convergence_thres);
estimator.setMaxIterations(estimator_max_iterations);
estimator.setMaxEvaluations(estimator_max_evaluations);
% Fetching optical measurement data
durations = (MJD - mjd_od) * 86400;
observableSatellite = ObservableSatellite(0);
for i = 1 : length(MJD)
orekitRaDec = AngularRaDec(groundStation, gcrf, initialDate.shiftedBy(durations(i)),...
[RA(i),DE(i)]*d2r, [sigma_ra,sigma_de], [weight_azi,weight_elv], observableSatellite);
estimator.addMeasurement(orekitRaDec);
end
% Fetching interSatelliteRange measurement data
MJD_range = SimuObsDataInterSatelliteRange.MJD;
for i = 10% : length(MJD_range)
[year, month, day, hr, min, sec] = invjday(MJD_range(i) + 2400000.5);
thisDate = AbsoluteDate(year, month, day, hr, min, sec, utc);
intersatellitesrange = InterSatellitesRange(satLocal, observableSatellite,...
false,thisDate,SimuObsDataInterSatelliteRange.interSatelliteRange(i),...
10.0,1.0);
estimator.addMeasurement(intersatellitesrange);
end
propagatorBuilderArray(2) = builderLocalSat; % propagator of the local sat of the interSatelliteRange measurement
%% Estimate the orbit
estimatedPropagatorArray = estimator.estimate();
With only optical measurements, the code ran well and the orbit determination result is okay. However, once I tried to add inter-satellite range measurements, it run error, with the following information:
Java exception occurred:
org.orekit.errors.OrekitException: no solar activity available at 2020-07-08T18:21:�Z, data available only in range [1998-09-01T00:00:00.000Z, 2041-10-01T00:00:00.000Z]
at org.orekit.models.earth.atmosphere.NRLMSISE00.getDensity(NRLMSISE00.java:1158)
at org.orekit.forces.drag.DragForce.acceleration(DragForce.java:80)
at org.orekit.forces.ForceModel.addContribution(ForceModel.java:109)
at org.orekit.propagation.numerical.NumericalPropagator$Main.computeDerivatives(NumericalPropagator.java:899)
at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:758)
at org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:134)
at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(AbstractIntegrator.java:265)
at org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator.initializeStep(AdaptiveStepsizeIntegrator.java:226)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:226)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:477)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:408)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:140)
at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(AbstractBatchLSModel.java:319)
at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:615)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:163)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:435)
Why is the solar activity data not available since the observation time is on 2020-07-08? By the way, the time was shown as ’2020-07-08T18:21:�Z‘, and I have know idea why the second information is missed.
To avoid this error, I removed the drag force from the code:
% propagatorBuilder.addForceModel(dragForce)
but Matlab show another error
Java exception occurred:
org.orekit.errors.OrekitException: NaN appears during integration near time �
at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:511)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:408)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:140)
at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(AbstractBatchLSModel.java:319)
at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:615)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:163)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:435)
Caused by: org.hipparchus.exception.MathIllegalStateException: NaN appears during integration near time �
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:269)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:477)
... 7 more
I’m confused of these results. It seems that something went wrong about ‘time’, but I cannot figure it out now. Do you have any idea?
Thank you very much for any help.
Xuefeng Tao.