Earth orbit not quite right :(

[selden]

As long as we're concentrating on orbital accuracy...

It seems that the Earth's orbit definition in Celestia 1.2.5pre7 still has some errors creeping in somewhere.

I was trying to generate a high-quality trajectory for Cassini to watch its encounters with various objects. Using xyz values for Cassini from Horizons with a sample interval of about 15 minutes, Celestia shows the time of closest approach to Earth to be at 1999 08 18 03:38:13 UTC at a distance of 92,684 km over the Arabian penninsula, instead of 3:28 at an altitude of 1,171 km over the eastern South Pacific

I double checked Horizons' trajectory by asking for the elliptical orbit elements for 3:28. Horizons reported the perigee to be 7,549 km, which corresonds to an altitude of 1,171 km.

*sigh*

[chris]

As long as we're concentrating on orbital accuracy...

It seems that the Earth's orbit definition in Celestia 1.2.5pre7 still has some errors creeping in somewhere.

I was trying to generate a high-quality trajectory for Cassini to watch its encounters with various objects. Using xyz values for Cassini from Horizons with a sample interval of about 15 minutes, Celestia shows the time of closest approach to Earth to be at 1999 08 18 03:38:13 UTC at a distance of 92,684 km over the Arabian penninsula, instead of 3:28 at an altitude of 1,171 km over the eastern South Pacific

I double checked Horizons' trajectory by asking for the elliptical orbit elements for 3:28. Horizons reported the perigee to be 7,549 km, which corresonds to an altitude of 1,171 km.

*sigh*

I just compared the Earth's position against values from Horizons. I created an .xyz file for Earth spanning 1 Jan 1998 to 1 Jan 2004 with a time step of one day. The Horizons trajectory matched the VSOP87 orbit used by Celestia to within a few kilometers at the sample points (interpolation error caused up to a 2000 kilometer difference exactly between sample points.) And Cassini never got closer than 90000km to either the VSOP87 Earth or the Horizons Earth. Something else is happening here, and I don't know what . . .

--Chris

[billybob884]

well that can't be too good. You know your in trouble when the maker of the program has no idea whats going on

[selden]

It's Cassini's xyz positions.

I can't tell if it's Horizons returning the "wrong" values or if Celestia is using them wrong, though. I put wrong in quotes because maybe I'm asking for the wrong coordinate system.

When I use Horizon's values for Cassini's sun-centered elliptical elements for the time of closest approach, Celestia puts Cassini at the correct location.

This is really strange, though, since I used the same query methods for the orbital parameters for the asteroid Annefrank and it was quite reasonable: the elliptical and xyz parameters matched exactly at the epoch specified.

[selden]

Well, I'm sorry to say that it does seem to be something not quite right in Celestia, as I guess Chris is aware.

It seems suggestive to me that the xyz positional error puts the spacecraft very close to where the orbit paths are drawn. Could the same bug be affecting all three of them? (xyz position, xyz orbit display, elliptical orbit display).

Sorry if I'm beating a dead horse, but I thought a picture of the problem might help.

When I compare Horizon's xyz coordinates for Cassini and the Earth for the epoch of the flyby, the distance between them is quite reasonable: about 3,500 km. At first this surprised me: "What? that's undergound!" until I realized that the distance is to the barycenter of the Earth-Moon system, not to the center of the Earth. Cassini's perigee was between the Earth and the Moon.

Anyhow, below s a thumbnail linking to a screenshot by Celestia illustrating the problem. Cassini's position when using elliptical orbital elements (cassini-es) is on one side of the planet while the xyz coordinate is 100,000 km away on the other side (cassini-xyz).

In the picture below, he horizontal blue line is the orbit path drawn for Cassini's elliptical orbit, the red line is the orbit path drawn for its xyz trajectory, and the diagonal blue line is the orbit path drawn for the earth.
Note that cassini-es ia above and very close to the crescent earth, while cassini-xyz is below it.

[url=http://www.lns.cornell.edu/~seb/celestia/cassini-orbits.jpg]

[/url]

[selden]

This may not be the problem ,but...

While looking through the code, I came across these lines in samporbit.cpp:

Code: Select all

struct Sample
{
    double t;
    float x, y, z;
};


and later in the routine SampledOrbit::addSample

Code: Select all

    Sample samp;
    samp.x = (float) x;
    samp.y = (float) y;
    samp.z = (float) z;
    samp.t = t;


Everywhere else the position coordinates are doubles.

I really don't know C++ (I still program in Fortran 77), but it looks to me like in some places the routine(?) Point3D is explicitly being passed doubles, but in others it's being passed elements of this single-precision structure.

Could the position difficulties be caused by this loss of precision?

[chris]

This may not be the problem ,but...

While looking through the code, I came across these lines in samporbit.cpp:

Code: Select all

struct Sample
{
    double t;
    float x, y, z;
};


Could the position difficulties be caused by this loss of precision?

The use of floats instead of doubles in sampled trajectories does limit the precision. However, the position of Cassini near Earth should be off by no more than about 12km (mean distance of Earth from Sun / 2^23, where 23 is the number of bits in the mantissa of a single precision float.) There will eventually be several types of sampled trajectory, some with single precision and other with double precision samples. It's just an accuracy vs memory usage tradeoff . . .

--Chris

[selden]

Looking through some other modules, I see some other numeric constructs that look somewhat dubious. Like I said, I'm not familiar with how the C++ compiler handles them, but the Fortran compilers I'm familiar with would not do what you might expect.

For example, in orbit.cpp I found the lines

Code: Select all

    double dt = 1.0 / 1440.0; // 1 minute
    Point3d p0 = orbit->positionAtTime(t0);
    Point3d p1 = orbit->positionAtTime(t1);
    Vec3d v0 = (orbit->positionAtTime(t0 + dt) - p0) / (86400 * dt);
    Vec3d v1 = (orbit->positionAtTime(t1 + dt) - p1) / (86400 * dt);


I have to be concerned that the operation "1.0/1440.0" would be performed in single precision by the compiler and the resulting single precision constant would then be promoted to double precision. I think you really want this calculation to be performed in double precision.

I'll admit that I'd like Fridger's comments on this. I'm sure he's had to worry about such things a lot more than I have.

[chris]

I have to be concerned that the operation "1.0/1440.0" would be performed in single precision by the compiler and the resulting single precision constant would then be promoted to double precision. I think you really want this calculation to be performed in double precision.

C++ is very clear that 1.0/1440.0 is a double precision calculation. It would be single precision only if it was written 1.0f/1440.0f. And 1/1440 will produce zero, since without the decimal point, the values are treated as integers. Writing Celestia has taught me to be extremely careful with numeric precision . . . Not that I don't still screw up from time to time

--Chris

[Rassilon]

Thats what happens when you build software...You will always find a faster better method after its coded

I agree that doubles would be good to use throughout the calculations most specifically those RA and Dec coordinates...This would also reduce the gaps/repeating patterns in alot of the nebula/galaxy models inc the star clusters formed with the program I wrote...No biggie I suppose...get the rest fixed first

[selden]

OK.

I have to believe, though, that *somewhere* there's some equivalent assumption about the meaning of a numeric operation is mistaken and that's what's biting us.

For example, could there be some circumstances in which there's some "offset" in the elliptical and high precision orbit calculations so that those results are the ones placing objects in the wrong places while the "nearby" xyz positions are the correct ones?

I've gotten burned often enough when I've assumed "Oh, that [insert your favorite blatently obvious mistake here] couldn't possibly be what's happening" when in fact is was...

[Calculus]

Don't you think this is because the elliptical elements of an orbit have to be calculated with respect to the earth-moon barycenter instead of the earth center ?

[selden]

Paul,

That's something I'm worried about. Unfortunately, that could only account for about a 3,500 km difference.

What we're seeing is that Cassini's xyz coordinates are placing it about 100,000 km away from the position calculated using sun-centered elliptical orbit elements.

Hmmm. What is the actual origin for the xyz coordinate system?

Grant did discover that the rotational elements for the various satellites are referenced to Earth's ecliptic, not to that of the planet they're orbiting. Could an equivalent coordinate transform have been overlooked here?

[selden]

The Earth center/barycenter offset is fine.

I just asked Horizons for the xyz position of earth's barycenter at the time of Cassini's flyby. When I goto that location, Celestia says my altitude is -1,493 km, exactly as it should be.

ARRGH! Why isn't it 100,000 km off???

[selden]

The problem only appears when the xyz trajectory file is "long."

When I provide a very short xyz file --- five lines of jd,x,y,z --- Cassini is placed just about exactly where it should be. I say "just about" because it's still 100 km away from the elliptical position. Hopefully that's just a result of the linear interpolation being used -- Cassini was moving like the proverbial bat-out-of-hell at the time of its Earth flyby.

[p.s. Yup: the minimum distance between the short xyz file and the solar eliptical orbit for Cassini is only 1km. I really can't get excited about that.

So what about long xyz files causes Celestia to get confused???]

[selden]

Problem resolved.

It's the format of the input file.
The problem is not in the position calculations within Celestia.
The problem is in the I/O library that it uses.

The cassini.xyz file that Chris has been providing was written containing xyz values in nnn.nnn floating point format.

The files I've been using for testing this evening were written containing xyz values in exponential format: n.nnnn E+mm

This exponential format works fine: Cassini is placed exactly where it should be.

Floating point format fails: Cassini is placed way out in left field.

A complete Cassini xyz trajectory using exponential format is now available in
http://www.lns.cornell.edu/~seb/celestia/cassini-all.zip
This zip archive also contains an appropriate .ssc file to invoke it
along with a readme.txt file describing the contents.

Additional information is available on my "spacecraft orbits" Web page at
http://www.lns.cornell.edu/~seb/celestia/spacecraft.html

This certainly is a relief.
Hopefully Chris will sleep better tonight

[chris]

The problem only appears when the xyz trajectory file is "long."

When I provide a very short xyz file --- five lines of jd,x,y,z --- Cassini is placed just about exactly where it should be. I say "just about" because it's still 100 km away from the elliptical position. Hopefully that's just a result of the linear interpolation being used -- Cassini was moving like the proverbial bat-out-of-hell at the time of its Earth flyby.

[p.s. Yup: the minimum distance between the short xyz file and the solar eliptical orbit for Cassini is only 1km. I really can't get excited about that.

So what about long xyz files causes Celestia to get confused???]

I'm pretty sure long xyz files are not the problem . . .

I created a new Cassini trajectory that uses a one day sample interval over most of the trajectory, and a 15 minute interval in the several days around the time of the Earth flyby. It works great, with Cassini passing about 1100km above the surface of the Earth (a ~400km discrepancy, but keep reading . . .)

I'm inclined to blame all the problems on interpolation error--Cassini was over one million kilometers from Earth by the end of the day that it passed 1500km over the surface. It's speed relative to the Earth was roughly 17km/s. The one day sample interval just didn't cut it. Even with a 15 minute sample interval, the errors can be considerable, as Cassini traveled about 15000 km in that time.

My new Cassini trajectory is here:

http://www.shatters.net/~claurel/trajectories/cassini.xyz

--Chris

[chris]

Problem resolved.

It's the format of the input file.
The problem is not in the position calculations within Celestia.
The problem is in the I/O library that it uses.

The cassini.xyz file that Chris has been providing was written containing xyz values in nnn.nnn floating point format.

The files I've been using for testing this evening were written containing xyz values in exponential format: n.nnnn E+mm

This exponential format works fine: Cassini is placed exactly where it should be.

Floating point format fails: Cassini is placed way out in left field.

I don't think that this was the problem . . . The trajectory file I linked to above does not use exponential format for its values and it works fine. At least, I hope it works fine . . . could you give it a try Selden?

--Chris

[selden]

Chris,

Sorry: I may be doing something wrong, but I still see Cassini 100,000 km away when using the xyz file in your link above. Your file above is dated October 16th, 2002. That's the one I've been using all along.

I just looked in your trajectories directory: maybe you meant to link to cassini2.xyz?
Your cassini2.xyz does seem OK!
But not quite identical to my exponential format trajectory -- I made a point of filling in more points around the times of the flybys.

So... maybe a trajectory update would fix the NEO trajectories, too!

[t00fri]

Looking through some other modules, I see some other numeric constructs that look somewhat dubious. Like I said, I'm not familiar with how the C++ compiler handles them, but the Fortran compilers I'm familiar with would not do what you might expect.

For example, in orbit.cpp I found the lines

Code: Select all

    double dt = 1.0 / 1440.0; // 1 minute
    Point3d p0 = orbit->positionAtTime(t0);
    Point3d p1 = orbit->positionAtTime(t1);
    Vec3d v0 = (orbit->positionAtTime(t0 + dt) - p0) / (86400 * dt);
    Vec3d v1 = (orbit->positionAtTime(t1 + dt) - p1) / (86400 * dt);


I have to be concerned that the operation "1.0/1440.0" would be performed in single precision by the compiler and the resulting single precision constant would then be promoted to double precision. I think you really want this calculation to be performed in double precision.

I'll admit that I'd like Fridger's comments on this. I'm sure he's had to worry about such things a lot more than I have. :)

Since I was cut off from the forum for >2 days (?), here are somewhat late comments.

I agree with Chris.

In Fortran, 1.0/2.0 would indeed mean a Real*4 division of two float numbers and to do it in Real*8, one has to write 1.D0/2.D0.
C++ has a different syntax, however: 1.0/2.0 is double by default and 1.0f/2.0f is float.

But even in a genuine double calculation, in principle, there may be cases in this game involving very large numbers, where the precision becomes insufficient, e.g. let p1, p2 be two almost equal, but very large positive doubles, then (p1+p2)/(p1-p2) may become insufficiently accurate...

Bye Fridger