Bias on measurement angles computed from CPF

[nrol]

Hello everyone

I’m currently trying to compare some of my angular measurements to a reference from the ILRS but I have some unexpected bias that seems to be present so I’m trying to figure if it comes from a mistake somewhere in my code or if it is really due to the measures.

At first, I read the CPF file :

	public static Map<String,CPF.CPFEphemeris> readCPFFile(String filepath) {
		
		// Number of samples to use (does not change much the results)
		int interpolationSamples = 20;
		
		// Frames
		final CelestialBodies celestialBodies = DataContext.getDefault().getCelestialBodies();
		Frames frames = Frames.of(TimeScalesFactory.getTimeScales(), celestialBodies);
		
		// Configure the CPF Parser
		CPFParser cpfParser = new CPFParser(Constants.IERS2010_EARTH_MU,
                interpolationSamples,
                IERSConventions.IERS_2010,
                TimeScalesFactory.getUTC(),
                frames);
		
		// Parse of the ILRS file
		CPF cpfData = cpfParser.parse(new DataSource(filepath));
		
		// Get the satellite Data
		return cpfData.getSatellites();
	}

With that, I can now get a propagator and a :

// Ephéméride lié au satellite recherché
CPFEphemeris ephemeris = mapCpfData.get(satName);
			
// Propagateur lié au satellite
BoundedPropagator propagator = ephemeris.getPropagator();

And now I compute the angles :

// Propagation at the desired date
SpacecraftState state = propagator.propagate(epoch);

// Position of the station in ITRF
final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
	Constants.IERS2010_EARTH_FLATTENING,
	FramesFactory.getITRF(ITRFVersion.ITRF_2020, IERSConventions.IERS_2010, false));
final TopocentricFrame topo = new TopocentricFrame(earth, coordStation, "topo");
final GroundStation station = new GroundStation(topo);

// Position of the station in ICRF
Vector3D stationIcrf = station.getBaseFrame().getPVCoordinates(epoch, FramesFactory.getICRF()).getPosition();

// Vector station -> satellite
Vector3D diffPosition = state.getOrbit().getPVCoordinates(FramesFactory.getICRF()).getPosition().subtract(stationIcrf);
		
// Refraction model
EarthITU453AtmosphereRefraction refractionModel = new EarthITU453AtmosphereRefraction(coordStation.getAltitude());
		
// Measurement angles
double rightAscension = Math.atan2(diffPosition.getY(), diffPosition.getX());
double declination = Math.asin(diffPosition.getZ() / diffPosition.getNorm());
double azimuth = topo.getAzimuth(state.getPVCoordinates().getPosition(), state.getFrame(), state.getDate());
double elevation = topo.getElevation(state.getPVCoordinates().getPosition(), state.getFrame(), state.getDate());
double refractionAngle = refractionModel.getRefraction(elevation);

Is there anywhere I did some mistake ?
Since the biais is mostly on the right ascension, it might be due to some frame incoherence ?
The CPF I use is this one :
ilrs_data.cpf (39.1 KB)

If you have any idea fr improvement, that would be great !

Have a nice day!

[Serrof]

Hey there,

you ought to define an AngularRaDecBuilder (and same thing for the Azimuth-Elevation/Altitude angles), as from a quick look I believe you are missing the light time delay.
Also how are your measurements obtained? Depending on that, you might need to discard the refraction and apply other corrections instead.

Best,
Romain.

[nrol]

Thanks for the reply!
Good point on the AngularRaDecBuilder, I thought of the light time delay but didn’t think of this way to take it into account. I’ll look at those classes! Would you happend to have an example on hand, by any chance ?
My measurements are optic measurements, computed first in RaDec (compared to stars) and then computed in AzEl with and without refraction. I computed the refraction since I didn’t find in the documentation if the elevation given by Orekit was the true or apparent elevation
So RaDec angles are the values I’m comparing at first and AzEl are a bonus/validation.
Which other corrections would you apply ?

Thanks and have a nice day !

[nrol]

Ok so I managed to use the AngularRaDecBuilder and AngularAzElBuilder but I still have 2 problems :

1. The timestep of the generator is fixed while I would like to generate on specific dates corresponding to my real measurements. Is it possible? The eventbasedScheduler asks specifically not to use a DateDetector, for example…
2. The RaDec angles I get are completly wrong while the AzEl are fine and I cannot find where I went wrong in the code. It might be a wrong choice of frame but I don’t see how/where… Any idea ?

The current generation code :

// Ground Station
final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
		Constants.IERS2010_EARTH_FLATTENING,
		FramesFactory.getITRF(ITRFVersion.ITRF_2020, IERSConventions.IERS_2010, false));
final TopocentricFrame topo = new TopocentricFrame(earth, coordStation, "topo");
final GroundStation station = new GroundStation(topo);

// Measurement configuration
double[] raDecError = new double[] {0.0, 0.0};
double[] raDecBaseWeight = new double[] {1.0, 1.0};
double[] azElError = new double[] {0.0, 0.0};
double[] azElBaseWeight = new double[] {1.0, 1.0};

// Time step and dates for the measurements
// TODO : find a way to get the angles only at the dates from the list and not at every step...
double step = 1.0;
final AbsoluteDate t0 = listOfDates.get(0);
final AbsoluteDate t1 = listOfDates.get(listOfDates.size() -1);

// Observable satellite
ObservableSatellite satellite = new ObservableSatellite(0);

// Creation of the angle builders
AngularRaDecBuilder radecBuilder = new AngularRaDecBuilder(null, station, FramesFactory.getICRF(), raDecError, raDecBaseWeight, satellite);
AngularAzElBuilder azelBuilder = new AngularAzElBuilder(null, station, azElError, azElBaseWeight, satellite);

// Generator
Generator generator = new Generator();
generator.addPropagator(propagator);
generator.addScheduler(new ContinuousScheduler<>(radecBuilder, new FixedStepSelector(step, null)));
generator.addScheduler(new ContinuousScheduler<>(azelBuilder, new FixedStepSelector(step, null)));

// Measurement generation
final SortedSet<ObservedMeasurement<?>> measurements = generator.generate(t0, t1);

// Mesures by type
final TreeSet<ObservedMeasurement<?>> radecSet = new TreeSet<ObservedMeasurement<?>>(measurements.stream().filter(meas -> meas instanceof AngularRaDec).collect(Collectors.toSet()));
final TreeSet<ObservedMeasurement<?>> azelSet = new TreeSet<ObservedMeasurement<?>>(measurements.stream().filter(meas -> meas instanceof AngularAzEl).collect(Collectors.toSet()));

// Iterators over the treesets
Iterator<ObservedMeasurement<?>> radecIterator = radecSet.iterator();
Iterator<ObservedMeasurement<?>> azelIterator = azelSet.iterator();

// Refraction model
EarthITU453AtmosphereRefraction refractionModel = new EarthITU453AtmosphereRefraction(coordStation.getAltitude());

And to get the values, I use the following code :

// Ecriture des données
while (radecIterator.hasNext() && azelIterator.hasNext()) {
	
	// Current data
	AngularRaDec radec = (AngularRaDec) radecIterator.next();
	AngularAzEl azel = (AngularAzEl) azelIterator.next();
	
	// Refraction angle
	double refraction = refractionModel.getRefraction(azel.getObservedValue()[1]);
	
	// Angle values
	double ra = (Math.toDegrees(radec.getObservedValue()[0]) + 360.0) % 360.0;
	double dec = (Math.toDegrees(radec.getObservedValue()[1]) + 360.0) % 360.0;
	double az = (Math.toDegrees(azel.getObservedValue()[0]) + 360.0) % 360.0;
	double el = (Math.toDegrees(azel.getObservedValue()[1]) + 360.0) % 360.0;
	double apparent_el = (Math.toDegrees(azel.getObservedValue()[1] + refraction) + 360.0) % 360.0;
	

Do you see anything wrong in my code ?

Thank you in advance !

[nrol]

I found a lead on the 2nd problem : if I change the FramesFactory.getICRF() to a FramesFactory.getGCRF() in the AngularRaDecBuilder I get much better RADEC measurements.
In order to be sure about what the different frames represent, is it true to say:

• ICRF is the sun-centered inertial frame for ICRS
• GCRF is the equivalent but Earth-centered inertial frame for ICRS
• EME2000 is the “old paradigm” frame equivalent and close to GCRF

Is the difference between ICRF and GCRF only the point of origin ? Reading on the Orekit page about frames, it is said:

All celestial bodies are linked to their own body-centered inertial frame, just as the Earth is linked to EME2000 and GCRF. Since Orekit provides implementations of the main solar system celestial bodies, it also provides body-centered frames for these bodies, one inertially oriented and one body oriented. The orientations of these frames are compliant with IAU poles and prime meridians definitions. The predefined frames are the Sun, the Moon, the eight planets and the Pluto dwarf planet. In addition to these real bodies, two points are supported for convenience as if they were real bodies: the solar system barycenter and the Earth-Moon barycenter ; in these cases, the associated frames are aligned with EME2000. One important case is the solar system barycenter, as its associated frame is the ICRF.

Which seems strange, ICRF feels more like a new paradigm frame the the classical one, no?

Anyway, I’m now looking for a way to solve the 1st problem (generating measurements on a list of specific dates)!

[luc]

@nrol

ICRF is the sun-centered inertial frame for ICRS

Yes

@nrol

GCRF is the equivalent but Earth-centered inertial frame for ICRS

Yes

@nrol

EME2000 is the “old paradigm” frame equivalent and close to GCRF

Yes

@nrol

Is the difference between ICRF and GCRF only the point of origin ?

Yes

@nrol

the associated frames are aligned with EME2000.

Aaargggh. It seems to be an error to me. They are aligned with ICRF (and hence GCRF).

Could you open an issue so we don’t forget to fix the documentation?

[MaximeJ]

Hi @nrol,

@nrol

I found a lead on the 2nd problem : if I change the FramesFactory.getICRF() to a FramesFactory.getGCRF() in the AngularRaDecBuilder I get much better RADEC measurements.

I think it’s a bug in AngularRaDec. Line 171 should probably use transformVector instead of transformPosition (the first one apply only the rotation while the second applies rotation + translation).
Could you please open a bug for this one too ? (two bugs in a post !! Thanks)

@nrol

Anyway, I’m now looking for a way to solve the 1st problem (generating measurements on a list of specific dates)!

Two solutions in my opinion:

1. You already have the reference measurements as Orekit ObservedMeasurements
   → Use a ResidualHandler with your propagator (see this post for an example)
   Afterwards you get a Map of (Observed, Estimated) and you can compare them.

2. You don’t have a list or it’s too difficult to obtain.
   → Use the list of dates in the handler instead of a list of ObservedMeasurement
   In the ResidualHandler before line 90 you need to initialize the observed variable.
   Initialize it like you would do a measurement (with sat, frame, weights, sigmas etc.) and put [0,0] for the observed values.
   When observed.estimate(…) is called, Orekit will compute the estimated (i.e. theoretical) measurement. However the residuals in this case will be meaningless.

[nrol]

Thanks for your answers!
@luc, the issue is open.
and @MaximeJ, the other issue is also open.

@MaximeJ

Two solutions in my opinion:

1. You already have the reference measurements as Orekit ObservedMeasurements
   → Use a ResidualHandler with your propagator (see this post for an example)
   Afterwards you get a Map of (Observed, Estimated) and you can compare them.

2. You don’t have a list or it’s too difficult to obtain.
   → Use the list of dates in the handler instead of a list of ObservedMeasurement
   In the ResidualHandler before line 90 you need to initialize the observed variable.
   Initialize it like you would do a measurement (with sat, frame, weights, sigmas etc.) and put [0,0] for the observed values.
   When observed.estimate(…) is called, Orekit will compute the estimated (i.e. theoretical) measurement. However the residuals in this case will be meaningless.

I’m building my listOfDates by reading a TDM file (in raw format) so it should be possible to build directly the ObservedMeasurements. Is there a way to do it easily ? For now, I have a TDM variable from which I get the dates by checking the epoch of the different observations in tdm.getSegments().get(0).getData().getObservations() but since I have 3 measures for each date (RA, DEC, Magnitude), I have to check everytime that the date is not already in my list so it’s not efficient.

2nd question regarding your ResidualsHandler example: how would I use it in my case ? As I understand, it implements OrekitStepHandler so I should probably add it to the propagator via propagator.setStepHandler(residualsHandler) but then, how will it work in order to get the map of evaluations?

(I’m a bit lost on how the handlers work, sorry )

[MaximeJ]

Thank you for opening the issues @nrol !

@nrol

2nd question regarding your ResidualsHandler example: how would I use it in my case ? As I understand, it implements OrekitStepHandler so I should probably add it to the propagator via propagator.setStepHandler(residualsHandler) but then, how will it work in order to get the map of evaluations?

Exactly, you initialize the handler with the ObservedMeasurement list.
ResidualsHandler myHandler = new ResidualsHandler(observedMeasurements)
Add it to the propagator:
propagator.setStepHandler(myHandler)
Propagate
propagator.propagate(...)
Then afterwards, retrieve the map that was filled during propagation:
myHandler.getEvaluations()

As for reading the TDM and creating the list of observations… I’m not super familiar with the "new " CCSDS API so I can’t answer to your question very quickly.

[nrol]

Ok, that’s exactly what I did but I have an error about the prime-meridian-offset.

// Residuals Handler
ResidualsHandler residualsHandler = new ResidualsHandler(tdmMeasurements);
propagator.setStepHandler(residualsHandler);
SpacecraftState state = propagator.propagate(t0, t1);
Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> measurementMap = residualsHandler.getEvaluations();

The error :

date de référence non initialisée pour le paramètre prime-meridian-offset

It might be due to how I created the ObservedMeasurements or the GroundStation, maybe ?

public static List<ObservedMeasurement<?>> getObservedMeasurementsFromTdm(String tdmfilepath, GeodeticPoint coordStation, int numberOfMeasurementTypes) {
	// TDM Parser
	TdmParser parser = new ParserBuilder().buildTdmParser();

	// Read the file
	Tdm tdm = parser.parseMessage(new DataSource(tdmfilepath));
	
	// List of observations
	List<Observation> observations = tdm.getSegments().get(0).getData().getObservations();

	// Measurement list init
	List<ObservedMeasurement<?>> measurements = new ArrayList<>();
	
	// Ground Station
	final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
			Constants.IERS2010_EARTH_FLATTENING,
			FramesFactory.getITRF(ITRFVersion.ITRF_2020, IERSConventions.IERS_2010, false));
	final TopocentricFrame topo = new TopocentricFrame(earth, coordStation, "topo");
	final GroundStation station = new GroundStation(topo);
	
	// Observable satellite
	ObservableSatellite satellite = new ObservableSatellite(0);

	// Loop in the observations
	int ind = 0;
	while (ind < observations.size()-2) {
		if (tdm.getSegments().get(0).getMetadata().getAngleType().equals(AngleType.RADEC)) {
			measurements.add(new AngularRaDec(station, FramesFactory.getGCRF(), 
					observations.get(ind).getEpoch(), new double[] {observations.get(ind).getMeasurement(), observations.get(ind+1).getMeasurement()}, 
					new double[] {0.0, 0.0}, new double[] {1.0, 1.0}, satellite));
		} else {
			measurements.add(new AngularAzEl(station, 
					observations.get(ind).getEpoch(), new double[] {observations.get(ind).getMeasurement(), observations.get(ind+1).getMeasurement()}, 
					new double[] {0.0, 0.0}, new double[] {1.0, 1.0}, satellite));
		}
		ind += numberOfMeasurementTypes;
	}

	return measurements;
}

[MaximeJ]

@nrol

The error :

date de référence non initialisée pour le paramètre prime-meridian-offset

It might be due to how I created the ObservedMeasurements or the GroundStation, maybe ?

Ah yes… This one is a classic (and a pain). It’s because you’re using OD related features (measurement estimation) even though you’re not doing an OD.
Try

station.getPrimeMeridianOffsetDriver().setReferenceDate(anyDate);
station.getPolarOffsetXDriver().setReferenceDate(anyDate);
station.getPolarOffsetYDriver().setReferenceDate(anyDate);

After ground station creation.
You can put any dummy date here since you’re not estimating the EOP.

Thus should fix your issue.

[nrol]

It did! As I saw in another topic here (another of your answer, by the way), I added the same thing as in the TDOAMeasurementCreator :

// Drivers init
for (ParameterDriver driver : Arrays.asList(station.getClockOffsetDriver(),
		station.getEastOffsetDriver(),
		station.getNorthOffsetDriver(),
		station.getZenithOffsetDriver(),
		station.getPrimeMeridianOffsetDriver(),
		station.getPrimeMeridianDriftDriver(),
		station.getPolarOffsetXDriver(),
		station.getPolarDriftXDriver(),
		station.getPolarOffsetYDriver(),
		station.getPolarDriftYDriver())) {
	if (driver.getReferenceDate() == null) {
		driver.setReferenceDate(observations.get(0).getEpoch());
	}
}

It might be an overkill but it seems to work !
Now I’ll go compare the new results with my measurements to see if the data is better but it already seems so !

Thanks a lot, all of you !