Satellite Prediction and Tracking

Orekit usage
1

Hi I am new to OREKIT and learning how to do the calculation correctly. My application is about making prediction of a known satellite using up-to-date TLE over a period of time.

I think I get most of the stuff accurately, except for the part to calculate the Ra Dec of the bird which is complicated. I am seeking an advice if there is a better/more efficient way to do it. Also I found my Ra Dec results are a bit off comparing to the referenced data on https://www.n2yo.com/ website.

I put the test result and code for reference. Any feedback would be really appreciated.

Cheers,
Tony

2024-01-04 19:06:28.373: Starting @UTC 2024-01-04T03:38:21.000Z
2024-01-04 19:06:28.400: Inertial Frame
2024-01-04 19:06:35.544: Earth Frame, time to load: 7144
2024-01-04 19:06:35.545: Earth Body Shape
2024-01-04 19:06:35.546: Ground Station
2024-01-04 19:06:35.550: Station Frame
2024-01-04 19:06:35.555: TLE Setup
2024-01-04 19:06:35.597: TLE Propagation
2024-01-04 19:06:35.597: TLE, time to load: 47
2024-01-04 19:06:35.597: Ref Orbit 2024-01-03T14:02:22.564896Z
2024-01-04 19:06:35.643: lat: 42.301535081271325 lon: 64.62980739846931 alt: 633.8707207112983 azi: 354.22590740439483 ele: -43.20642383798637 ra: 218.55983672907203 dec: 1.0085541200270927 @2024-01-04T03:38:21.000Z
2024-01-04 19:06:35.645: lat: 42.23994306030694 lon: 64.63212753393859 alt: 633.8583575165108 azi: 354.21836707589466 ele: -43.239736561193666 ra: 218.56710804737838 dec: 0.9749277580124084 @2024-01-04T03:38:22.000Z
2024-01-04 19:06:35.646: lat: 42.178349963844255 lon: 64.63443506107998 alt: 633.8460054357508 azi: 354.21083219693526 ele: -43.27304690317423 ra: 218.5743694801105 dec: 0.941303988479511 @2024-01-04T03:38:23.000Z
2024-01-04 19:06:35.647: TLE, time to calculate: 50
2024-01-04 19:06:35.647: Completed

N2YO Results

{"info":{"satname":"IRIDIUM 179","satid":56730,"transactionscount":2},"positions":[{"satlatitude":42.23989589,"satlongitude":64.63640979,"sataltitude":633.86,"azimuth":354.22,"elevation":-43.24,"ra":218.55996888,"dec":1.00631395,"timestamp":1704339501,"eclipsed":false},{"satlatitude":42.17830313,"satlongitude":64.6387172,"sataltitude":633.85,"azimuth":354.21,"elevation":-43.27,"ra":218.56723853,"dec":0.97268697,"timestamp":1704339502,"eclipsed":false}]}

Here is what I have in the code

		long start, end;
				
		// Initial date - 2024-01-04T03:38:21.000Z
		AbsoluteDate initialDate = new AbsoluteDate(2024, 1, 4, 3, 38, 21, TimeScalesFactory.getUTC()); 

		log("Starting @UTC " + initialDate.toString());		

		//AbsoluteDate initialDate = nowToAbsoluteDate();
		AbsoluteDate finalDate = initialDate.shiftedBy(2.0);
				
		// Inertial frame - Earth Mean Equator J2000
		// https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Tutorials/pdf/individual_docs/17_frames_and_coordinate_systems.pdf
		Frame inertialFrame = FramesFactory.getEME2000();

		log("Inertial Frame");

		start = System.currentTimeMillis();
		
		// Earth frame - International Terrestrial Reference Frame (ITRF)
		FactoryManagedFrame earthFrame  = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
		
		end = System.currentTimeMillis();

		log("Earth Frame, time to load: %s", end-start);
		
		// Earth - from ITRF
		BodyShape earth = new OneAxisEllipsoid(
				Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
				Constants.WGS84_EARTH_FLATTENING, 
				earthFrame);
		/**/

		log("Earth Body Shape");
		
		// ground station
		GeodeticPoint groundStation = new GeodeticPoint(
				FastMath.toRadians(45.743397777778), 
				FastMath.toRadians(-123.18619086111), 0);

		log("Ground Station");
		
		// station frame
		TopocentricFrame groundStationFrame = new TopocentricFrame(earth, groundStation, "Ground Station Topo Frame");
		
		log("Station Frame");
		
		// satellite frame
		// Frame satelliteFrame = FramesFactory.getTEME();
		
		start = System.currentTimeMillis();
			
		TLE inputTle = new TLE(
				//"1 56730U 23068W   23365.61045716  .00001544  00000-0  19991-3 0  9995",
	            //"2 56730  86.6779  46.9572 0002221  94.9932 265.1544 14.80201516 33441"
				"1 56730U 23068W   24003.58498339  .00001299  00000-0  16768-3 0  9990",
				"2 56730  86.6779  45.7229 0002244  95.2820 264.8658 14.80209222 33884"
				);

		log("TLE Setup");
		
		// Retrieve the correct propagator that was used to generate the TLE (SGP4 or SDP4)
	    Propagator tlePropagator = TLEPropagator.selectExtrapolator(inputTle);

	    log("TLE Propagation");
		
		end = System.currentTimeMillis();

		log("TLE, time to load: %s", end-start);
	    
	    // Initial reference Orbit for the builder
	    // Orbit referenceOrbit = tlePropagator.getInitialState().getOrbit();

	    log("Ref Orbit " + inputTle.getDate().toString());
	    
	    AbsoluteDate desiredEpoch = initialDate;
	    
	    start = System.currentTimeMillis();
	    
	    Generator generator = new Generator();
	    ObservableSatellite satellite = generator.addPropagator(tlePropagator);	    
        AngularRaDecBuilder builder = new AngularRaDecBuilder(
        		null, 
        		new GroundStation(groundStationFrame), 
        		inertialFrame,
                new double[]{0., 0.}, 
                new double[]{1., 1.}, satellite);
	    
        builder.init(desiredEpoch, desiredEpoch);       
	    
	    while (desiredEpoch.compareTo(finalDate) <= 0.0)
	    {	    	
	    	// Satellite state at a desired epoch
	    	SpacecraftState satelliteState = tlePropagator.propagate(desiredEpoch);
	    	TimeStampedPVCoordinates pv = satelliteState.getPVCoordinates(earthFrame);
	    	Vector3D satPosition = pv.getPosition();	    	
	    	    	
	    	// Satellite position relative to Earth
	    	GeodeticPoint point = earth.transform(satPosition, earthFrame, desiredEpoch);	    	
	    	double lat = FastMath.toDegrees(point.getLatitude());
	    	double lon = FastMath.toDegrees(point.getLongitude());
	    	double alt = point.getAltitude();
	    	   	
	    	// Satellite position relative to the ground station
	    	double azi = FastMath.toDegrees(groundStationFrame.getAzimuth(satPosition, earthFrame, desiredEpoch));
	    	double ele = FastMath.toDegrees(groundStationFrame.getElevation(satPosition, earthFrame, desiredEpoch));
	        
	    	// Satellite Ra Dec calculation
	        double[] baseRaDec = builder.build(desiredEpoch, Collections.singletonMap(satellite, createFakeInterpolator(satelliteState))).getObservedValue();

	        final StaticTransform transform = inertialFrame.getStaticTransformTo(inertialFrame, desiredEpoch);
	        final Vector3D transformedLoS = transform.transformVector(new Vector3D(baseRaDec[0], baseRaDec[1]));
	        double[] rd = new double[] {
	        		FastMath.toDegrees(MathUtils.normalizeAngle(transformedLoS.getAlpha(), FastMath.PI)),
	                FastMath.toDegrees(transformedLoS.getDelta())};
	    	
	    	double ras = rd[0];
	    	double dec = rd[1];	    	       
	    	
	    	log("lat: %s lon: %s alt: %s azi: %s ele: %s ra: %s dec: %s @%s",
	    			lat,
	    			lon,
	    			alt / 1000,
	    			azi,
	    			ele,
	    			ras,
	    			dec,
	    			desiredEpoch
	    			);
	    	desiredEpoch = desiredEpoch.shiftedBy(1.0);
	    }	    	
		
		end = System.currentTimeMillis();

		log("TLE, time to calculate: %s", end-start);

	    log("Completed");
	    
	}
	
	public static OrekitStepInterpolator createFakeInterpolator(SpacecraftState spacecraftState) {
        return new OrekitStepInterpolator() {
            public OrekitStepInterpolator restrictStep(SpacecraftState newPreviousState, SpacecraftState newCurrentState) { return null; }
            public boolean isPreviousStateInterpolated() { return false; }
            public boolean isForward() { return true; }
            public boolean isCurrentStateInterpolated() { return false; }
            public SpacecraftState getPreviousState() { return spacecraftState; }
            public SpacecraftState getInterpolatedState(AbsoluteDate date) { return spacecraftState; }
            public SpacecraftState getCurrentState() { return spacecraftState; }
        };		
	}
2

Hello,

A few remarks/questions already:

  • what do you mean exactly by off w.r.t. n2yo?
  • are we sure they use altitude of 0?
  • be aware that as it stands you have your RADEC with light time delay (as they come from Orekit measurement model) whilst your AZEL does not (you use a pure geometric conversion from position)
  • the part where you use the StaticTransform is useless since initial and final frames are identical.

Cheers,
Romain

3

Hi Romain,

Thanks for the reply.

  • RA DEC off w.r.t n2yo is because I use their API to poll for the info using the same coords with altitude 0. I feed both systems with this station coords:
    Lat: 45.743397777778
    Lon: -123.18619086111
    Alt: 0

My results

ra: 218.55983672907203 dec: 1.0085541200270927

vs n2yo results

“ra”:218.55996888, “dec”:1.00631395

  • I suspect the distance from the station to the satellite must be taken into account as you pointed out instead of a pure geometric conversion.

  • About the static transform I will remove that part.

Regards,
Tony

4

Hi there,

The elevation is negative, that is the satellite is below horizon. Can that be the reason?

5

The elevation is negative, that is the satellite is below horizon. Can that be the reason?

Thats not the question.
My question was why there was such an error between orekit vs. n2yo
ra: 218.55983672907203 vs 218.55996888
I also notice that the result from n2yo seems to be one second ahead of mine which seems to indicate that they probably try to compensate the round trip time since the data returned from them is via a http request over the Internet. They only allow to predict the satellite from the current time plus the extra seconds via their API

https://www.n2yo.com/api/

Get satellite positions
Retrieve the future positions of any satellite as groundtrack (latitude, longitude) to
display orbits on maps. Also return the satellite’s azimuth and elevation
with respect to the observer location.
Each element in the response array is one second of calculation.
First element is calculated for current UTC >time.

I was expecting to see if someone knows how accurate orekit prediction is comparing to some other known systems.

6

Hi @pdxtigre,

If you want to take into account light time delay in your RADEC builder you should probably replace:

public SpacecraftState getInterpolatedState(AbsoluteDate date) { return spacecraftState; }

With

public SpacecraftState getInterpolatedState(AbsoluteDate date) { return spacecraftState.shiftedBy(date.durationFrom(spacecraftState.getDate())); }

In method createFakeInterpolator (see this discussion).

Regards,
Maxime