Ranging Data Generation

Orekit Python Wrapper usage
1

Dear All,

I’m writing a script to generate Ranging data for a Geostationary mission, in order to later evaluate my orbit determination performance. So far I have been successful to generate data for Range and RangeRate, but getting stuck with generating Angular AzEl. Below is my code:

import orekit
vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()
from orekit import JArray, JArray_double

# Measurements Builders
from org.orekit.estimation.measurements.generation import Generator
from org.orekit.estimation.measurements.generation import RangeBuilder
from org.orekit.estimation.measurements.generation import RangeRateBuilder
from org.orekit.estimation.measurements.generation import AngularAzElBuilder

# Modifiers - Measurements Perturbations
from org.orekit.estimation.measurements.modifiers import RangeTroposphericDelayModifier
from org.orekit.estimation.measurements.modifiers import RangeRateTroposphericDelayModifier

from org.orekit.estimation.measurements.modifiers import RangeIonosphericDelayModifier
from org.orekit.estimation.measurements.modifiers import RangeRateIonosphericDelayModifier

from org.orekit.estimation.measurements.modifiers import AngularTroposphericDelayModifier
from org.orekit.estimation.measurements.modifiers import AngularIonosphericDelayModifier

from org.orekit.estimation.measurements.modifiers import TurnAroundRangeIonosphericDelayModifier
from org.orekit.estimation.measurements.modifiers import TurnAroundRangeTroposphericDelayModifier

from org.orekit.estimation.measurements import Range
from org.orekit.estimation.measurements import RangeRate
from org.orekit.estimation.measurements import AngularAzEl
from org.orekit.estimation.measurements import GroundStation
from org.orekit.estimation.measurements import ObservableSatellite

from utils import *
from Propagator import *
from MissionDefinitions import *

import sys
import numpy as np
from datetime import datetime
import matplotlib.pyplot as plt
from math import degrees, radians

def generate_measurements(propagator, station, station_name, altitude, elevAngle, meas_type, t0, duration,
sigma          = None,
base_weight    = 1.0,
two_way        = True,
withRefraction = True,
step           = 3600.0,
seed           = 0):

    tf = t0.shiftedBy(3600.0*duration)

    noise_source = noise_generator(seed, sigma = sigma)
    generator  = Generator()
    satellite = generator.addPropagator(propagator)

    if meas_type == 'RANGE':
        builder = RangeBuilder(noise_source, station, two_way, sigma, base_weight, satellite)
    elif meas_type == 'RANGERATE':
            builder = RangeRateBuilder(noise_source, station, two_way, sigma, base_weight, satellite)
    elif meas_type == 'AZEL':
            builder = AngularAzElBuilder(noise_source, station, [sigma, sigma], [base_weight, base_weight], satellite)

    fixed_step_selector = FixedStepSelector(step, TimeScalesFactory.getUTC())
    if withRefraction:
            elevation_detector = ElevationDetector(station.getBaseFrame()).withConstantElevation(elevAngle).withRefraction(EarthITU453AtmosphereRefraction(altitude)).withHandler(ContinueOnEvent())
    else:
            elevation_detector = ElevationDetector(station.getBaseFrame()).withConstantElevation(elevAngle).withHandler(ContinueOnEvent())
    scheduler = EventBasedScheduler(builder, fixed_step_selector, propagator, elevation_detector, SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE)
    generator.addScheduler(scheduler)

    measurements = generator.generate(t0, tf)
    for measObject in measurements:
            if meas_type == 'RANGE':
                    meas = Range.cast_(measObject)
            elif meas_type == 'RANGERATE':
                    meas = RangeRate.cast_(measObject)
            elif meas_type == 'AZEL':
                    meas = AngularAzEl.cast_(measObject)
            else:
                    raise ValueError("unrecognized measurement type")

            if two_way:
                    way_str = "TWOWAY"
            else:
                    way_str = "ONEWAY"
            print("{}   {} {}       {}        {}".format(meas.getDate().toString(), way_str, meas_type, station_name, meas.getObservedValue()[0]))

if __name__ == '__main__':

utc = TimeScalesFactory.getUTC()

# Frame and Earth definitions
itrf93 = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth  = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf93)

# Beginning Of Life (BOL)
date = datetime(2022, 2, 1, 10, 30, 0)  # (YYYY,M,D,H,M,S)
# Convert to OREKIT Date
epochDate   = AbsoluteDate(date.year,date.month,date.day,date.hour,date.minute,float(date.second),utc)
initialDate = epochDate

# Initial State
gravitySpacecraft = Spacecraft("GravitySpace", 21.86, 0.29731552, 2.2, 1.8, 11111)
gravityOrbit      = Orbit(814, 786, 98.55, 90.0, 5.1917, 0.0567634)

# Setup numerical propagator
gravityProp = numericalPropagator(initialDate, gravityOrbit.sma(), gravityOrbit.eccentricity(),
	gravityOrbit.inclination, gravityOrbit.aop, gravityOrbit.raan, gravityOrbit.lv,
	gravitySpacecraft.mass, gravitySpacecraft.cross_section, gravitySpacecraft.cr, 70, 70)

# Define SSC Ground Stations
alaska      = GroundSegment("Alaska", +64.8000, -147.5000, 6.1905, 5.0, 295.1, 0.55, 1013.25, 22222)
puntaArenas = GroundSegment("PuntaArenas", -52.9400, -70.8600, 0.0340, 5.0, 295.1, 0.55, 1013.25, 33333)
esrange     = GroundSegment("Esrange", +67.8800, +21.0500, 0.5300, 5.0, 295.1, 0.55, 1013.25, 44444)
dongara     = GroundSegment("Dongara", -29.0500, +115.3500, 0.0340, 5.0, 295.1, 0.55, 1013.25, 55555)

stations    = (alaska, puntaArenas, esrange, dongara)
#stations, stationNames = orekitSpaceInventorStations(stations)

# Measurements Errors - Applicable to all ground stations from SSC
sscSigma = StationError(0.080, 15.000, 0.040) # [deg, deg, m, m.s-1]

# Ranging Data Generation
base_weight = 1.0
seed = int(sys.argv[1])
step = 60*float(sys.argv[2])
if sys.argv[3] == '1':
	rangeTwoWay = True
else:
	rangeTwoWay = False
if sys.argv[4] == '1':
	rateTwoWay = True
else:
	rateTwoWay = False

for station in stations:

	generate_measurements(gravityProp, station.groundStation(earth), station.name, station.altitude, station.elevationAngle,
		'RANGE', initialDate, 24.0*10, sigma=sscSigma.rangeSigma, two_way=rangeTwoWay, withRefraction=True, step=step, seed=seed)
	generate_measurements(gravityProp, station.groundStation(earth), station.name, station.altitude, station.elevationAngle,
		'RANGERATE', initialDate, 24.0*10, sigma=sscSigma.rangeRateSigma, two_way=rangeTwoWay, withRefraction=True, step=step, seed=seed)
	generate_measurements(gravityProp, station.groundStation(earth), station.name, station.altitude, station.elevationAngle,
		'AZEL', initialDate, 24.0, sigma=sscSigma.azelSigma, two_way=rangeTwoWay, withRefraction=True, step=step, seed=seed)

There are some functions that I’ve defined and that are not shown in the code, such as the noise generation and the GroundStation, but it’s not where I’m getting the errors. As I mentioned, the computation for Range and RangeRate work as expected. But for the AngularAzEl I always get the same error:

org.orekit.estimation.measurements.generation.AngularAzElBuilder@55322aab
Traceback (most recent call last):
    File "rangingDataSimulator.py", line 153, in <module>
        generate_measurements(gravityProp, station.groundStation(earth), station.name, station.altitude, station.elevationAngle,
    File "rangingDataSimulator.py", line 81, in generate_measurements
measurements = generator.generate(t0, tf)
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
Java stacktrace:
java.lang.ArrayIndexOutOfBoundsException: 1
    at org.orekit.estimation.measurements.generation.AngularAzElBuilder.build(AngularAzElBuilder.java:87)
    at org.orekit.estimation.measurements.generation.AngularAzElBuilder.build(AngularAzElBuilder.java:33)
    at org.orekit.estimation.measurements.generation.EventBasedScheduler.generate(EventBasedScheduler.java:124)
    at org.orekit.estimation.measurements.generation.Generator$GeneratorHandler.handleStep(Generator.java:130)
    at org.orekit.propagation.PropagatorsParallelizer$SinglePropagatorHandler.handleStep(PropagatorsParallelizer.java:322)
    at org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedStepHandler.handleStep(AbstractIntegratedPropagator.java:860)
    at org.hipparchus.ode.AbstractIntegrator.acceptStep(AbstractIntegrator.java:429)
    at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:290)
    at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:469)
    at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:414)
    at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:141)
    at org.orekit.estimation.measurements.generation.Generator.generate(Generator.java:95)

I don’t get errors when defining the “builder” as AngularAzEl, but when I call:

measurements = generator.generate(t0, tf)

I have already tried defining as a tuple, like:

az, el = generator.generate(t0, tf)

since it’s actually two measurements, but still no success.

I would be most grateful if someone could help me figuring out this one.

Many thanks in advance.

Regards,
/Leonardo.

2

Hi Leonardo,

I am not able to run the code directly as there are some external dependencies, but the ArrayIndexOutOfBoundsException error could be a lead.

Maybe as a first step, try to explicitly give the lists as JArray_double([value1, value2])?

Regards
/Petrus

1 Like
3

Dear Petrus,

Thank you for the prompt answer.

I am attaching the complete code with the functions I’ve defined in the other scripts (most definitey some libraries are repeated, since I just added everything together, but it runs now).

I tried defining the sigma and weights for the AngularAzEl as:

sigmaAzEl  = JArray_double([sigma, sigma])
weightAzEl = JArray_double([base_weight, base_weight])
builder = AngularAzElBuilder(noise_source, station, sigmaAzEl, weightAzEl, satellite)

But then I get the error of wrong inputs to the function. Was that your suggestion?

Thanks again.
/Leonardo.

orekitShow.py (14.8 KB)

4

Hi,

I think your covariance matrix needs to be two dimensional,

covariance = MatrixUtils.createRealDiagonalMatrix([float(sigma * sigma),float(sigma * sigma)])

Regards
/Petrus

1 Like
5

Hello Petrus,

Thank you very much! The code works now :smiley:, and I gotta say that I would have taken a relatively long time to find that one mistake.

Thanks a lot for the prompt support.

Best regards,
/Leonardo.