Help with IODGooding

I have three observations, instantiated as AngularRaDec objects in the topocentric J2000 frame (all from the same site for initial testing, but eventually I will need to support different ground-stations across observations). I am trying to use IODGooding for an initial IOD of a track in GEO, but I’m getting clearly wrong values back (e.g. negative semi-major axes with eccentricities > 100). I think it might be an issue with my frame or using the zenith vectors for the observer locations, but I’m hoping you can help. I tried a simpler example with the observers at the center of Earth like the testing code, and that worked for me.

I’m using the python wrapper, but I think this is just an issue with how I’m using the method rather than anything python specific. Code below. I’ve omitted the imports/constants for brevity, but can include if you think that might be the source of the issue.

mu = Constants.IERS2010_EARTH_MU
IOD = IodGooding(eme, mu)

Transform observer positions from J2000 topocentric to EME

t1 = obs1.getStation().getBaseFrame().getTransformTo(eme, obs1.getDate())
obs_pos1 = t1.transformVector(obs1.getStation().getBaseFrame().zenith)
t2 = obs2.getStation().getBaseFrame().getTransformTo(eme, obs2.getDate())
obs_pos2 = t2.transformVector(obs2.getStation().getBaseFrame().zenith)
t3 = obs3.getStation().getBaseFrame().getTransformTo(eme, obs3.getDate())
obs_pos3 = t3.transformVector(obs3.getStation().getBaseFrame().zenith)

Transform observation RaDecs to unit vectors

x1 = math.cos(obs1.getObservedValue()[0])*math.cos(obs1.getObservedValue()[1])
y1 = math.sin(obs1.getObservedValue()[0])*math.cos(obs1.getObservedValue()[1])
z1 = math.sin(obs1.getObservedValue()[1])
los_1 = Vector3D(x1, y1, z1)

x2 = math.cos(obs2.getObservedValue()[0])*math.cos(obs2.getObservedValue()[1])
y2 = math.sin(obs2.getObservedValue()[0])*math.cos(obs2.getObservedValue()[1])
z2 = math.sin(obs2.getObservedValue()[1])
los_2 = Vector3D(x2, y2, z2)

x3 = math.cos(obs3.getObservedValue()[0])*math.cos(obs3.getObservedValue()[1])
y3 = math.sin(obs3.getObservedValue()[0])*math.cos(obs3.getObservedValue()[1])
z3 = math.sin(obs3.getObservedValue()[1])
los_3 = Vector3D(x3, y3, z3)

date_1 = obs1.getDate()
date_2 = obs2.getDate()

date_3 = obs3.getDate()

geo_dist = 37000000.0 #meters

range_1 = geo_dist
range_3 = geo_dist

estimate = IOD.estimate(obs_pos1, obs_pos2, obs_pos3, los_1, date_1, los_2, date_2, los_3, date_3, range_1, range_3)

Could someone help point me to what I’m doing wrong?

Hi @mlifson

Welcome to the Orekit forum.

I don’t see an error on your code but I can give you some suggestions:

  • You can use the constructor Vector3D(alpha, delta) (where alpha is getObservedValue()[0] and delta is getObservedValue()[1]) to build the satellite position:

       position1 = Vector3D(obs1.getObservedValue()[0], obs1.getObservedValue()[1])
    
  • Make sure your angles are in radians.

  • Compute the los vector using:

       los_1 = position_1.normalize()
    
  • Compute range_1 value using:

       range_1 = position_1.getNorm()
    
  • When calling the estimate method be sure to convert all the vectors to the same inertial frame you used to create the class

I am very far to be an expert of IOD methods. Maybe someone can give you better recommendations. My suggestions are inspired by the Orekit test class of IOD Gooding method.

If my suggestions do not help you. Can you provide us a working test in order to find the issue?

Regards,
Bryan

Hi @mlifson,

Welcome to the Orekit forum!

I think there may be a mistake in your code when you are transforming the zenith vector into EME frame.
I would transform the underlying geodetic point of the ground station instead of its zenith vector.
For this, try replacing:

t1 = obs1.getStation().getBaseFrame().getTransformTo(eme, obs1.getDate())
obs_pos1 = t1.transformVector(obs1.getStation().getBaseFrame().zenith)

With:

obs_pos1 = obs1.getStation().getBaseFrame().getPVCoordinates(obs1.getDate(), eme).getPosition()

Like Bryan, I’m not an expert of IOD methods.
But I hope our comments will help you debugging your code.

Maxime

Hello,

Thank you both for your suggestions. I believe I have implemented both of them, but the code still does not generate a reasonable OD. I have attached a version of the code and a sample data set. I don’t have authorization to share my actual data set, so I have used another satellite/ground site as test data. This ground site/satellite pair errors and won’t run the OD, but I think the problem is likely similar.

Access Data MIT.csv (161.4 KB) OD_script_sample.py (5.2 KB)

I appreciate your help!

Miles

Hi,

  1. The latitude and longitude arguments are swapped in the GeodeticPoint constructor, latitude should be first (it happens to me quite often too):

     sitella = GeodeticPoint(site_lat, site_lon, site_elev)
    
  2. You used the IODGooding::estimate method without the nRev and direction arguments, therefore there must be no more than half a revolution between start and final date (https://www.orekit.org/site-orekit-development/apidocs/org/orekit/estimation/iod/IodGooding.html). As a dirty fix, I hardcoded some values which correspond to <AbsoluteDate: 2014-08-27T23:58:52.816>, <AbsoluteDate: 2014-08-28T04:08:52.816> and <AbsoluteDate: 2014-08-28T08:18:52.816>:

    obs_1_index = 0
    obs_2_index = 250
    obs_3_index = 500
    
  3. The variables range_1 and range_3 equal 1.0, because there is no range information in the vectors position_1 and position_3. These variables should instead be:

     range_1 = float(1e3*AER_data.iloc[obs_1_index]["Range (km)"])
     range_3 = float(1e3*AER_data.iloc[obs_3_index]["Range (km)"])
    

With these three changes, the estimation gets a little bit better :wink:

Esimate: Keplerian parameters: {a: 4.3370460159959204E7; e: 0.06512061963690145; i: 3.145738954932238; pa: 52.42655895744585; raan: -131.86014692789726; v: 9.455609345584236;}
Truth: Keplerian parameters: {a: 4.2165631616E7; e: 0.001117; i: 4.388; pa: 185.542; raan: 324.234; v: 190.642;}
IOD Error:
 	a: -1204828.5439592078 
	e: -0.06400361963690146 
	i: 0.02168154540569811 
	pa: 2.3233027314369634 
	raan: 7.960344562966587 
	nu: 3.1622990767242385

Cheers
Clément

Clément,

Thank you for pointing out these errors. I am able to replicate your results for the MIT ground-station after making these corrections.

Hours between first and last observation: 8.333333279995713
Esimate: Keplerian parameters: {a: 4.3370460159959204E7; e: 0.06512061963690145; i: 3.145738954932238; pa: 52.42655895744585; raan: -131.86014692789726; v: 9.455609345584236;}
Truth: Keplerian parameters: {a: 4.2165631616E7; e: 0.001117; i: 4.388; pa: 185.542; raan: 324.234; v: 190.642;}
IOD Error:
a: -1204828.5439592078
e: -0.06400361963690146
i: 0.02168154540569811
pa: 2.3233027314369634
raan: 7.960344562966587
nu: 3.1622990767242385

However, I still get bad results for other ground-station locations. See below for a ground-station in Phoenix, AZ generated with the same method as the MIT site. I also get issues for San Francisco and my actual data.

Hours between first and last observation: 8.333333279995713
Esimate: Keplerian parameters: {a: 3.2325328703900866E7; e: 0.4903417768138178; i: 36.60730423351098; pa: -91.51426852343015; raan: -96.61311701825664; v: 179.2653413922996;}
Truth: Keplerian parameters: {a: 4.2165631616E7; e: 0.001117; i: 4.388; pa: 185.542; raan: 324.234; v: 190.642;}
IOD Error:
a: 9840302.91209913
e: -0.4892247768138178
i: -0.5623329415765146
pa: 4.835544099023385
raan: 7.345167839494438
nu: 0.19856015057972565

OD_script_sample.py (5.9 KB) Access Data Phoenix.csv (161.4 KB)

Is there anything else I am doing wrong?

You have a typo in the Phoenix station longitude: site_lon = float(-pi/112.076) # rad should be site_lon = float(-pi/180*112.076) # rad. Using the np.deg2rad function is recommended though. And saving the station coordinates in a file avoids the risk of a mistake when commenting/uncommenting station coordinates in your script.

Other than that, I don’t see any mistake in the code, but I get the same results as you for the Phoenix access data. Either there something wrong in this data, but the azimuth seems to make sense, or something else, but I am far from being an expert in IOD. Is it simulated data or actual flight data?