Point Inside Ellipsoid Error when Using LevenbergMarquardtOptimizer and BatchLSEstimator, but Data Look Good

rhooper · November 1, 2024, 2:49pm

Hello. When using the BatchLSEstimator with the LevenbergMarquardtOptimizer, I get the org.orekit.errors.OrekitException: point is inside ellipsoid exception after 23 iterations. The data I am using is simulated observation data from 3 ObservationStations, where I pulled the azimuth, elevation, and range if the satellite was visible to that station at each time step. The data is organized into tracks, where each track is made up of the consecutive data from each station that could see the satellite at each time step, as long as at least one station could see the satellite. I then took that data and created AngularAzEl and Range measurements.

The satellite I am simulating is the ENVISAT which has an orbital altitude of about 770km. I dug into the data using the TopocentricFrame.pointAtDistance(azimuth,elevation,range).getAltitude() and in MATLAB using their aer2ned, then ned2ecef functions and both found that the altitude stays around 770km.

Here are the parameters I used to make the stations:

if name == "Eglin":
    self._latitude = FastMath.toRadians(30.57) # [rad]
    self._longitude = FastMath.toRadians(-86.21) # [rad]
    self._altitude = 34.7 # [m]
    self._minElevation = FastMath.toRadians(1.0) # Minimum viewing angle
    self._maxRange = 13210.0e5 # Maximum observation distance, [m]
    self._azSigma = FastMath.toRadians(0.0154) # [rad]
    self._elSigma = FastMath.toRadians(0.0147) # [rad]
    self._rangeSigma = 32.1 # [m]
elif name == "Kaena Point":
    self._latitude = FastMath.toRadians(21.57) # [rad]
    self._longitude = FastMath.toRadians(-158.27) # [rad]
    self._altitude = 300.2 # [m]
    self._minElevation = FastMath.toRadians(0.0) # Minimum viewing angle
    self._maxRange = 6380.0e3 # Maximum observation distance, [m]
    self._azSigma = FastMath.toRadians(0.0224) # [rad]
    self._elSigma = FastMath.toRadians(0.0139) # [rad]
    self._rangeSigma = 92.5 # [m]
elif name == "Clear":
    self._latitude = FastMath.toRadians(64.29) # [rad]
    self._longitude = FastMath.toRadians(-149.19) # [rad]
    self._altitude = 213.3 # [m]
    self._minElevation = FastMath.toRadians(1.0) # Minimum viewing angle
    self._maxRange = 4910.0e3 # Maximum observation distance, [m]
    self._azSigma = FastMath.toRadians(0.0791) # [rad]
    self._elSigma = FastMath.toRadians(0.0240) # [rad]
    self._rangeSigma = 62.5 # [m]

I am running the estimation on 12 hours worth of tracks and only if that 12 hours has at least 100 measurements (50 AngularAzEl and 50 Range) in it.

I’ve done the estimation with and without simulated random error and got the same results. Can someone help me answer why I am getting the inside ellipsoid exception?

Here are the relevant blocks of code:

self.observableSatellite = ObservableSatellite(0)
# Create measurements
# Read azimuth and elevation, then add error
self._tracks[indexMeasurement].measurements.append(AngularAzEl(
                station._groundStationForEstimation,
                AbsoluteDate(row["date"],Common.utc()),
                [row["azimuth[rad]"]+random.gauss(0.0,1.0)*station._azSigma*sigmaErrorOn,row["elevation[rad]"]+random.gauss(0.0,1.0)*station._elSigma*sigmaErrorOn],
                [station._azSigma,station._elSigma],
                [1.0,1.0],
                self.observableSatellite
            ))

# Read range and add error
two_way = False
self._tracks[indexMeasurement].measurements.append(Range(
                station._groundStationForEstimation,
                two_way,
                AbsoluteDate(row["date"],Common.utc()),
                row["range[m]"]*random.gauss(0.0,1.0)*station._rangeSigma*sigmaErrorOn,
                station._rangeSigma,
                1.0,
                self.observableSatellite
            ))
dfEstimate = self.estimationOneState(track,dfEstimate)

def estimationOneState(self, track: 'Track', dfEstimate: pd.DataFrame):
    RSO = self._RSO
    initialOrbitGuess = CartesianOrbit(track.timeStampedTruthPVs[0],Common.inertialFrame,Common.muEarth) # Take the first time of each track as the moment of observation
    initialStateGuess = SpacecraftState(initialOrbitGuess,self._RSO.mass)

    # Use the same dynamics model, but with different initial dragCoefficient and reflectionCoeffient
    estimationPropagatorBuilder = Common.propagatorBuilder(
        initialStateGuess,RSO.area,RSO.dragCoefficientEstimationAndPrediction,RSO.reflectionCoefficientEstimationAndPrediction,self._atmosphereModel
    )

    # Construct Batch Least Squares Estimator
    optimizer = LevenbergMarquardtOptimizer()
    estimator = BatchLSEstimator(optimizer,[estimationPropagatorBuilder])
    estimator.setParametersConvergenceThreshold(self._ESTIMATION_TOLERANCE)
    estimator.setMaxIterations(self._MAXITERATIONS)
    estimator.setMaxEvaluations(self._MAXEVALUATIONS)

    # Add other supported estimated variables
    driver: ParameterDriver
    for driver in estimator.getPropagatorParametersDrivers(False).getDrivers():
        if driver.getName() == "drag coefficient":
            driver.setSelected(self._estimateDragCoefficient)
        if driver.getName() == "reflection coefficient":
            driver.setSelected(self._estimateReflectionCoefficient)
        if driver.getName() == "cross section":
            driver.setSelected(False)

    # Add the measurement observations (including the error) within the previous estimationInterval of the current track to the estimator.
    measurementNumber: int = 0
    if self._estimationInterval <= 0:
        previousTracks: List['SimulateEstimation.Track'] = self._tracks[0:self._tracks.index(track)+1]

        # When the length of estimationInterval is reached, the calculation will start again
        if previousTracks[0].timeStampedTruthPVs[0].getDate().compareTo(
            track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)) <= 0:
            oneTrack: SimulateEstimation.Track
            for oneTrack in previousTracks:
                if oneTrack.timeStampedTruthPVs[len(oneTrack.timeStampedTruthPVs)-1].getDate().compareTo(track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)) < 0:
                    continue
                oneMeasurement: ObservedMeasurement
                for oneMeasurement in oneTrack.measurements:
                    estimator.addMeasurement(oneMeasurement)
                    measurementNumber += 1
    elif self._estimationInterval > 0:
        laterTracks: List[SimulateEstimation.Track] = self._tracks[self._tracks.index(track):len(self._tracks)]
        
        # When the length of estimationInterval is reached, the calculation will start again
        if laterTracks[-1].timeStampedTruthPVs[len(laterTracks[-1].timeStampedTruthPVs)-1].getDate().compareTo(
            track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)
        ) >= 0:
            oneTrack: SimulateEstimation.Track
            for oneTrack in laterTracks:
                if oneTrack.timeStampedTruthPVs[0].getDate().compareTo(track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)) > 0:
                    break
                oneMeasurement: ObservedMeasurement
                for oneMeasurement in oneTrack.measurements:
                    estimator.addMeasurement(oneMeasurement)
                    measurementNumber += 1

    # If there are more than the maximum number of observations, then integration begins
    if measurementNumber >= self._minTrackMeasurement:
        tempLogString = ""

        # Run an estimate
        ######## Orekit 12.0 returns a Propagator list for estimate(), so I'm guessing tat grabbing [0] is correct ########
        estimated: Orbit = estimator.estimate()[0].getInitialState().getOrbit()

        # Output data - pt = position truth, pe = position estimation, dp = difference of position 
        pt: Vector3D
        vt: Vector3D
        pe: Vector3D
        ve: Vector3D
        dp: Vector3D
        dv: Vector3D

        tempEstimateData = []
        
        # Set truth values
        pt = track.timeStampedTruthPVs[0].getPosition()
        vt = track.timeStampedTruthPVs[0].getVelocity()
        tempEstimateData = tempEstimateData + [pt.getX(),pt.getY(),pt.getZ()]
        tempEstimateData = tempEstimateData + [vt.getX(),vt.getY(),vt.getZ()]

        # Set estimation values
        pe = estimated.getPVCoordinates().getPosition()
        ve = estimated.getPVCoordinates().getVelocity()
        tempEstimateData = tempEstimateData + [pe.getX(),pe.getY(),pe.getZ()]
        tempEstimateData = tempEstimateData + [ve.getX(),ve.getY(),ve.getZ()]

        # Set truth minus estimation values
        dp = pt.add(-1,pe)
        dv = vt.add(-1,ve)
        tempEstimateData = tempEstimateData + [dp.getX(),dp.getY(),dp.getZ()]
        tempEstimateData = tempEstimateData + [dv.getX(),dv.getY(),dv.getZ()]

        # The azimuth corresponding to the maximum elevation moment on the current track, elevation, range [rad, m]
        maxElevation = 0
        maxElevationIndex = 0
        m: ObservedMeasurement
        for m in track.measurements:
            if len(m.getObservedValue()) == 1:
                continue
            if maxElevation < m.getObservedValue()[1]:
                maxElevation = m.getObservedValue()[1]
                maxElevationIndex = track.measurements.index(m)

        tempEstimateData = tempEstimateData + [track.measurements[maxElevationIndex].getObservedValue()[0],
                                                maxElevation,
                                                track.measurements[maxElevationIndex+1].getObservedValue()[0]]
        
        # Time + index
        tempEstimateData = tempEstimateData + [f"{estimated.getDate()}",f"{self._tracks.index(track)}"]

        # Cd_true + Cd_estimated + Cr_true + Cr_estimated
        pd: ParameterDriver
        for pd in estimator.getPropagatorParametersDrivers(False).getDrivers():
            if pd.getName() == "drag coefficient":
                tempEstimateData = tempEstimateData + [RSO.dragCoefficient,pd.getValue()]
            if pd.getName() == "reflection coefficient":
                tempEstimateData = tempEstimateData + [RSO.reflectionCoefficient,pd.getValue()]
        
        # Add row to output DataFrame
        dfEstimate.loc[len(dfEstimate)] = tempEstimateData

        # Output Screen
        tempLogString = tempLogString + f"Result on track {self._tracks.index(track)}:\n"
        tempLogString = tempLogString + f"{measurementNumber} measurements on track {self._tracks.index(track)}.\n"
        tempLogString = tempLogString + f"Truth:     \n"
        tempLogString = tempLogString + f"{track.timeStampedTruthPVs[0]}\n"
        tempLogString = tempLogString + f"drag coefficient = {RSO.dragCoefficient}\n"
        tempLogString = tempLogString + f"Guess:     \n"
        tempLogString = tempLogString + f"{initialOrbitGuess.getPVCoordinates()}\n"
        tempLogString = tempLogString + f"drag coefficient = {RSO.dragCoefficientEstimationAndPrediction}\n"
        tempLogString = tempLogString + f"Estimated: \n"
        tempLogString = tempLogString + f"{estimated.getPVCoordinates()}\n"
        pd: ParameterDriver
        for pd in estimator.getPropagatorParametersDrivers(True).getDrivers():
            tempLogString = tempLogString + f"{pd.getName()} = {pd.getValue()}\n"
        tempLogString = tempLogString + f"\n{estimator.getIterationsCount()} iterations, and {estimator.getEvaluationsCount()} evaluations.\n"
        
        # Output log: Output the same results as on the screen to a log file
        tempLogString = tempLogString + f"{estimated.getDate().durationFrom(RSO.initialDate)} | {dp.getNorm()} | {dp.getNorm()/pt.getNorm()*100}% | {dv.getNorm()} | {dv.getNorm()/vt.getNorm()*100}%.\n"
        print(tempLogString)
        with open(self._estimationLogFileName,'a') as f:
            f.write(tempLogString)
    else:
        tempOut = f"# Track {self._tracks.index(track)} is not estimated, because there are not enough measurements on this track."
        print(tempOut)

    return dfEstimate
class Track:
    """
    Stores the observation data of each track, including the PV of each point and the corresponding 
    measurements
    """
    def __init__(self):
        self.timeStampedTruthPVs: List[TimeStampedPVCoordinates] = []
        self.stationIDs: List[int] = []
        self.measurements: List[ObservedMeasurement] = []

MaximeJ · November 4, 2024, 7:33pm

Hi @rhooper,

Sorry for the delay.

It is an unusual exception for the orbit determination, and 23 iterations is a lot so there’s definitely something going wrong here.
The exception comes from the solar radiation pressure model. Could you try deactivating it and see what you get?
It seems like your satellite is heading for reentry though. It would be a good thing for Envisat but I doubt that it is true.

Guess it must be the angles then.
Some pointers that may help:

Everything in Orekit is in SI units, so angles are in radians. But given the name of the angles’ columns in your dictionary of measurements, I doubt it is the issue,
Azimuth are counted clockwise in Orekit from the North direction of the topocentric frame (see here)

Also, are you aware that Orekit has built-in functions for measurements’ generation?
They are all located in the package estimation.measurements.generation and can help you produce the simulated data (with or without noise).
There are many examples in the corresponding test package and also in the class MeasurementGenerator which is called by the tutorial PerformanceTesting.

Hope this helps,
Maxime

rhooper · November 4, 2024, 8:26pm

Hello Maxime,

Thank you for the response. Here is what I get when I turn off the solar radiation pressure model:

org.orekit.errors.OrekitException: Infinite value appears during computation of atmospheric density in NRLMSISE00 model
	at org.orekit.models.earth.atmosphere.NRLMSISE00$Output.densu(NRLMSISE00.java:2622)
	at org.orekit.models.earth.atmosphere.NRLMSISE00$Output.gts7(NRLMSISE00.java:1570)
	at org.orekit.models.earth.atmosphere.NRLMSISE00$Output.gtd7(NRLMSISE00.java:1820)
	at org.orekit.models.earth.atmosphere.NRLMSISE00$Output.gtd7d(NRLMSISE00.java:1928)
	at org.orekit.models.earth.atmosphere.NRLMSISE00.getDensity(NRLMSISE00.java:1181)
	at org.orekit.forces.drag.DragForce.acceleration(DragForce.java:90)
	at org.orekit.forces.ForceModel.addContribution(ForceModel.java:120)
	at org.orekit.propagation.numerical.NumericalPropagator$Main.computeDerivatives(NumericalPropagator.java:979)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:810)
	at org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:134)
	at org.hipparchus.ode.nonstiff.ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(ExplicitRungeKuttaIntegrator.java:131)
	at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:241)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:509)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:440)
	at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:157)
	at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(AbstractBatchLSModel.java:299)
	at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
	at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:645)
	at org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer.optimize(LevenbergMarquardtOptimizer.java:448)
	at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:459)

This happens after 23 iterations as well. I also tried turning both the atmosphere and solar forces off and it didn’t converge.

I have seen those measurement builders before. What I am doing though is taking someone’s code from 2018 and translating it from Java to Python. I’m not experienced in Orekit enough yet to know where would have been good to deviate from the old code unless the new implementations in 12.2 have forced me to.

MaximeJ · November 5, 2024, 8:18am

Hi Roger,

The new exception is again typical of a reentry. Your satellite is crashing down on Earth for some reason.
We added this exception a few months back but now I realize it isn’t explicit enough. We would at least need to output the current altitude to make it more understandable.
Could you please open a bug report on the forge regarding that point?

Is there a way you could share with us a runnable piece of code reproducing the error?

Do the same simulation work in Java?

Cheers,
Maxime

rhooper · November 5, 2024, 1:38pm

Yes, the original script comes from a professor I’m working with who wrote two papers using the data he produced from the Java version (both published in 2018).

Here’s SimulateEstimation, the class that I’m currently trying to run and getting the error:

# Original Author: Dr. Hao Peng, Embry-Riddle Aeronautical University
# Translation from Java to Python by: Roger Hooper, EV42, MSFC, NASA
# Most comments originally in Chinese, translated to English using Google Translate

from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.optim.nonlinear.vector.leastsquares import LeastSquaresOptimizer,LevenbergMarquardtOptimizer
from org.hipparchus.util import FastMath
from org.orekit.estimation.leastsquares import BatchLSEstimator
from org.orekit.estimation.measurements import ObservedMeasurement,AngularAzEl,ObservableSatellite,Range
from org.orekit.orbits import CartesianOrbit,Orbit
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.conversion import NumericalPropagatorBuilder
from org.orekit.time import AbsoluteDate
from org.orekit.utils import PVCoordinates,ParameterDriver,TimeStampedPVCoordinates
from Common import Common
from ResidentSpaceObject import ResidentSpaceObject
from typing import List
from ObserveStation import ObserveStation
from TimeUsed import TimeUsed
import random
import pandas as pd
import multiprocessing as mp

class SimulateEstimation:
    """
    Perform orbit determination simulation on observation data files
    The format of observation files is given by SimulateObservation, which is usually used after it
    """
    _OAF = Common.outputAccuracyFormat # Unified control output accuracy
    _observeFileName = "" # Observation data file
    _estimationFileName = ""
    _estimationLogFileName = ""
    _estimateDragCoefficient = True
    _estimateReflectionCoefficient = True
    _ESTIMATION_TOLERANCE = 1.0e-3
    _MAXITERATIONS = 30
    _MAXEVALUATIONS = 10*_MAXITERATIONS
    _minTrackMeasurement = 100 # The minimum number of observations required per segment when performing estimation
    _estimationInterval = 0.0 # [s] How long the track is from the current point
    _tracks: List['Track'] = []

    _RSO: ResidentSpaceObject # Simulated Resident Space Object (RSO)

    _atmosphereModel = "" # Atmosphere model for extimation and prediction procedures

    def run(self, stations: List[ObserveStation], RSO: ResidentSpaceObject, estimationInterval: float,
            observeFileName: str, estimationFileName: str, coreNumber: int, sigmaErrorOn: float,
            estimateDragCoefficient: bool, estimateReflectionCoefficient: bool, atmosphereModel: str):
        """
        Simulate the process of orbit determination of a given RSO using measurments from given stations, divided 
        into different tracks
        """
        self._RSO = RSO
        self._estimationInterval = estimationInterval
        self._observeFileName = observeFileName
        self._estimationFileName = estimationFileName
        self._estimationLogFileName = estimationFileName[0:-4] + "_Log.txt"
        self._estimateDragCoefficient = estimateDragCoefficient
        self._estimateReflectionCoefficient = estimateReflectionCoefficient
        self._atmosphereModel = atmosphereModel
        self.observableSatellite = ObservableSatellite(0)

        # If only one observation station is used, modify the minimum number of observations
        if estimationInterval == 0:
            self._minTrackMeasurement = 30
            print("Only one track is being used for each estimation!")
            print("SimulateEstimation._minTrackEstimation is set to 30!")
            print("This mode should not usually be used!")

        # Timing
        TimeUsed.timeIt("Start Estimation...")

        # Set the data file location
        # SetUpOrekitDataFile.set();

        # Read all observation data
        # Add simulated Gauss error to each observation
        # self._tracks = [self.Track()]
        if Common.testingRandomSeed != 0:
            random.seed(Common.testingRandomSeed)

        # Read all observations into tracks
        dfObserve = pd.read_csv(self._observeFileName)
        cols = ["x_true[m]","y_true[m]","z_true[m]","vx_true[m]","vy_true[m]","vz_true[m]",
                "x_est[m]","y_est[m]","z_est[m]","vx_est[m]","vy_est[m]","vz_est[m]",
                "x_err[m]","y_err[m]","z_err[m]","vx_err[m]","vy_err[m]","vz_err[m]",
                "observation_at_max_elevation","max_elevation[deg]","observation_after_max_elevation",
                "date","track_index",
                "drag_coefficient_RSO","drag_coefficient_estimated",
                "reflection_coefficient_RSO","reflection_coefficient_estimated"]
        dfEstimate = pd.DataFrame(columns=cols)
        indexMeasurement = dfObserve.iloc[0]["track"]
        self._tracks.append(self.Track())
        for index, row in dfObserve.iterrows():
            if indexMeasurement != row["track"]:
                self._tracks.append(self.Track())
                indexMeasurement = row["track"]
            
            # Read PV to provide an initial guess for subsequent calculations
            self._tracks[indexMeasurement].timeStampedTruthPVs.append(TimeStampedPVCoordinates(
                AbsoluteDate(row["date"],Common.utc()),
                Vector3D(row["x[m]"],row["y[m]"],row["z[m]"]),
                Vector3D(row["vx[m/s]"],row["vy[m/s]"],row["vz[m/s]"])
            ))

            # Read the station corresponding to the current observation
            tempInt = row["stationindex"]
            self._tracks[indexMeasurement].stationIDs.append(tempInt)
            station = stations[tempInt]

            # Read azimuth and elevation, then add error
            self._tracks[indexMeasurement].measurements.append(AngularAzEl(
                station._groundStationForEstimation,
                AbsoluteDate(row["date"],Common.utc()),
                [row["azimuth[rad]"]+random.gauss(0.0,1.0)*station._azSigma*sigmaErrorOn,row["elevation[rad]"]+random.gauss(0.0,1.0)*station._elSigma*sigmaErrorOn],
                [station._azSigma,station._elSigma],
                [1.0,1.0],
                self.observableSatellite ######## I'M NOT SURE ABOUT THIS. 0 IS THE PROPAGATOR INDEX. THERE IS ONLY ONE SATELLITE, THUS ONLY ONE PROPAGATOR ########
            ))

            # Read range and add error
            two_way = False
            self._tracks[indexMeasurement].measurements.append(Range(
                station._groundStationForEstimation,
                two_way,
                AbsoluteDate(row["date"],Common.utc()),
                row["range[m]"]*random.gauss(0.0,1.0)*station._rangeSigma*sigmaErrorOn,
                station._rangeSigma,
                1.0,
                self.observableSatellite ######## I'M NOT SURE ABOUT THIS. 0 IS THE PROPAGATOR INDEX. THERE IS ONLY ONE SATELLITE, THUS ONLY ONE PROPAGATOR ########
            ))

        print(f"{len(self._tracks)} tracks in total.")
        for index,track in enumerate(self._tracks):
            print(f"Track {index} has {len(track.timeStampedTruthPVs)} observations and {len(track.measurements)} measurements.")
        TimeUsed.timeIt("Finished reading observations and adding noise.")

        # Orbit determination simulation
        if coreNumber > 1:
            dfTemp = pd.DataFrame(columns=cols)
            grid = [(track,dfTemp) for track in self._tracks]
            pool = mp.Pool(coreNumber)
            calc = pool.starmap(self.estimationOneState,grid)
            pool.close()
            pool.join()
            dfEstimate = dfEstimate.append(calc)
            dfEstimate = dfEstimate.reset_index(drop=True)
        else:
            for track in self._tracks:
                dfEstimate = self.estimationOneState(track,dfEstimate)
    
    def estimationOneState(self, track: 'Track', dfEstimate: pd.DataFrame):
        RSO = self._RSO
        initialOrbitGuess = CartesianOrbit(track.timeStampedTruthPVs[0],Common.inertialFrame,Common.muEarth) # Take the first time of each track as the moment of observation
        initialStateGuess = SpacecraftState(initialOrbitGuess,self._RSO.mass)

        # Use the same dynamics model, but with different initial dragCoefficient and reflectionCoeffient
        estimationPropagatorBuilder = Common.propagatorBuilder(
            initialStateGuess,RSO.area,RSO.dragCoefficientEstimationAndPrediction,RSO.reflectionCoefficientEstimationAndPrediction,self._atmosphereModel
        )

        # Construct Batch Least Squares Estimator
        optimizer = LevenbergMarquardtOptimizer()
        estimator = BatchLSEstimator(optimizer,[estimationPropagatorBuilder])
        estimator.setParametersConvergenceThreshold(self._ESTIMATION_TOLERANCE)
        estimator.setMaxIterations(self._MAXITERATIONS)
        estimator.setMaxEvaluations(self._MAXEVALUATIONS)

        # Add other supported estimated variables
        driver: ParameterDriver
        for driver in estimator.getPropagatorParametersDrivers(False).getDrivers():
            if driver.getName() == "drag coefficient":
                driver.setSelected(self._estimateDragCoefficient)
            if driver.getName() == "reflection coefficient":
                driver.setSelected(self._estimateReflectionCoefficient)
            if driver.getName() == "cross section":
                driver.setSelected(False)

        # Add the measurement observations (including the error) within the previous estimationInterval of the current track to the estimator.
        measurementNumber: int = 0
        if self._estimationInterval <= 0:
            previousTracks: List['SimulateEstimation.Track'] = self._tracks[0:self._tracks.index(track)+1]

            # When the length of estimationInterval is reached, the calculation will start again
            if previousTracks[0].timeStampedTruthPVs[0].getDate().compareTo(
                track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)) <= 0:
                oneTrack: SimulateEstimation.Track
                for oneTrack in previousTracks:
                    if oneTrack.timeStampedTruthPVs[len(oneTrack.timeStampedTruthPVs)-1].getDate().compareTo(track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)) < 0:
                        continue
                    oneMeasurement: ObservedMeasurement
                    for oneMeasurement in oneTrack.measurements:
                        estimator.addMeasurement(oneMeasurement)
                        measurementNumber += 1
        elif self._estimationInterval > 0:
            laterTracks: List[SimulateEstimation.Track] = self._tracks[self._tracks.index(track):len(self._tracks)]
            
            # When the length of estimationInterval is reached, the calculation will start again
            if laterTracks[-1].timeStampedTruthPVs[len(laterTracks[-1].timeStampedTruthPVs)-1].getDate().compareTo(
                track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)
            ) >= 0:
                oneTrack: SimulateEstimation.Track
                for oneTrack in laterTracks:
                    if oneTrack.timeStampedTruthPVs[0].getDate().compareTo(track.timeStampedTruthPVs[0].getDate().shiftedBy(self._estimationInterval)) > 0:
                        break
                    oneMeasurement: ObservedMeasurement
                    for oneMeasurement in oneTrack.measurements:
                        estimator.addMeasurement(oneMeasurement)
                        measurementNumber += 1

        # If there are more than the maximum number of observations, then integration begins
        if measurementNumber >= self._minTrackMeasurement:
            tempLogString = ""

            # Run an estimate
            ######## Orekit 12.0 returns a Propagator list for estimate(), so I'm guessing tat grabbing [0] is correct ########
            estimated: Orbit = estimator.estimate()[0].getInitialState().getOrbit()

            # Output data - pt = position truth, pe = position estimation, dp = difference of position 
            pt: Vector3D
            vt: Vector3D
            pe: Vector3D
            ve: Vector3D
            dp: Vector3D
            dv: Vector3D

            tempEstimateData = []
            
            # Set truth values
            pt = track.timeStampedTruthPVs[0].getPosition()
            vt = track.timeStampedTruthPVs[0].getVelocity()
            tempEstimateData = tempEstimateData + [pt.getX(),pt.getY(),pt.getZ()]
            tempEstimateData = tempEstimateData + [vt.getX(),vt.getY(),vt.getZ()]

            # Set estimation values
            pe = estimated.getPVCoordinates().getPosition()
            ve = estimated.getPVCoordinates().getVelocity()
            tempEstimateData = tempEstimateData + [pe.getX(),pe.getY(),pe.getZ()]
            tempEstimateData = tempEstimateData + [ve.getX(),ve.getY(),ve.getZ()]

            # Set truth minus estimation values
            dp = pt.add(-1,pe)
            dv = vt.add(-1,ve)
            tempEstimateData = tempEstimateData + [dp.getX(),dp.getY(),dp.getZ()]
            tempEstimateData = tempEstimateData + [dv.getX(),dv.getY(),dv.getZ()]

            # The azimuth corresponding to the maximum elevation moment on the current track, elevation, range [rad, m]
            maxElevation = 0
            maxElevationIndex = 0
            m: ObservedMeasurement
            for m in track.measurements:
                if len(m.getObservedValue()) == 1:
                    continue
                if maxElevation < m.getObservedValue()[1]:
                    maxElevation = m.getObservedValue()[1]
                    maxElevationIndex = track.measurements.index(m)

            tempEstimateData = tempEstimateData + [track.measurements[maxElevationIndex].getObservedValue()[0],
                                                    maxElevation,
                                                    track.measurements[maxElevationIndex+1].getObservedValue()[0]]
            
            # Time + index
            tempEstimateData = tempEstimateData + [f"{estimated.getDate()}",f"{self._tracks.index(track)}"]

            # Cd_true + Cd_estimated + Cr_true + Cr_estimated
            pd: ParameterDriver
            for pd in estimator.getPropagatorParametersDrivers(False).getDrivers():
                if pd.getName() == "drag coefficient":
                    tempEstimateData = tempEstimateData + [RSO.dragCoefficient,pd.getValue()]
                if pd.getName() == "reflection coefficient":
                    tempEstimateData = tempEstimateData + [RSO.reflectionCoefficient,pd.getValue()]
            
            # Add row to output DataFrame
            dfEstimate.loc[len(dfEstimate)] = tempEstimateData

            # Output Screen
            tempLogString = tempLogString + f"Result on track {self._tracks.index(track)}:\n"
            tempLogString = tempLogString + f"{measurementNumber} measurements on track {self._tracks.index(track)}.\n"
            tempLogString = tempLogString + f"Truth:     \n"
            tempLogString = tempLogString + f"{track.timeStampedTruthPVs[0]}\n"
            tempLogString = tempLogString + f"drag coefficient = {RSO.dragCoefficient}\n"
            tempLogString = tempLogString + f"Guess:     \n"
            tempLogString = tempLogString + f"{initialOrbitGuess.getPVCoordinates()}\n"
            tempLogString = tempLogString + f"drag coefficient = {RSO.dragCoefficientEstimationAndPrediction}\n"
            tempLogString = tempLogString + f"Estimated: \n"
            tempLogString = tempLogString + f"{estimated.getPVCoordinates()}\n"
            pd: ParameterDriver
            for pd in estimator.getPropagatorParametersDrivers(True).getDrivers():
                tempLogString = tempLogString + f"{pd.getName()} = {pd.getValue()}\n"
            tempLogString = tempLogString + f"\n{estimator.getIterationsCount()} iterations, and {estimator.getEvaluationsCount()} evaluations.\n"
            
            # Output log: Output the same results as on the screen to a log file
            tempLogString = tempLogString + f"{estimated.getDate().durationFrom(RSO.initialDate)} | {dp.getNorm()} | {dp.getNorm()/pt.getNorm()*100}% | {dv.getNorm()} | {dv.getNorm()/vt.getNorm()*100}%.\n"
            print(tempLogString)
            with open(self._estimationLogFileName,'a') as f:
                f.write(tempLogString)
        else:
            tempOut = f"# Track {self._tracks.index(track)} is not estimated, because there are not enough measurements on this track."
            print(tempOut)

        return dfEstimate


    class Track:
        """
        Stores the observation data of each track, including the PV of each point and the corresponding 
        measurements
        """
        def __init__(self):
            self.timeStampedTruthPVs: List[TimeStampedPVCoordinates] = []
            self.stationIDs: List[int] = []
            self.measurements: List[ObservedMeasurement] = []

You can comment out the lines that use “TimeUsed.” It’s not important to getting it to run.

Here’s the “Common” class:

# Original Author: Dr. Hao Peng, Embry-Riddle Aeronautical University
# Translation from Java to Python by: Roger Hooper, EV42, MSFC, NASA
# Most comments originally in Chinese, translated to English using Google Translate

# The same numerical integrator was used throughout the study (excluding the mechanical model)

import SetUpOrekitDataFile
from org.orekit.bodies import CelestialBodyFactory,OneAxisEllipsoid
from org.orekit.data import DataProvidersManager,DataContext
from org.orekit.errors import OrekitException
from org.orekit.forces import ForceModel
from org.orekit.forces.drag import DragForce,IsotropicDrag
from org.orekit.models.earth.atmosphere import Atmosphere,DTM2000,NRLMSISE00,SimpleExponentialAtmosphere,PythonNRLMSISE00InputParameters
from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation,MarshallSolarActivityFutureEstimationLoader
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel,ThirdBodyAttraction
from org.orekit.forces.gravity.potential import GravityFieldFactory,NormalizedSphericalHarmonicsProvider
from org.orekit.forces.radiation import IsotropicRadiationSingleCoefficient,SolarRadiationPressure
from org.orekit.frames import Frame,FramesFactory,FactoryManagedFrame
from org.orekit.orbits import CartesianOrbit,Orbit,OrbitType,PositionAngleType
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.conversion import DormandPrince853IntegratorBuilder,NumericalPropagatorBuilder,ODEIntegratorBuilder
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.time import AbsoluteDate,TimeScale,TimeScalesFactory
from org.orekit.utils import Constants,IERSConventions,PVCoordinates,PVCoordinatesProvider
from typing import List

class Common():

    minStep = 0.001 # [s]
    maxStep = 90.0 # [s] After experimentation, this is not too slow and can capture most situations
    positionTolerance = 0.1 # [m]
    inertialFrame: FactoryManagedFrame = FramesFactory.getEME2000() # J2000
    muEarth = Constants.EIGEN5C_EARTH_MU # The Earth's gravitational constant (to be consistent with the mechanical model used later, the difference between the models is very small)
    radiusEarth = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
    flatteningEarth = Constants.WGS84_EARTH_FLATTENING
    areaMin = 0.001 # [m^2]
    areaMax = 100000
    outputAccuracyFormat = "+20.15E"
    testingRandomSeed = 19881226

    # Integrator Parameters
    # def __init__(self):
        # Integrator parameters
        # self.minStep = 0.001 # [s]
        # self.maxStep = 90.0 # [s] After experimentation, this is not too slow and can capture most situations
        # self.positionTolerance = 0.1 # [m]

        # Common model parameters
        # self.inertialFrame = FramesFactory.getEME2000() # J2000
        # self.muEarth = Constants.EIGEN5C_EARTH_MU # The Earth's gravitational constant (to be consistent with the mechanical model used later, the difference between the models is very small)
        # self.radiusEarth = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
        # self.flatteningEarth = Constants.WGS84_EARTH_FLATTENING

        # Spacecraft parameters
        # self.areaMin = 0.001 # [m^2]
        # self.areaMax = 100000 # [m^2]

        # Global output Accuracy
        # self.outputAccuracyFormat = "%+20.15E"

        # In order to ensure the repeatability of the calculation, the random number seed used
        # self.testingRandomSeed = 19881226

    # Earth shape model
    @staticmethod
    def earthShape() -> OneAxisEllipsoid:
        # SetUpOrekitDataFile.set()
        # shape = OneAxisEllipsoid(self.radiusEarth,self.flatteningEarth,FramesFactory.getITRF(IERSConventions.IERS_2010,True))
        shape = OneAxisEllipsoid(Common.radiusEarth,Common.flatteningEarth,FramesFactory.getITRF(IERSConventions.IERS_2010,True))
        return shape
    
    # Establish a common integrator
    @staticmethod
    def integrator() -> ODEIntegratorBuilder:
        builder = DormandPrince853IntegratorBuilder(Common.minStep,Common.maxStep,Common.positionTolerance)
        return builder
    
    # Establish a common propagator
    @staticmethod
    def propagatorForObservation(spacecraftState: SpacecraftState, spacecraftArea, spacecraftDragCoefficient, spacecraftReflectionCoefficient) -> NumericalPropagator:
        np = NumericalPropagator(Common.integrator().buildIntegrator(CartesianOrbit(spacecraftState.getOrbit()),OrbitType.CARTESIAN))
        np.setOrbitType(OrbitType.CARTESIAN)
        np.setPositionAngleType(PositionAngleType.TRUE) # The position angle used in the integration should have nothing to do with the angle in Orbit
        np.setInitialState(spacecraftState) # Set the initial state of the integrator, including quality, etc.

        # Add mechanical models required for simulation
        np.addForceModel(Common._harmonicForce(40,40))
        for iiForce in Common._thirdBodyForces(10):
            np.addForceModel(iiForce)
        np.addForceModel(Common._myAtmosphereDragForce(spacecraftArea,spacecraftDragCoefficient,"DTM2000"))
        np.addForceModel(Common._solarRadiationPressureForce(spacecraftArea,spacecraftReflectionCoefficient))
        return np

    # Integrator used for estimation and prediction, because different mechanical models can be set here to introduce model errors
    @staticmethod
    def propagatorForEstimationAndPrediction(spacecraftState: SpacecraftState, spacecraftArea: float, spacecraftDragCoefficient: float, spacecraftReflectionCoefficient: float, atmosphereModel: str) -> NumericalPropagator:
        np = NumericalPropagator(
            Common.integrator().buildIntegrator(CartesianOrbit(spacecraftState.getOrbit()),OrbitType.CARTESIAN))
        np.setOrbitType(OrbitType.CARTESIAN)
        np.setPositionAngleType(PositionAngleType.TRUE) # The position angle used in the integration should have nothing to do with the angle in Orbit
        np.setInitialState(spacecraftState) # Set the initial state of the integrator, including quality, etc.
        # Add the mechanical model required for simulation <------------------------------------------------------- The mechanical model used is set to be different from that used in simulation
        np.addForceModel(Common._harmonicForce(10,10)) # <------------------------------------- The mechanical model used is set to be different from that used in the simulation
        for iiForce in Common._thirdBodyForces(3): # <------------------------------------- The mechanical model used is set to be different from that used in the simulation
            np.addForceModel(iiForce)
        np.addForceModel(Common._myAtmosphereDragForce(spacecraftArea,spacecraftDragCoefficient,atmosphereModel))
        np.addForceModel(Common._solarRadiationPressureForce(spacecraftArea,spacecraftReflectionCoefficient))
        return np
    
    @staticmethod
    def propagatorBuilder(spacecraftState: SpacecraftState, spacecraftArea: float, spacecraftDragCoefficient: float, spacecraftReflectionCoefficient: float, atmosphereModel: str) -> NumericalPropagatorBuilder:
        """
        * Specifies the builder of the integrator of the track. It is needed to construct the estimator during estimation.
        *
        * Usually, do not use this builder when constructing an integrator.
        * I don't quite understand the significance of setting it in Orekit. I have only seen it used in the estimator.
        """
        # Use propagator() to ensure that the constructed integrators are exactly the same
        np = Common.propagatorForEstimationAndPrediction(spacecraftState,spacecraftArea,spacecraftDragCoefficient,spacecraftReflectionCoefficient,atmosphereModel)
        # Constructs the same initial track type as the propagator() default
        orbitWithConsistentType = np.getOrbitType().convertType(spacecraftState.getOrbit())
        # Builder
        npBuilder = NumericalPropagatorBuilder(orbitWithConsistentType,Common.integrator(),np.getPositionAngleType(),Common.positionTolerance)
        # Add spacecraft mass
        npBuilder.setMass(spacecraftState.getMass())
        # Add mechanical models
        for iiForce in np.getAllForceModels():
            npBuilder.addForceModel(iiForce)
        
        return npBuilder

    @staticmethod
    def _harmonicForce(degree: int, order: int) -> ForceModel:
        harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree,order) # Spherical harmonics
        attractionModel = HolmesFeatherstoneAttractionModel(Common.earthShape().getBodyFrame(),harmonicsGravityProvider)
        return attractionModel

    @staticmethod
    def _thirdBodyForces(number: int) -> List[ForceModel]:
        tempForces = []
        tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getSun()))
        tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getMoon()))
        tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getJupiter()))
        if number == 10:
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getMercury()))
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getVenus()))
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getMars()))
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getSaturn()))
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getNeptune()))
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getUranus()))
            tempForces.append(ThirdBodyAttraction(CelestialBodyFactory.getPluto()))

        return tempForces
    
    @classmethod
    def _myAtmosphereDragForce(self, spacecraftArea: float, spacecraftDragCoefficient: float, atmosphereModel: str) -> DragForce:
        # Create desired atmosphere model
        myAtmosphereModel = 0
        case = atmosphereModel.upper()
        if case == "DTM2000":
            msafe = MarshallSolarActivityFutureEstimation("\\p{Alpha}\\p{Lower}\\p{Lower}\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}(?:f|F)10(?:[-_]prd)?\\.(?:txt|TXT)",
                                                          MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE) # Provide solar activity forecast data
            DM = DataContext.getDefault().getDataProvidersManager()
            ################### Changing from msafe to msafeLoader #################
            ########################################################################
            msafeLoader = MarshallSolarActivityFutureEstimationLoader(MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
            DM.feed(msafe.getSupportedNames(),msafeLoader) # Load data
            ########################################################################
            myAtmosphereModel = DTM2000(msafe, CelestialBodyFactory.getSun(), self.earthShape()) # Set DTM2000 as atmospheric model

        elif case == "NRLMSISE00" or case == "NRLMSISE-00":
            # Atmosphere model NRLMSISE-00
            myAtmosphereModel = NRLMSISE00(SimpleParametersForNRLMSISE00(),CelestialBodyFactory.getSun(),self.earthShape())

        elif case == "EXPONENTIAL":
            # Exponential atmosphere model
            # Parameters selected from: http://scipp.ucsc.edu/outreach/balloon/glost/environment3.html
            # (Not sure if this is correct)
            myAtmosphereModel = SimpleExponentialAtmosphere(self.earthShape(),1.225e3,0,8.42e3)
            print("CHECK THE SIMPLE EXPONENTIAL ATMOSPHERE MODEL!")
        else:
            print(f"\n!!! Invalid or unsupported atmosphere model: {atmosphereModel}")
            raise(OrekitException())
        
        spacecraftShapeForDrag = IsotropicDrag(spacecraftArea,spacecraftDragCoefficient)
        dragForce = DragForce(myAtmosphereModel,spacecraftShapeForDrag)
        return dragForce
    
    def _solarRadiationPressureForce(spacecraftArea: float, spacecraftReflectionCoefficient: float) -> ForceModel:
        sun = PVCoordinatesProvider.cast_(CelestialBodyFactory.getSun())
        solarRadiationPressure = SolarRadiationPressure(sun,
                                                    Common.earthShape(),
                                                    IsotropicRadiationSingleCoefficient(spacecraftArea,spacecraftReflectionCoefficient))
        return solarRadiationPressure

    def utc() -> TimeScale:
        """
        Common time scale, UTC, for all simulations
        """
        # SetUpOrekitDataFile.set()
        return TimeScalesFactory.getUTC()

class SimpleParametersForNRLMSISE00(PythonNRLMSISE00InputParameters):
    def __init__(self):
        super(SimpleParametersForNRLMSISE00,self).__init__()

    def finalize(self):
        pass

    def getAp(self,date: AbsoluteDate):# -> list[float]:
        return [4.0,100.0,100.0,100.0,100.0,100.0,100.0]
    
    def getAverageFlux(self,date: AbsoluteDate) -> float:
        return 150.0
    
    def getDailyFlux(self,date: AbsoluteDate) -> float:
        return 150.0

    def getMaxDate(self) -> AbsoluteDate:
        absDate = AbsoluteDate(2018,12,31,TimeScalesFactory.getUTC())
        return absDate
    
    def getMinDate(self) -> AbsoluteDate:
        absDate = AbsoluteDate(2017,1,1,TimeScalesFactory.getUTC())
        return absDate
    
    def pythonDecRef(self):
        return super().pythonDecRef()

rhooper · November 5, 2024, 1:39pm

Here’s the ResidentSpaceObject class:

# Original Author: Dr. Hao Peng, Embry-Riddle Aeronautical University
# Translation from Java to Python by: Roger Hooper, EV42, MSFC, NASA
# Most comments originally in Chinese, translated to English using Google Translate

# Used to store the data of the resident space object required for simulation
# * mass, surface area, reflection coefficient, atmospheric drag coefficient...

# Orbit information is not included because orbit can be specified arbitrarily, and in our research, the rest of the properties of RSO are basically unchanged

from org.hipparchus.util import FastMath
from org.orekit.orbits import Orbit,CartesianOrbit,KeplerianOrbit,PositionAngleType
from org.orekit.propagation.analytical.tle import TLE,TLEPropagator
from org.orekit.time import AbsoluteDate
from org.orekit.utils import PVCoordinates,PVCoordinatesProvider
from Common import Common

class ResidentSpaceObject():
    def __init__(self,name: str, mass: float=0, area: float=0, dragCoefficient: float=2.2, reflectionCoefficient: float=1.25,
                 initialDate: AbsoluteDate=None, initialOrbit: CartesianOrbit=None):
        # RSO Attribute Parameters
        self.name = name
        self.NORADID = 0
        self.mass = mass
        self.area = area
        self.dragCoefficient = dragCoefficient
        self.dragCoefficientEstimationAndPrediction = 2.2
        self.reflectionCoefficient = reflectionCoefficient
        self.reflectionCoefficientEstimationAndPrediction = 1.25

        # Orbital Parameters
        if initialDate:
            self.initialDate = initialDate
        if initialOrbit:
            self.initialOrbit = initialOrbit
            self.a = initialOrbit.getA()
            self.e = initialOrbit.getE()
            self.i = initialOrbit.getI()
            self.omega = KeplerianOrbit(initialOrbit).getPerigeeArgument()
            self.raan = KeplerianOrbit(initialOrbit).getRightAscensionOfAscendingNode()
            self.positionAngle = PositionAngleType.TRUE
            self.angle = KeplerianOrbit(initialOrbit).getAnomaly(self.positionAngle)

        # Setting up the Orekit database
        # SetUpOrekitDataFile.set()

        if name.upper() == "ENVISAT":
            # TRUE parameters
			# 2017/08/27 13:05:00 extracted from leo-spctrk.txt of PreviSat 3.5 on MacOS 
			# 0 ENVISAT
			# 1 27386U 02009A   17238.87724141 +.00000000 +00000-0 +13804-4 0  9992
			# 2 27386 098.2254 281.1813 0001340 081.3267 278.8065 14.37901343811233
            self.tleEnvisat = TLE("1 27386U 02009A   17238.87724141 +.00000000 +00000-0 +13804-4 0  9992",
                             "2 27386 098.2254 281.1813 0001340 081.3267 278.8065 14.37901343811233")
            self.initialDate = self.tleEnvisat.getDate()
            extrapolator = TLEPropagator.selectExtrapolator(self.tleEnvisat)
            pvProvider: PVCoordinatesProvider = PVCoordinatesProvider.cast_(extrapolator) # Cast the extrapolator to PVCoordinatesProvider to gain access to getPVCoordinates(AbsoluteDate,Frame) overload
            self.initialOrbit = CartesianOrbit(pvProvider.getPVCoordinates(self.initialDate,Common.inertialFrame),
                                               Common.inertialFrame,self.initialDate,Common.muEarth) # This has been fixed: TEMEE to EME2000 conversion
            self.mass = 8211.0 # [kg]
            self.area = 16.0 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            self.dragCoefficient = 2.2
            self.dragCoefficientEstimationAndPrediction = 2.2 # The drag coefficient used during observation
            self.reflectionCoefficient = 1.25 # A generic reflection coefficient given by Vallado.

            self.a = self.initialOrbit.getA()
            self.e = self.initialOrbit.getE()
            self.i = self.initialOrbit.getI()
            self.omega = KeplerianOrbit(self.initialOrbit).getPerigeeArgument()
            self.raan = KeplerianOrbit(self.initialOrbit).getRightAscensionOfAscendingNode()
            self.positionAngle = PositionAngleType.TRUE
            self.angle = KeplerianOrbit(self.initialOrbit).getAnomaly(self.positionAngle)

        # All data below taken from https://heavens-above.com/orbit.aspx?satid=xxxxx
        else:
            if name == "paper31-ISS":
                self.tmpTLE = TLE("1 25544U 98067A   18247.55883065  .00016717  00000-0  10270-3 0  9043",
                             "2 25544  51.6419 343.1595 0005693 129.7955 230.3700 15.53945776 10850")
                self.mass = 419725.0 # [kg]
                self.area = 0.2*(72.8*108.5) # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-LARETS":
                self.tmpTLE = TLE("1 27944U 03042F   18246.50353587 -.00000051  00000-0 -81033-6 0  9991",
                                  "2 27944  98.0375  15.2582 0011429 310.4899  49.5307 14.63172440797557")
                self.mass = 23.28 # [kg]
                self.area = FastMath.PI*0.12*0.12 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-STELLA":
                self.tmpTLE = TLE("1 22824U 93061B   18246.18488426 -.00000034  00000-0  46876-5 0  9993",
                                  "2 22824  98.8801 213.4116 0006099 190.9595   0.5798 14.27374063299202")
                self.mass = 48.0 # [kg]
                self.area = FastMath.PI*0.12*0.12 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-BLITS":
                self.tmpTLE = TLE("1 35871U 09049G   18244.20890810 -.00000034 +00000-0 +31883-5 0  9999",
                                  "2 35871 098.4104 246.8687 0002524 317.6258 042.4725 14.21785510464033")
                self.mass = 8211.0 # [kg]
                self.area = 16.0 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-HAIYANG-2A":
                self.tmpTLE = TLE("1 37781U 11043A   18246.55461026  .00000375  00000-0  41240-3 0  9994", 
                                  "2 37781  99.3243 255.4140 0000608  55.7211 304.4041 13.78717189355072")
                self.mass = 1500.0 # [kg]
                self.area = 38.94 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-SWARM-A":
                self.tmpTLE = TLE("1 39452U 13067B   18246.18535172  .00000921  00000-0  19667-4 0  9998",
                                  "2 39452  87.3506   8.9093 0003256  83.5516 276.6111 15.43304605268435")
                self.mass = 369.0 # [kg]
                self.area = 1.275 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-STARLETTE":
                self.tmpTLE = TLE("1 07646U 75010A   18246.47129393 -.00000140  00000-0  36309-5 0  9997",
                                  "2 07646  49.8233 168.7493 0205743  29.0956 332.1221 13.82305008201168")
                self.mass = 47.5 # [kg]
                self.area = FastMath.PI*0.12*0.12 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-BEACON-C":
                self.tmpTLE = TLE("1 01328U 65032A   18245.54589041 -.00000086  00000-0 -11973-4 0  9997",
                                  "2 01328  41.1869 305.8765 0254535 269.3200  87.8385 13.38575203607532")
                self.mass = 60.0 # [kg]
                self.area = 2.2651 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-LARES":
                self.tmpTLE = TLE("1 38077U 12006A   18246.11501872  .00000018  00000-0  35121-3 0  9994",
                                  "2 38077  69.4934 110.1713 0010914 217.0386 142.9890 12.54931564300341")
                self.mass = 400.0 # [kg]
                self.area = FastMath.PI*0.18*0.18 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-AJISAT":
                self.tmpTLE = TLE("1 16908U 86061A   18246.53469451 -.00000079  00000-0  12839-3 0  9996",
                                  "2 16908  50.0070 251.3190 0011379 142.4877  36.6732 12.44488521126648")
                self.mass = 685.0 # [kg]
                self.area = FastMath.PI*1.075*1.075 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-LAGEOS-1":
                self.tmpTLE = TLE("1 08820U 76039A   18246.23521883  .00000024  00000-0  00000+0 0  9991",
                                  "2 08820 109.8466 285.1654 0044696 191.1608 198.5007  6.38664811731936")
                self.mass = 406.965 # [kg]
                self.area = FastMath.PI*0.3*0.3 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-LAGEOS-2":
                self.tmpTLE = TLE("1 22195U 92070B   18245.69377977 -.00000009  00000-0  00000+0 0  9991",
                                  "2 22195  52.6459 264.3746 0138125  23.9743 179.0148  6.47294245611474")
                self.mass = 405.38 # [kg]
                self.area = FastMath.PI*0.3*0.3 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-ETALON-1":
                self.tmpTLE = TLE("1 19751U 89001C   18246.15301736  .00000010  00000-0  00000+0 0  9998",
                                  "2 19751  64.1481 169.8510 0021847 209.7562  32.3517  2.13156267230810")
                self.mass = 1415.0 # [kg]
                self.area = FastMath.PI*0.647*0.647 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            elif name == "paper31-ETALON-2":
                self.tmpTLE = TLE("1 20026U 89039C   18244.73162994 +.00000004 +00000-0 +00000-0 0  9993",
                                  "2 20026 065.2522 047.4358 0017224 223.1424 046.6927 02.13204417227839")
                self.mass = 1415.0 # [kg]
                self.area = FastMath.PI*0.647*0.647 # [m^2] Calculate the cross-sectional area for atmospheric drag and light pressure
            else:
                raise(Exception(f"No such RSO: {name}"))
            self.NORADID = self.tmpTLE.getSatelliteNumber()
            self.initialDate = self.tmpTLE.getDate()
            self.initialOrbit = CartesianOrbit(TLEPropagator.selectExtrapolator(self.tmpTLE).getPVCoordinates(self.initialDate,Common.inertialFrame),
                                               Common.inertialFrame,self.initialDate,Common.muEarth)
            self.a = initialOrbit.getA()
            self.e = initialOrbit.getE()
            self.i = initialOrbit.getI()
            self.omega = KeplerianOrbit(self.initialOrbit).getPerigeeArgument()
            self.raan = KeplerianOrbit(self.initialOrbit).getRightAscensionOfAscendingNode()
            self.positionAngle = PositionAngleType.TRUE
            self.angle = KeplerianOrbit(self.initialOrbit).getAnomaly(self.positionAngle)

    def newChangeA(self,dA: float):
        """
        Change semi-major axis of RSO
        param: dA - [m]
        """
        newRSO = ResidentSpaceObject(self.name+"_dA"+f"{dA/1e3:.0f}", # New name
                                     self.mass,self.area,self.dragCoefficient,self.reflectionCoefficient,
                                     self.initialDate,
                                     CartesianOrbit(
                                         KeplerianOrbit(self.a + dA, # new value here
                                                        self.e,
                                                        FastMath.acos(-1.99096871e-7*FastMath.sqrt(398600.440e9)/1.5/1.7555e25*FastMath.pow(self.a + dA,7.0/2.0)),
                                                        self.omega,self.raan,
                                                        self.angle,self.positionAngle,
                                                        Common.inertialFrame,
                                                        self.initialDate,Common.muEarth) # Temporary Keplerian track
                                     ))
        return newRSO
    
    def newChangeA(self,dRAAN: float):
        """
        Change right ascension of ascending node of RSO
        param: dRAAN - [rad]
        """
        newRSO = ResidentSpaceObject(self.name+"_dRAAN"+f"{FastMath.toDegrees(dRAAN):.0f}", # New name
                                     self.mass,self.area,self.dragCoefficient,self.reflectionCoefficient,
                                     self.initialDate,
                                     CartesianOrbit(
                                         KeplerianOrbit(self.a,self.e,self.i,
                                                        self.omega,
                                                        self.raan + dRAAN, # new value here
                                                        self.angle,self.positionAngle,
                                                        Common.inertialFrame,
                                                        self.initialDate,Common.muEarth) # Temporary Keplerian track
                                     ))
        return newRSO
    
    def setMass(self, mass: float):
        self.mass = mass

    def setArea(self, area: float):
        self.area = area

    def setDragCoefficient(self, dragCoefficient: float):
        self.dragCoefficient = dragCoefficient

    def setReflectionCoefficient(self, reflectionCoefficient: float):
        self.reflectionCoefficient = reflectionCoefficient

Here’s the ObserveStation class:

# Original Author: Dr. Hao Peng, Embry-Riddle Aeronautical University, 11-15-2016
# Translation from Java to Python by: Roger Hooper, EV42, MSFC, NASA, 9-20-2024
# Several comments originally in Chinese, translated to English using Google Translate

from org.hipparchus.util import FastMath
from org.orekit.bodies import GeodeticPoint
from org.orekit.errors import OrekitException
from org.orekit.estimation.measurements import GroundStation
from org.orekit.frames import TopocentricFrame
from Common import Common

class ObserveStation:
    def __init__(self, name: str):
        """
        Use observatory data from the following literature:

        Hill, K., Sabol, C. & Alfriend, K.T., 2012. Comparison of Covariance Based Track Association Approaches Using Simulated Radar Data. 
	    The Journal of the Astronautical Sciences, 59(1C2), pp.281C300. Available at: http://link.springer.com/10.1007/s40295-013-0018-1.
        """
        common = Common() # Instantiate Common and its attributes
        self._name = name
        
        if name == "Eglin":
            self._latitude = FastMath.toRadians(30.57) # [rad]
            self._longitude = FastMath.toRadians(-86.21) # [rad]
            self._altitude = 34.7 # [m]
            self._minElevation = FastMath.toRadians(1.0) # Minimum viewing angle
            self._maxRange = 13210.0e5 # Maximum observation distance, [m]
            self._azSigma = FastMath.toRadians(0.0154) # [rad]
            self._elSigma = FastMath.toRadians(0.0147) # [rad]
            self._rangeSigma = 32.1 # [m]
        elif name == "Kaena Point":
            self._latitude = FastMath.toRadians(21.57) # [rad]
            self._longitude = FastMath.toRadians(-158.27) # [rad]
            self._altitude = 300.2 # [m]
            self._minElevation = FastMath.toRadians(0.0) # Minimum viewing angle
            self._maxRange = 6380.0e3 # Maximum observation distance, [m]
            self._azSigma = FastMath.toRadians(0.0224) # [rad]
            self._elSigma = FastMath.toRadians(0.0139) # [rad]
            self._rangeSigma = 92.5 # [m]
        elif name == "Clear":
            self._latitude = FastMath.toRadians(64.29) # [rad]
            self._longitude = FastMath.toRadians(-149.19) # [rad]
            self._altitude = 213.3 # [m]
            self._minElevation = FastMath.toRadians(1.0) # Minimum viewing angle
            self._maxRange = 4910.0e3 # Maximum observation distance, [m]
            self._azSigma = FastMath.toRadians(0.0791) # [rad]
            self._elSigma = FastMath.toRadians(0.0240) # [rad]
            self._rangeSigma = 62.5 # [m]
        else:
            print(f"No such observation station: {name}")
            raise(OrekitException())
        
        self._geodeticPoint: GeodeticPoint = GeodeticPoint(self._latitude,self._longitude,self._altitude) # Geometric points of observation stations
        self._topocentricFrame: TopocentricFrame = TopocentricFrame(common.earthShape(),self._geodeticPoint,name) # The coordinate system of the observation station (the coordinate system actually used for calculation)
        self._groundStationForEstimation: GroundStation = GroundStation(self._topocentricFrame)

        # Close the estimated parameters of the station
        drivers = [self._groundStationForEstimation.getEastOffsetDriver(),
        self._groundStationForEstimation.getNorthOffsetDriver(),
        self._groundStationForEstimation.getZenithOffsetDriver(),
        self._groundStationForEstimation.getPrimeMeridianOffsetDriver()]
        for driver in drivers:
            # driver.setReferenceDate()
            driver.setSelected(False)

rhooper · November 5, 2024, 1:49pm

If you need it, here’s the first 9 tracks of data, which is what was used to produce the error:

,x[m],y[m],z[m],vx[m/s],vy[m/s],vz[m/s],azimuth[rad],elevation[rad],range[m],date,stationindex,track
0,-902374.555616142,-639881.030163178,7048824.829466281,-1484.9948915835032,7308.192433028367,472.3922924705768,0.3659750349880249,0.0610238490542397,2881240.221548581,2017-08-26T21:27:13.657824Z,1,0
1,-989642.9308211908,-200419.6298103351,7063265.403292177,-1422.9983929917971,7335.722529205721,8.803443316195668,0.2829874810477769,0.1282484370218327,2531486.6479531303,2017-08-26T21:28:13.657824Z,1,0
2,-1073022.5468289023,239829.13413332915,7049880.587195259,-1355.4118124042247,7334.4293720317555,-454.8181694156617,0.1708986971616245,0.2027592460890204,2206076.0020055734,2017-08-26T21:29:13.657824Z,1,0
3,-1152185.8106272046,679135.4465979201,7008722.852708407,-1282.49902352265,7304.317504244945,-916.6585257873236,0.0162075878989685,0.2834180452777631,1919073.0676831985,2017-08-26T21:30:13.657824Z,1,0
4,-1226821.576321234,1115773.0966435308,6939953.342897174,-1204.5442139762029,7245.502407378679,-1374.910374074767,6.085570374747363,0.3629765687962213,1691573.6806041582,2017-08-26T21:31:13.657824Z,1,0
5,-1296636.3348092912,1548026.1123516378,6843841.276625583,-1121.8509448367734,7158.209826425679,-1827.779415481821,5.804147596542931,0.4224034999209775,1551066.9418047443,2017-08-26T21:32:13.657824Z,1,0
6,-1361355.3484390436,1974195.3695270407,6720762.989167661,-1034.7411857492457,7042.775578913478,-2273.49052675213,5.479649247712456,0.4355176870484315,1522431.1161303734,2017-08-26T21:33:13.657824Z,1,0
7,-1420723.7204720487,2392605.17884686,6571200.573335108,-943.5541343871154,6899.644769246908,-2710.294727157773,5.17183818743846,0.3947519316561246,1611727.956312161,2017-08-26T21:34:13.657824Z,1,0
8,-1474507.3962201504,2801609.796061519,6395740.114385724,-848.6451731239442,6729.370223046121,-3136.475600247696,4.925472776184298,0.321193116115922,1801162.984455339,2017-08-26T21:35:13.657824Z,1,0
9,-1522494.1044519783,3199599.8462815173,6195069.546223336,-750.3848286978912,6532.611193826113,-3550.3557696280977,4.744477888443639,0.2392812180352601,2062672.3781960723,2017-08-26T21:36:13.657824Z,1,0
10,-1564494.2244736496,3585008.654716302,5969976.096402147,-649.157331156752,6310.131316025975,-3950.303971488127,4.6137595944906415,0.1612963679677057,2371647.45405773,2017-08-26T21:37:13.657824Z,1,0
11,-1600341.562780278,3956318.424188544,5721343.32346206,-545.3591245086104,6062.795729587452,-4334.74142867125,4.5178647105861085,0.0905496015477771,2710856.522162978,2017-08-26T21:38:13.657824Z,1,0
12,-1629894.0440601627,4312066.24181349,5450147.789704978,-439.397457233904,5791.568317703928,-4702.147805211405,4.445647359093886,0.0267784091369952,3069123.776063215,2017-08-26T21:39:13.657824Z,1,0
13,-696132.3917920114,-1597661.6836635729,6919285.775826914,-1608.8874860623707,7144.6294272525365,1484.8056841484413,0.5767766998748823,0.0271085758392875,3078298.970882717,2017-08-26T23:05:13.657824Z,1,1
14,-791234.4272978703,-1166125.5652030874,6994685.17540184,-1560.1427152299364,7235.1944431043585,1027.6876184282794,0.4718900159429863,0.0790713732902835,2781465.359862721,2017-08-26T23:06:13.657824Z,1,1
15,-883226.9841031705,-730007.1275438869,7042527.847717881,-1505.2720484615083,7297.32424611296,566.5479221786715,0.341469047947979,0.13073892789768,2519484.7202356146,2017-08-26T23:07:13.657824Z,1,1
16,-971748.8341233756,-291020.2354571408,7062626.590447333,-1444.4901917860554,7330.778750304818,103.19173371831216,0.1808142819563975,0.1784542344667629,2305512.22354457,2017-08-26T23:08:13.657824Z,1,1
17,-1056452.3178484223,149110.15254156175,7054902.796639451,-1378.034600819371,7335.42847009629,-360.56737212206497,6.272115808631328,0.2159547014694693,2154882.4551582383,2017-08-26T23:09:13.657824Z,1,1
18,-1137004.683678958,588654.646319555,7019386.748807686,-1306.164572873355,7311.255152549618,-822.9145035034211,6.056161894574287,0.235471703722036,2082046.63911852,2017-08-26T23:10:13.657824Z,1,1
19,-1213089.365302822,1025886.0897737368,6956217.4958604695,-1229.1601389788552,7258.351478610991,-1282.0400380882654,5.833014617537068,0.2316837533633272,2095422.489149412,2017-08-26T23:11:13.657824Z,1,1
20,-1284407.1958436847,1459086.2131847716,6865642.34774274,-1147.3211122492917,7176.920953935372,-1736.1460167625914,5.623061407695267,0.2056855367067862,2193363.6332462365,2017-08-26T23:12:13.657824Z,1,1
21,-1350677.5697817078,1886552.2698125204,6748015.970794838,-1060.9662270410136,7067.277346484064,-2183.452983979974,5.440438629743981,0.1641021786284968,2365003.9355249074,2017-08-26T23:13:13.657824Z,1,1
22,-1411639.5439605252,2306603.6197148818,6603799.078658338,-970.4319296697556,6929.843583098113,-2622.206549344014,5.289352114681659,0.1145361915948446,2595123.511453379,2017-08-26T23:14:13.657824Z,1,1
23,-1467052.859521723,2717588.2412932534,6433556.737384236,-876.0710530276593,6765.15043777943,-3050.683746716737,5.167204542144472,0.062440720406348,2868847.905399054,2017-08-26T23:15:13.657824Z,1,1
24,1035758.3725410592,-6777341.149011067,2007401.180302729,-1402.9294352852942,1896.1806495695496,7089.883287241839,2.0798208763548702,0.0327217215443939,3018318.2898378912,2017-08-27T00:28:13.657824Z,0,2
25,949598.5960003112,-6650301.808246001,2428551.7146481,-1468.120367417118,2337.0751457782403,6943.847548618184,1.9733490707794543,0.0862308066989948,2718675.262773524,2017-08-27T00:29:13.657824Z,0,2
26,859699.5790224236,-6497076.210545166,2840114.5956086465,-1527.5295525340882,2768.7669695765394,6770.400945861687,1.840570299209516,0.140092162000122,2452903.892919917,2017-08-27T00:30:13.657824Z,0,2
27,766415.4115960462,-6318268.118415394,3240465.360745417,-1580.9214496297827,3189.5429545898987,6570.234549850246,1.6759520354519055,0.1906107771659906,2234538.984782204,2017-08-27T00:31:13.657824Z,0,2
28,670113.5931200476,-6114582.775515221,3628024.242504981,-1628.0847930843065,3597.734695126401,6344.146674571956,1.477620642951078,0.2310669581244641,2079736.2940382576,2017-08-27T00:32:13.657824Z,0,2
29,571173.5420930161,-5886824.008527225,4001262.5132505097,-1668.8337278966128,3991.725736304247,6093.039700047162,1.2526413042832536,0.2527008286541543,2004119.2005095668,2017-08-27T00:33:13.657824Z,0,2
30,469985.0460937348,-5635890.8910913095,4358708.601031738,-1703.0085154255464,4369.958754130538,5817.915598082253,1.0194089417894423,0.2492356832168402,2017071.156917401,2017-08-27T00:34:13.657824Z,0,2
31,366946.6928697641,-5362774.010987807,4698953.912793815,-1730.475552910944,4730.941451754122,5519.870895243036,0.8006910810399775,0.2217762565472205,2117045.172626764,2017-08-27T00:35:13.657824Z,0,2
32,262464.2941778935,-5068551.418691527,5020658.376146743,-1751.127790182959,5073.251599065046,5200.092453116222,0.6120169805910727,0.1778503359676841,2292428.0539943567,2017-08-27T00:36:13.657824Z,0,2
33,156949.26335870908,-4754384.245395161,5322555.735309669,-1764.885691258149,5395.543187118886,4859.8530716394025,0.45749319578652,0.1259440246254154,2527039.6175031946,2017-08-27T00:37:13.657824Z,0,2
34,50816.95089040574,-4421511.93837082,5603458.525706097,-1771.6975260478846,5696.552560919848,4500.505167242737,0.3337614644763843,0.0719198717426268,2805318.2007556264,2017-08-27T00:38:13.657824Z,0,2
35,-55515.01661700645,-4071247.180640486,5862262.648013403,-1771.539118476421,5975.102826933691,4123.474123665547,0.2349847028205159,0.0187784067269895,3114677.4204661897,2017-08-27T00:39:13.657824Z,0,2
36,-474499.38148148014,-2524776.816123608,6657552.851032172,-1701.6764237862965,6844.412190977422,2469.2379830069253,0.8944008308811532,0.0533651207375533,2922050.0955155264,2017-08-27T00:43:13.657824Z,1,3
37,-575600.4982531521,-2109419.425801342,6792493.007486952,-1667.2573665082602,6996.294892953707,2027.2981373167445,0.7705994290620741,0.1019110512879034,2660238.841713846,2017-08-27T00:44:13.657824Z,1,3
38,-674439.2943266158,-1685771.7088266837,6900669.7850972125,-1626.2905943855342,7120.667185131708,1577.4171643255888,0.620014469217215,0.1469126462143614,2443202.567557712,2017-08-27T00:45:13.657824Z,1,3
39,-770627.6304188956,-1255499.271253761,6981659.517603364,-1578.9376687167298,7217.047717337939,1121.3594664241934,0.4415361501936363,0.1833306870805827,2284760.206330688,2017-08-27T00:46:13.657824Z,1,3
40,-863787.7920690558,-820293.3438970811,7035145.119092663,-1525.3844164671473,7285.063853192443,660.9122124283751,0.2402734433515899,0.2047733486757058,2198339.100957037,2017-08-27T00:47:13.657824Z,1,3
41,-953553.932733588,-381864.22727634886,7060917.252681859,-1465.8405585506614,7324.452758803254,197.8787450159006,0.0292694498415535,0.2062203309072847,2192847.652509128,2017-08-27T00:48:13.657824Z,1,3
42,-1039573.49149188,58065.32433841149,7058875.095478988,-1400.5390552130734,7335.062222958881,-265.92828132423097,6.108988878352045,0.1872988938627319,2268902.4860237963,2017-08-27T00:49:13.657824Z,1,3
43,-1121508.5533999267,497766.752436364,7029026.697593444,-1329.734956919842,7316.851045692236,-728.6935084748583,5.927024412621047,0.1525604216602979,2418537.7085250407,2017-08-27T00:50:13.657824Z,1,3
44,-1199037.1439023367,935512.3627030192,6971488.944540124,-1253.7044119276943,7269.889199721346,-1188.6055860449972,5.772722162078488,0.1084115331238563,2628670.287545812,2017-08-27T00:51:13.657824Z,1,3
45,-1271854.4685124017,1369581.9829371606,6886487.117204737,-1172.743794836988,7194.357361826329,-1643.863727909301,5.645605439965319,0.060152320636757,2885355.805105638,2017-08-27T00:52:13.657824Z,1,3
46,1291088.2378023218,-7001270.202073585,612880.3972467138,-1158.5291999456044,440.6556807318183,7367.292324939345,3.139153130620718,0.0227507437357251,3073005.6178752864,2017-08-27T02:05:13.657824Z,0,4
47,1219080.8124051562,-6961067.029374576,1053417.5434487006,-1240.9326321689691,899.0210945578057,7312.445745150037,3.1744905520122195,0.0928072097988228,2679026.52711356,2017-08-27T02:06:13.657824Z,0,4
48,1142273.8390821344,-6893457.369213049,1489796.328720865,-1318.4608659942803,1353.9020408789736,7228.725206578596,3.2223598553659945,0.1754438356104341,2291525.1484356057,2017-08-27T02:07:13.657824Z,0,4
49,1060969.216082859,-6798704.707840329,1920293.5764759025,-1390.8044172929772,1803.4875433982008,7116.458961532833,3.291229686695216,0.2771762664956642,1916411.2711825185,2017-08-27T02:08:13.657824Z,0,4
50,975486.7885782912,-6677180.577761494,2343209.307800057,-1457.6747281634457,2245.988345892875,6976.09172466172,3.39886330285238,0.4082133429532308,1564589.920248781,2017-08-27T02:09:13.657824Z,0,4
51,886163.0563133617,-6529363.01637246,2756873.6739540803,-1518.8054524197512,2679.6449162450044,6808.182729680557,3.5871408791738304,0.5804241906987899,1257650.754528684,2017-08-27T02:10:13.657824Z,0,4
52,793349.7952398912,-6355834.53563023,3159653.7393729542,-1573.9539499837915,3102.7355878291014,6613.402495655254,3.958220175075326,0.7800040735389753,1037815.3599717998,2017-08-27T02:11:13.657824Z,0,4
53,697412.6053434683,-6157279.639054243,3549960.030457077,-1622.9021231020258,3513.583263264571,6392.528571834212,4.617332607066203,0.873543317847827,967358.046634182,2017-08-27T02:12:13.657824Z,0,4
54,598729.4273988152,-5934481.980231351,3926252.854525414,-1665.4568143900435,3910.561024395399,6146.442349171753,5.221633818678467,0.7389280937719603,1076138.868457097,2017-08-27T02:13:13.657824Z,0,4
55,497689.0200597209,-5688321.156686255,4287048.448333874,-1701.4507783987567,4292.098906760995,5876.126132960576,5.540015642040721,0.5407337434768812,1320106.959558586,2017-08-27T02:14:13.657824Z,0,4
56,394689.377110779,-5419769.070968061,4630924.89907213,-1730.743476376043,4656.691143422365,5582.658360006769,5.705560785755184,0.3787579190740337,1639167.359855118,2017-08-27T02:15:13.657824Z,0,4
57,290136.1218363009,-5129885.900611596,4956527.743688935,-1753.2211252125269,5002.9018072246845,5267.208059682576,5.803466208976598,0.2556017058823361,1996607.2828733823,2017-08-27T02:16:13.657824Z,0,4
58,184440.897419514,-4819815.764226855,5262575.260763804,-1768.7969253721556,5329.369637859422,4931.030199456197,5.86812214112732,0.1590756553062487,2373975.350089655,2017-08-27T02:17:13.657824Z,0,4
59,78019.73234417598,-4490782.079205271,5547863.489696167,-1777.4115751945642,5634.81357570153,4575.460951864871,5.914486768243303,0.0798951357615428,2761842.624193181,2017-08-27T02:18:13.657824Z,0,4
60,-28708.61964691432,-4144082.57940612,5811270.928724926,-1779.0335523493602,5918.038310298036,4201.911670721593,1.549388208483289,0.0208939507094058,3107852.7351296367,2017-08-27T02:19:13.657824Z,1,4
61,-135324.3460914462,-3781084.030867239,6051762.843843219,-1773.6590498617772,6177.938725843701,3811.8621296230217,1.4617295432574915,0.0788463312075637,2775760.4460570426,2017-08-27T02:20:13.657824Z,1,4
62,-241408.22957013932,-3403216.7116976106,6268395.172210495,-1761.3116449506203,6413.503386992397,3406.8538252827057,1.3491575486734315,0.1401844854699125,2469529.2381889294,2017-08-27T02:21:13.657824Z,1,4
63,-346543.2795919147,-3011968.706099175,6460318.030939804,-1742.041934449891,6623.817254501834,2988.483577638898,1.2033243071936757,0.2029220108808491,2201518.703235484,2017-08-27T02:22:13.657824Z,1,4
64,-450316.3433888437,-2608880.0494746505,6626778.862171168,-1715.9271854222588,6808.0640971766925,2558.3978335725856,1.0161088163252945,0.2614847538716138,1988581.200222024,2017-08-27T02:23:13.657824Z,1,4
65,-552319.6957056168,-2195536.7152596377,6767125.245069744,-1683.0709993165449,6965.529384250727,2118.287051569815,0.7857871779202454,0.3050928055606676,1850857.802658276,2017-08-27T02:24:13.657824Z,1,4
66,-652152.6112406626,-1773564.426366744,6880807.3500010185,-1643.6030724182706,7095.603083369336,1669.8793870726931,0.5265476075925113,0.3205605073689299,1806273.93193709,2017-08-27T02:25:13.657824Z,1,4
67,-749422.9232849006,-1344622.313946575,6967380.000568332,-1597.6789330830154,7197.7819674630455,1214.933967231636,0.2691573669749056,0.3019163182531947,1861740.228120193,2017-08-27T02:26:13.657824Z,1,4
68,-843748.5579469486,-910396.4555036852,7026504.319013008,-1545.4793373993352,7271.671225171997,755.2337650738064,0.0432063261557927,0.2563458349842889,2008800.773938134,2017-08-27T02:27:13.657824Z,1,4
69,-934759.0187874832,-472593.3459019199,7057948.953419115,-1487.2092678976858,7316.985209614383,292.5789933531545,6.144043299065267,0.1971301514516607,2228895.905828837,2017-08-27T02:28:13.657824Z,1,4
70,-1022096.8109940996,-32933.31773298224,7061590.931342294,-1423.0969926618598,7333.548539179951,-171.21898327516902,6.002450486688392,0.1345060701680578,2502079.6819803235,2017-08-27T02:29:13.657824Z,1,4
71,-1105418.8136282992,406856.1084756657,7037416.1244415445,-1353.3932168316214,7321.297054345917,-634.3446026338348,5.8932691935970025,0.0735739412087846,2811993.478592078,2017-08-27T02:30:13.657824Z,1,4
72,212022.46554953928,-4884738.805693457,5201285.973611156,-1772.357566666575,5261.967489391264,5001.867035638744,2.2129652821558814,0.0784253304862942,2771117.8155375062,2017-08-27T03:57:13.657824Z,1,5
73,105333.58328064214,-4559617.905933,5490941.566124916,-1782.7719580709115,5571.83204542471,4650.154521530702,2.1596762708907087,0.1566509953122595,2388453.759614359,2017-08-27T03:58:13.657824Z,1,5
74,-1769.8827173165864,-4216558.487687281,5758938.907255516,-1786.172929102013,5859.723168716728,4280.163894660507,2.0829386068924083,0.2504691804945179,2020018.0331769928,2017-08-27T03:59:13.657824Z,1,5
75,-108866.66414034407,-3856913.1500431504,6004224.935821293,-1782.549604078076,6124.516157809523,3893.36051421188,1.965302220694761,0.3663165530660167,1677210.2786905405,2017-08-27T04:00:13.657824Z,1,5
76,-215535.65810176337,-3482099.2226186246,6225836.460672719,-1771.9183221344877,6365.178638368745,3491.273648499556,1.7713355433183982,0.508471641859623,1381154.6593953993,2017-08-27T04:01:13.657824Z,1,5
77,-321357.5525983564,-3093593.124368558,6422903.813204285,-1754.32243412685,6580.773879853085,3075.4905026101383,1.4361041793723783,0.6581714942685504,1169548.7829644324,2017-08-27T04:02:13.657824Z,1,5
78,-425916.4486066632,-2692924.52056412,6594654.15032476,-1729.832366896483,6770.464374538865,2647.650326834313,0.9249347017533808,0.730225777931231,1093625.638522032,2017-08-27T04:03:13.657824Z,1,5
79,-528801.4834416929,-2281670.260570207,6740414.37283641,-1698.545528950891,6933.515313584184,2209.437421632439,0.4221614042825011,0.6504816839290172,1180176.994417761,2017-08-27T04:04:13.657824Z,1,5
80,-629608.4397667112,-1861448.1428753985,6859613.603210989,-1660.5859395566606,7069.29678653761,1762.5736464769295,0.0982662448037246,0.5002040693035756,1399090.4942288236,2017-08-27T04:05:13.657824Z,1,5
81,-727941.3226018763,-1433910.58698835,6951785.227138666,-1616.1035003126065,7177.284771154173,1308.8115450721716,6.194217608878329,0.3600247151664408,1699304.3818685485,2017-08-27T04:06:13.657824Z,1,5
82,-823413.883937054,-1000738.2577441444,7016568.562276012,-1565.2730831938857,7257.062009770266,849.9287848105339,6.08040426141701,0.2460298806330614,2044402.6017182893,2017-08-27T04:07:13.657824Z,1,5
83,-915651.1034792658,-563633.6059784305,7053710.203161258,-1508.2940083407775,7308.319925575672,387.7225335858176,6.006127513171457,0.1535412531214787,2414133.995407463,2017-08-27T04:08:13.657824Z,1,5
84,-1004290.6399005662,-124314.29101069672,7063064.987044321,-1445.389440500653,7330.860324618609,-75.99749921318268,5.954680993120701,0.0762303118538199,2797555.597033437,2017-08-27T04:09:13.657824Z,1,5
85,1265547.820245828,-6977719.477265529,872801.9563443242,-1209.429804232511,708.2334223978027,7338.723166996189,2.120336593748216,0.0259835847465827,3056922.461831155,2017-08-27T05:26:13.657824Z,2,6
86,1190537.994997691,-6921516.15677078,1311112.71017648,-1290.0785783275849,1164.6044572568371,7266.825510789536,2.012317278655713,0.0772096124907912,2764513.397870727,2017-08-27T05:27:13.657824Z,2,6
87,1110839.847565106,-6838057.8005315615,1744246.9988725658,-1365.6552007380783,1616.433192558015,7166.230808097901,1.879181327041193,0.1281240869951137,2506382.3335967534,2017-08-27T05:28:13.657824Z,2,6
88,1026766.8591618232,-6727670.995925953,2170494.329419328,-1435.8581425721902,2061.921694027763,7037.336054026964,1.7163405667800444,0.1751320062867401,2295516.8799180877,2017-08-27T05:29:13.657824Z,2,6
89,938649.9517031468,-6590789.457644702,2588171.4872718006,-1500.4077902255706,2499.2983013860094,6880.653133242695,1.5228651247964058,0.2121085934091417,2147041.917307486,2017-08-27T05:30:13.657824Z,2,6
90,846836.1380127529,-6427952.223977516,2995629.3505860786,-1559.047552344195,2926.8250502412457,6696.8061257573,1.305802563537245,0.2314711307786684,2075228.1548060656,2017-08-27T05:31:13.657824Z,2,6
91,751687.1120846664,-6239801.426852599,3391259.542563408,-1611.5448161259108,3342.80463753509,6486.528684736204,1.0816897853951382,0.2279966233286829,2088469.0617197736,2017-08-27T05:32:13.657824Z,2,6
92,653577.7792113245,-6027079.633758762,3773500.9199770456,-1657.6919882706577,3745.587757503848,6250.661027354617,0.87055574465919,0.2026726595191541,2185323.81842237,2017-08-27T05:33:13.657824Z,2,6
93,552894.7265938474,-5790626.746016362,4140845.851669694,-1697.3073904492292,4133.580447458918,5990.145913424309,0.6864229334974267,0.1618978555834595,2355279.091005046,2017-08-27T05:34:13.657824Z,2,6
94,450034.648778719,-5531376.472937954,4491846.2248697,-1730.235868993102,4505.250742070331,5706.0236382328185,2.9974071371065145,0.0284225278450466,3051017.2260597576,2017-08-27T05:35:13.657824Z,1,6
95,450034.648778719,-5531376.472937954,4491846.2248697,-1730.235868993102,4505.250742070331,5706.0236382328185,0.5336125138902292,0.1130991768512294,2583458.8532177475,2017-08-27T05:35:13.657824Z,2,6
96,345402.7436125809,-5250352.443667796,4825119.139717215,-1756.349191421572,4859.134122275767,5399.426861514974,3.0237958477008973,0.1002243118450129,2653231.786450115,2017-08-27T05:36:13.657824Z,1,6
97,345402.7436125809,-5250352.443667796,4825119.139717215,-1756.349191421572,4859.134122275767,5399.426861514974,0.409700099791901,0.0616582277770588,2855235.77025771,2017-08-27T05:36:13.657824Z,2,6
98,239411.0869087472,-4948664.005757899,5139352.309319224,-1775.5463750264,5193.838728358949,5071.575886227563,3.058837812801426,0.1858513542201758,2260358.0425871536,2017-08-27T05:37:13.657824Z,1,6
99,239411.0869087472,-4948664.005757899,5139352.309319224,-1775.5463750264,5193.838728358949,5071.575886227563,0.3096982857738736,0.0105120396096109,3158508.313991748,2017-08-27T05:37:13.657824Z,2,6
100,132476.99015963587,-4627501.706637615,5433309.165729655,-1787.7538989818306,5508.050559162531,4723.773466679243,3.1087596459601143,0.2932511074536035,1877499.53241409,2017-08-27T05:38:13.657824Z,1,6
101,25021.34984819226,-4288132.484320443,5705833.648031051,-1792.925808595147,5800.538015776506,4357.399590999438,3.187479604444827,0.4364389805332222,1514317.139631533,2017-08-27T05:39:13.657824Z,1,6
102,-82533.01247153309,-3931894.580492152,5955854.684196963,-1791.0439294885632,6070.156755973403,3973.9062845227177,3.3325954296497167,0.6381911949922403,1190970.7931994775,2017-08-27T05:40:13.657824Z,1,6
103,-189763.0219858464,-3560192.159474568,6182390.329814166,-1782.1179929200048,6315.854391615049,3574.811334192193,3.67412161240394,0.9099113003275378,950825.1131490784,2017-08-27T05:41:13.657824Z,1,6
104,-296247.0234626838,-3174489.675704932,6384551.517005492,-1766.1855736208165,6536.674226088032,3161.6917822670266,4.557143690878405,1.0763892752963748,867166.3005044294,2017-08-27T05:42:13.657824Z,1,6
105,-401566.4482062717,-2776306.026080383,6561545.417470762,-1743.311925595498,6731.758788948,2736.1774401288603,5.346935488324667,0.8663593502415249,981088.9094594636,2017-08-27T05:43:13.657824Z,1,6
106,-505307.4575463311,-2367208.5001097987,6712678.396910877,-1713.589288696372,6900.352865323504,2299.943730401087,5.640415084301143,0.603864418465801,1238825.2568329035,2017-08-27T05:44:13.657824Z,1,6
107,-607062.529341667,-1948806.5877659733,6837358.535334338,-1677.1358218162509,7041.805159423598,1854.704620079783,5.770964945777164,0.413605318721885,1570329.953242737,2017-08-27T05:45:13.657824Z,1,6
108,-706431.9944039163,-1522745.713779295,6935097.758334265,-1634.0951148974807,7155.569379227109,1402.206767399715,5.844228808300451,0.2776569319767175,1937098.129411721,2017-08-27T05:46:13.657824Z,1,6
109,-803025.5606034243,-1090700.883239488,7005513.632223813,-1584.6361109886086,7241.206182531113,944.223828718476,5.891920580285829,0.1746405296958884,2321332.043866805,2017-08-27T05:47:13.657824Z,1,6
110,-896463.8216180959,-654370.2037025016,7048330.785644498,-1528.9526036165796,7298.385056890773,482.5497340095021,5.926190286768624,0.0917856282744106,2714387.7624785127,2017-08-27T05:48:13.657824Z,1,6
111,-986379.7244193936,-215468.3214682575,7063381.921690346,-1467.2624393843878,7326.885367712362,18.991938677951687,5.952606077141354,0.021839343774465,3111624.5555944466,2017-08-27T05:49:13.657824Z,1,6
112,1479048.883477582,-6970660.441682713,-546115.9366797627,-920.0524659904634,-764.472013765114,7373.301407769612,3.183042871443845,0.0162624057846465,3113149.892252409,2017-08-27T07:03:13.657824Z,2,7
113,1420971.7592723158,-7002783.439930312,-102931.51004666496,-1015.2204684139496,-305.92424730870084,7394.654498907963,3.2245241644493974,0.0843468944961994,2723305.012730168,2017-08-27T07:04:13.657824Z,2,7
114,1357301.928719996,-7007345.166211559,340658.9627982064,-1106.4139271333383,153.93253466467888,7386.829535538975,3.2804616499393706,0.1637256729692457,2341115.2788168355,2017-08-27T07:05:13.657824Z,2,7
115,1288288.875920237,-6984321.940456451,782904.7751149649,-1193.2685719816343,613.2691180803854,7349.845588598687,3.36001569606212,0.2596323956297913,1972963.687220223,2017-08-27T07:06:13.657824Z,2,7
116,1214203.4261985722,-6933799.845588213,1222059.9368132183,-1275.437535717726,1070.257330993566,7283.840757441291,3.481262173488936,0.379336735134801,1630289.6029534484,2017-08-27T07:07:13.657824Z,2,7
117,1135336.6607927924,-6855974.419277707,1656390.2990358584,-1352.5927146585625,1523.0781894167078,7189.0714796884,3.682192173046416,0.5283431469807196,1334652.5872132997,2017-08-27T07:08:13.657824Z,2,7
118,1051998.7487414218,-6751149.86373253,2084180.61635218,-1424.4262005015937,1969.9298245948671,7065.91112691355,4.036129456373028,0.6868867156970874,1125242.56734765,2017-08-27T07:09:13.657824Z,2,7
119,964517.6913924944,-6619737.792871555,2503741.5019021723,-1490.6517168387145,2409.0348759275794,6914.847921536171,4.583623289368976,0.7582028766543661,1055557.2808241046,2017-08-27T07:10:13.657824Z,2,7
120,873237.9859506501,-6462255.5544901285,2913416.2598212725,-1551.0058698350963,2838.6474367135,6736.483067557064,5.106726852950873,0.6640345706284846,1151361.1105537135,2017-08-27T07:11:13.657824Z,2,7
121,778519.2199416595,-6279324.117162184,3311587.6094506863,-1605.2492077169393,3257.060740461698,6531.528781059332,5.430231833400738,0.5035045251780526,1378120.7050622553,2017-08-27T07:12:13.657824Z,2,7
122,680734.6047542145,-6071665.477493783,3696684.2541563977,-1653.1672229017513,3662.6152824968112,6300.80495684723,5.6151404585793,0.3595292129306994,1682920.8395609658,2017-08-27T07:13:13.657824Z,2,7
123,580269.4506768897,-5840099.595083704,4067187.203969052,-1694.5712400002417,4053.706254172032,6045.234367860908,5.728966764502592,0.2445191734359302,2029889.106894709,2017-08-27T07:14:13.657824Z,2,7
124,477519.60132379545,-5585540.92772508,4421635.790459676,-1729.2988796179286,4428.789375898683,5765.8375044370905,5.805338534989761,0.1519853082313573,2399635.470762946,2017-08-27T07:15:13.657824Z,2,7
125,372889.84335936286,-5308994.645411275,4758633.398063041,-1757.2145193667377,4786.386109062984,5463.728280539817,3.8262701709381095,0.024762764476744,3075395.9487462672,2017-08-27T07:16:13.657824Z,1,7
126,372889.84335936286,-5308994.645411275,4758633.398063041,-1757.2145193667377,4786.386109062984,5463.728280539817,5.860361984286388,0.0749625503291114,2781825.969614091,2017-08-27T07:16:13.657824Z,2,7
127,266792.2837934009,-5011552.517573171,5076852.935359145,-1778.209839637402,5125.089237244297,5140.10941423252,3.9220798151781167,0.0808351060371316,2757023.966156984,2017-08-27T07:17:13.657824Z,1,7
128,266792.2837934009,-5011552.517573171,5076852.935359145,-1778.209839637402,5125.089237244297,5140.10941423252,5.9023080226395335,0.0085699252590014,3170632.3427366503,2017-08-27T07:17:13.657824Z,2,7
129,159644.7012392035,-4694388.479888575,5375042.015527859,-1792.2041635567527,5443.567953889345,4796.267578989386,4.041998317003673,0.1397268733562395,2464519.1145635904,2017-08-27T07:18:13.657824Z,1,7
130,51868.87627040412,-4358753.885322175,5652027.849964534,-1799.144818592841,5740.573474111445,4433.568411855634,4.193422957714214,0.1993623741287082,2209439.650038583,2017-08-27T07:19:13.657824Z,1,7
131,-56111.09296953393,-4005972.420138218,5906721.796964833,-1799.0072049708133,6014.944162584598,4053.450058666,4.383126559836619,0.254471257220058,2007163.860467131,2017-08-27T07:20:13.657824Z,1,7
132,-163870.47207229072,-3637434.7596720606,6138123.523545324,-1791.7947076588855,6265.609331602045,3657.417042706239,4.611905124480507,0.2954482348236144,1875629.1309824567,2017-08-27T07:21:13.657824Z,1,7
133,-270985.565115019,-3254592.990533316,6345324.807040704,-1777.5384711867125,6491.59325540915,3247.0340503669763,4.866972736075409,0.3109034931269542,1830571.7099543056,2017-08-27T07:22:13.657824Z,1,7
134,-377035.3603238378,-2858954.817213111,6527512.942795382,-1756.2967104185566,6692.018119380259,2823.9193215633472,5.121670424868438,0.2954515921374597,1878275.7516322532,2017-08-27T07:23:13.657824Z,1,7
135,-481603.1414237222,-2452077.610286341,6683973.785295056,-1728.154396304758,6866.106747254169,2389.73894967643,5.349610475641237,0.2546986050826237,2011867.4029911412,2017-08-27T07:24:13.657824Z,1,7
136,-584278.0933550715,-2035562.271117792,6814094.426197263,-1693.2232339035702,7013.185714058148,1946.2001766508477,5.53849421212096,0.2000352998160708,2215304.496443138,2017-08-27T07:25:13.657824Z,1,7
137,-684656.8959305326,-1611046.9469183788,6917365.444462818,-1651.641249239982,7132.687127708159,1495.044225689177,5.689414460658134,0.1409467805017026,2470626.983367701,2017-08-27T07:26:13.657824Z,1,7
138,-782345.2839476723,-1180200.6588192363,6993382.7531783935,-1603.5721293293327,7224.150005028607,1038.039980983615,5.80922251694844,0.0826089757201759,2762564.907601641,2017-08-27T07:27:13.657824Z,1,7
139,-876959.5648680559,-744716.8417656654,7041849.070317031,-1549.2044871566245,7287.221762968344,576.9776317471247,5.905293712655265,0.027059205918703,3079687.9881675104,2017-08-27T07:28:13.657824Z,1,7
140,-137791.05970145788,-3713915.4873063737,6092933.548669801,-1801.0643966138225,6214.401968447623,3739.2729176128537,4.938085227182425,0.0249892647093509,3084793.606126934,2017-08-27T09:01:13.657824Z,1,8
141,-245513.245319119,-3333993.6415216844,6305131.679005999,-1788.4982089138737,6445.501302557409,3331.685105035034,5.083133130220539,0.0414572731659506,2988078.8473082925,2017-08-27T09:02:13.657824Z,1,8
142,-352270.2174535589,-2940963.004266342,6492474.447519068,-1768.901566668403,6651.220330737861,2911.033219148947,5.235193813643194,0.0496249183083077,2942026.7662604945,2017-08-27T09:03:13.657824Z,1,8
143,-457642.3922606952,-2536370.8225806537,6654227.499680898,-1742.3537677378984,6830.761760009273,2478.9744445619485,5.389264902594126,0.0486205021031354,2949088.698038944,2017-08-27T09:04:13.657824Z,1,8
144,-561215.7524340661,-2121809.0997582884,6789757.198633085,-1708.960991273609,6983.43134672443,2037.2077269278184,5.539916304011146,0.0385883118283948,3008767.7992689083,2017-08-27T09:05:13.657824Z,1,8
145,-662583.4455836449,-1698908.358624406,6898532.93074814,-1668.855791343127,7108.639591759544,1587.4671707684038,5.682474890840112,0.0206325140227236,3117735.451541536,2017-08-27T09:06:13.657824Z,1,8
146,-1724107.6100858022,5003997.409715902,4791873.719478032,-258.72432266558496,5125.85476751838,-5431.538680328174,0.9609316262517682,0.0190491284051834,3110488.063325463,2017-08-27T12:43:13.657824Z,0,9
147,-1736228.6738100254,5301500.602848212,4456742.618699767,-145.1756747080364,4787.656889501391,-5735.83077862828,1.0935979181462512,0.0512125075844868,2918547.3617994203,2017-08-27T12:44:13.657824Z,0,9
148,-1741516.9303709816,5578138.404000148,4104025.9446926746,-31.039088669316904,4430.569962532962,-6017.529875805579,1.2426315898261,0.0773481505322588,2771989.04778215,2017-08-27T12:45:13.657824Z,0,9
149,-1739950.6283683907,5832819.316399192,3735113.405911972,83.23425980194682,4056.001353873829,-6275.513479605449,1.4053713262294714,0.094681467181346,2679100.552915829,2017-08-27T12:46:13.657824Z,0,9
150,-1731535.131928014,6064538.480812722,3351459.280866178,197.1919965522488,3665.4306573007375,-6508.752693249968,1.576169022551674,0.1008533782530596,2646132.6048157285,2017-08-27T12:47:13.657824Z,0,9
151,-1716302.9353202507,6272381.82345711,2954576.642892917,310.3825002552121,3260.4033490178663,-6716.317322301596,1.747194139433475,0.0948723930547391,2675639.616086564,2017-08-27T12:48:13.657824Z,0,9
152,-1694313.5559580969,6455529.811602893,2546031.313977267,422.3569613819844,2842.524239397856,-6897.379913147118,1.9105047712847547,0.0775910641563364,2765708.897786181,2017-08-27T12:49:13.657824Z,0,9
153,-1665653.3063557334,6613260.813057279,2127435.588730189,532.6713104983378,2413.4506918844586,-7051.219368798915,2.0601894552272366,0.0512954623528591,2910575.4083868302,2017-08-27T12:50:13.657824Z,0,9
154,-1630434.9627376408,6744954.040858742,1700441.7668487977,640.8877123718137,1974.88577284828,-7177.223694620839,2.193310449334154,0.0187625995156722,3102217.1762043573,2017-08-27T12:51:13.657824Z,0,9

rhooper · November 5, 2024, 1:50pm

If you save the data as a csv, you can run the code using this:

from ObserveStation import ObserveStation
from ResidentSpaceObject import ResidentSpaceObject

stations: List[ObserveStation] = [ObserveStation("Eglin"),ObserveStation("Clear"),ObserveStation("Kaena Point")]
RSO: ResidentSpaceObject = ResidentSpaceObject("ENVISAT")
estimationInterval = -12.0*60.0*60.0
observationFilePath = "PATH/TO/THE/DATA/I/GAVE"
estimationFilePath = "PATH/TO/WHERE/THE/RESULTS/SHOULD/GO"
coreNumberForEstimation = 1
sigmaErrorScale = 1 # 1 if including error, 0 if not
flagEstimateDragCoefficient = True
flagEstimateReflectionCoefficient = False
atmospherModelForEstimationAndPrediction = 'NRLMSISE00'
SimulateEstimation().run(stations,RSO,estimationInterval,observationFilePath,estimationFilePath,coreNumberForEstimation,sigmaErrorScale,flagEstimateDragCoefficient,flagEstimateReflectionCoefficient,atmospherModelForEstimationAndPrediction)

MaximeJ · November 6, 2024, 5:21pm

Hi @rhooper,

That’s a lot of code, I can’t promise that I’ll have time to dig into it very soon.

I don’t know if you and your professor are aware but we have a section “Scientific Publications Using Orekit” where we list papers using Orekit (in the Publications section of the website).
Would you like to add the two papers you mentionned there?

Cheers,
Maxime

rhooper · November 6, 2024, 9:00pm

Yes it is a lot of code. I was expecting it to be a little too much to post in the forum. But the good news is I found the problem. I was multiplying my range measurements by the noise, instead of adding the noise. It’s here from my original post:

# Read range and add error
two_way = False
self._tracks[indexMeasurement].measurements.append(Range(
                station._groundStationForEstimation,
                two_way,
                AbsoluteDate(row["date"],Common.utc()),
          >>>>> row["range[m]"]*random.gauss(0.0,1.0)*station._rangeSigma*sigmaErrorOn,
                station._rangeSigma,
                1.0,
                self.observableSatellite
            ))

I’ll reply to the Publications thread and mention the professor’s papers in there.

Thank you for all your help.

MaximeJ · November 7, 2024, 10:18am

Hi @rhooper,

Nice to see that you managed to fix the issue! I didn’t catch it when I read your first post.

Great, thanks!

Cheers,
Maxime