Member Discussions

Keeping this possibly short and relatively simple, I'm having a measure of difficulty finding the information I'm looking for with respect to real-life astrodynamics and the equations for which I could use to translate certain physics into a way that both I and my codebase could make sense of.

The foremost problem I'm working on in this respect has to do with calculating orbital paths and moving an object along a path around a designated object or focal point. The primary example would be a planet orbiting a star.

Now, I've done a fair amount of research so far on the matter, and I've found a few things regarding keplar and newtonian laws and whatnot, but the equations I've found so far don't really translate as well into use for a game (so far, according to my own understanding).

What I'm trying to come out of from all this is a number of usable functions.

The first of which deals with moving a planet or other object from one coordinate location to another coordinate location, based on the fixed velocity/speed of that object (I know the velocities change, but for the purposes of the game they're staying constant) along the orbital path, based on a focal coordinate location, usually a star. Whether or not the orbital path shape is circular or eliptical doesn't make a whole lot of difference to me. From what I read, circular tends to be a far cry easier to deal with (since the orbital radius stays the same), so I'll probably lean towards that direction; although eliptical would go for a more realistic feel.

The second of which deals with checking to see if an object is in the orbital path of another object during it's movement; I.e, if the game is moving an asteroid along its orbital path and during the move it's coordinates cross a planet, the game would know whether or not to smack said asteroid into said planet.

All I'm really looking for here is some ideas by which I can be pointed in the right direction to finding a solution to these problems.

In tigermud, I created a coordinates 0-7 around the planet. Then I put in a sun and moon traveling through these coordinates each day. Imagine that coordinate 0 is sunrise, 1-2 is daytime, and 3 is sunset. 4-7 is night time, then recycling to zero at sunrise.

On the NASA Web site, I looked up the moon's orbital time (in seconds) around the Earth and the Earth's rotation (in seconds) and used it in the code.

In the end, I have the moon and sun moving through the 0-7 coordinates together, however since the moon orbit is slightly out of sync with the sun, as days pass the sun and moon get into different coordinates and you go through the lunar phases of full moon, third quarter, half, etc.

To do the phases I just had a simple table that told me if the sun and moon were in the same coordinate, then the moon was a new moon. If they were opposite each other, then it was a full moon, etc.

I first tried like you did with the physics, but in the end, people only see the location and phases of planets in my game. I figured it would be educational process to write all that code, but a waste of CPU and time since nobody would appreciate it.

Cheers!

Well, thanks for the input. That's one possible way of doing things. But unfortunately, it's really not an option I can explore for my game.

More or less, I'm trying to go for a near-full-depth universe system. The space side of my game started out as swr. Despite the inherent buggishness and inefficiency, it provided at least a decent framework for a good space system, in my view. I've reconstructed it twice now, and under the current system of universe I'm almost to where the entire in-game universe has a very explorable feel and depth to it (galaxies, starsystems, planets, atmospheres, etc etc). As such, when I looked at the swr way of assigning a planet a fixed coordinate, from which it never changes, I decided that it would be a great idea to give planets orbital paths so that the game can move planets and other bodies around their focal points (stars, etc) automaticly.

But this prevents the obvious mathematical problems. I won't be following ALL the laws of physics involved in creating this universe (fixed orbital radiuses, constant planetary velocities, circular orbits, etc). I just need to find a way to calculate where to move an orbiting body to, based on the coordinate it started from, the speed its going, and the focal point the object orbits. Think of the function I need as being the code equivalent of a protractor.

Going to be researching this until I find a solution, but in the meantime I would continue to appreciate any input people have on the matter. Again, thanks, AdamMil.

You probably want a tree of reference frames. Each frame would contain a bunch of other frames, and each of those contained frames would have coordinates in the base frame. For example, the top frame could be galaxy, which contains solar system frames, which contain planet frames, which contain moon frames. To compute the galactic coordinate of something you sum up offsets from the center of each frame from that object to the root of the tree. That is, the galactic coordinate offset of a moon is the offset from moon to the origin of the planet frame, plus the offset from the planet to the origin of the solar system frame, etc.

You might run into precision problems if you try to use a single set of galactic coordinates when talking about things at a low level. Thus if you are trying to resolve something going on inside of a spaceship, you might be better off just working in the coordinates relative to the origin of the spaceship's frame rather than going all the way back to galactic coordinates.

As far as converting a circular orbit to cartesian (x,y,z) coordinates, it works like this:

say that the orbital period is T
say that the planet was on the +x axis at time t0
say that the orbital radius is r
say that the current time is t

then if the orbit is in the (x,y) plane,

x = r*cos( 2*Pi*( t - t0 ) / T )
y = r*sin( 2*Pi*( t - t0 ) / T )

This gives you the offset from the center of the solar system frame to the planet. You can use the same scheme for moons. If you want the orbit out of the xy plane, there are formulas for that, but I think your players will thank you if you don't do that.

"You probably want a tree of reference frames. Each frame would contain a bunch of other frames, and each of those contained frames would have coordinates in the base frame. For example, the top frame could be galaxy, which contains solar system frames, which contain planet frames, which contain moon frames. To compute the galactic coordinate of something you sum up offsets from the center of each frame from that object to the root of the tree. That is, the galactic coordinate offset of a moon is the offset from moon to the origin of the planet frame, plus the offset from the planet to the origin of the solar system frame, etc.

You might run into precision problems if you try to use a single set of galactic coordinates when talking about things at a low level. Thus if you are trying to resolve something going on inside of a spaceship, you might be better off just working in the coordinates relative to the origin of the spaceship's frame rather than going all the way back to galactic coordinates."

No, I don't help with that part, sorry. I already have the universe constructed. The part I need help with is the equations to move an object along it's orbital path.

The equations you posted might do that, but I'm not exactly sure, and somewhat doubtful, as they don't really appear to.

So far from doing research on the matter, I keep running across a lot of information I don't need. Most/all of the sites regarding orbits keep throwing out a lot of time this, period that, mass this, acceleration that. That's all well and good for doing calculations in real life, but for the purposes of the game it bears little relevance, if any at all.

Restating what my mission is here in all this, is to take what will be the known information: starting coordinates of an object, the radius from the object to the focal point it orbits, and the speed at which the object orbits the focal point i.e. the near-exact distance its going to travel along its orbital path per related game tick, and turn that into a usable equation by which I can enter any given set of the known information and have it calculate the proper destination coordinates for the object along its path so that it can be moved accordingly.

As it stands, the closest I've come to this is partially what you posted just now (soft of), and some acticles regarding in-game physics on the gamedev.net site. Unfortunately, the one article there that seemed the most promising, was written by someone who decided they didn't need to provide any meaningful explanation to their equations.

In the research I've been doing so far on the matter, it seems that what I'm most closely looking for here is called uniform circular motion. But don't think of it in terms of planets and stars with gravity and friction and whatnot. Think of it in terms of one single 'cartessian' vertex point rotating around a fixed vertex point.

Again, all help is greatly appreciated.

May I suggest using vectors instead of coordinates?

Calculating the change in the x and y is easy, but applying the change to the x and y of a body already in motion means you'll have to calulate the change in the angle of the momental tangent. If you have a lib with a function to convert Cartesian to Radial and back, you can let the machine do your work for you.

Convert the present x,y position to radians. Multiply the speed of the object (stored as degrees) by the time since the last update (another hint to make this event-driven) and then just convert back to plot it.


HTH,
Keri

Well, I've done a little bit of exploration into vectors so far, and I'm a little bit at a loss as to how to make practical use of them for this instance so far.

'If you have a lib with a function to convert Cartesian to Radial and back, you can let the machine do your work for you.

Convert the present x,y position to radians. Multiply the speed of the object (stored as degrees) by the time since the last update (another hint to make this event-driven) and then just convert back to plot it.'

That almost sounds to me like it wants to make sense, but I'm not exactly sure how that would take into account what the x,y position is orbiting, how far away, what direction, etc etc.

And Kjartan, I gave your formulas and descriptions of what they do some further thought, and it's beginning to sound like they might actually do the trick. However, if you could given a little further detail and description it might help clear up some uncertainties about how exactly to use them. I tried doing some test calculations with them, and the results don't even seem remotely correct.

Also, the equations provided don't exactly seem to account for where the actual starting coordinates of the vertex in question is.

For example, I tried a test with the following scenarios, following what you stated as closely as possible. I decided to use a x/y plane for this, completely ignoring the z since that's going to overcomplicate things, and created a circle with a radius of 100, the center of which is the coordinate 0,0. I'm assuming for the purpose of this orbit, the orbit is going counter-clockwise (although I'm both assuming and hoping the equations are direction-independant), starting at the right-side point of the circle where it rests on the x axis, coordinate 100,0. So, moving onward:

'say that the orbital period is T'

From what I gather, the orbital period of the circle means the circumferance, i.e. the length around the circle it's going to take for one revolution. For a circle of radius 100, this came out to be approximately 628 (c = 2*pi*r).

'say that the planet was on the +x axis at time t0'

Now this one in particular is a little ambiguous. The best I can assume of that is that you mean the planet was ON the x axis at time t0, meaning that the planet would have a coordinate of 100(x),0(y) at t0. So for the purpose of this, we're going to assume the time 0 is exactly 0; i.e. the starting point for the orbit.

'say that the orbital radius is r'

Well, we've established the radius is 100, so it is thus.

'say that the current time is t'

For the purpose of this example, I decided to use the point where the circle touches the Y axis at coordinate 0,100 as the place where the current time would have the planet/object that is in orbit end up at. This equates to 1/4th of the orbital period, so 't' would end up being 157.

Ok. Now, for the values so far:

T = 628
t0 = 0
r = 100
t = 157

Entering these values into the equations provided, I end up with:

x = 100*cos( 2*Pi*( 157 - 0 ) / 628 )
y = 100*sin( 2*Pi*( 157 - 0 ) / 628 )

Now, calculating these equations I end up with the following values:

x = 99.96
y = 2.74

Which, of course, doesn't really seem to be all that helpful.

I'm willing to bet there's a pretty good chance I'm completely misunderstanding what you stated, so I'm hoping you'll try giving some further clarification so I can make sense of your equations.

By 'orbital period' I meant the time for one orbit, e.g. one year (for earth) or whatever. Then t is the elapsed time since the planet was on the +x axis. But actually, the way you did it is fine for the sake of this example, let's just say the orbit lasts 628 days and we are current 157 days into it. (Note: T and t have to be in the same units, both in days or both in seconds or whatever.)

The reason you aren't getting the expected answer is that I had assumed you were using the C sin() and cos() functions, which take their argument in radians. You may be using a calculator which takes the argument in degrees. In radians, the equations you have written out give

x = 100 * cos( 2 * Pi * ( 157 - 0 ) / 628 ) = 0
y = 100 * sin( 2 * Pi * ( 157 - 0 ) / 628 ) = 100

In degrees, you need to replace the 2*Pi in both of my formulas with 360 (basically that term is 'one circle', since there are 360 degrees in a circle and 2*Pi radians in a circle).

Ok, let me try to work through the exact problem you have. Let's say the orbit's center is at x=x0,y=y0 (previously I assumed this to be at 0,0), and let's say right now, which I will call time t1, the planet is at x=x1,y=y1. What we don't know is t0 and r. Here is how to compute them.

r is easy: it's sqrt((x0-x1)^2+(y0-y1)^2).

t0 is a little trickier. This parameter tells you where the planet is in its orbit. You may not care about t0; if you just set all the t0's in the whole universe to be zero, you would get a perfectly reasonable and consistent universe. But, each planet would necessarily be on the +x axis (like you thought, I meant 'x = r, y = 0' by that) when the clock was at zero.

If it's currently time t1 and the planet is at (x1,y1) (meaning x = x1, y = y1) then we can compute t0. There's a C function called atan2 you can use for this. Probably whatever language you are using has such a function; I know Java, C, C++, C# all do. Here's a C snippet to compute r and t0 from x0, y0, x1, y1, t1, and T:

r = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
t0 = t1-T*atan2(y1/r,x1/r)/(2.0*PI);

(Comment added by Kjartan on Tue Sep 20 11:31:18 2005)

Whoops, that last line should actually be:

t0 = t1 - T * atan2((y1-y0)/r,(x1-x0)/r)/(2.0*PI);

To reiterate what all of these variables are:
x0, y0 are the coordinates of the center of the orbit
x1, y1 are the coordinates of the planet right now
t1 is the time right now
T is the time it takes to complete one orbit

and the solution t0 is the t0 for use in my earlier formulas, which is the time at which the planet was on the +x axis at x = r, y = 0.

Ahh this is my favorite question! (Heh, mad at myself that I didn't notice the thread until today.)

Anyways, the easiest and best (in my opinion) way to have properly orbiting bodies is to not hard-code in perfectly circular orbits and such, but just let gravity do the work for you. I've got an orbital system running on my MUD based purely off of laws of gravitation / realistic collision equations. Telnet to americanairport.us port 3000 if you'd like to have a look. You can do a 'bodyinfo 4' to look at Earth and 'bodyinfo 1' to look at the Sun, and then 'orbitinfo 4 1' to see how Earth is orbiting the Sun.

Basically I'm just calculating the gravitational influence of all the other bodies on body A, and then figuring out it's acceleration, then updating it's velocity and position based off that acceleration. Then I move onto body B and do the same. It's not taking up as much CPU as you'd think (but load will increase exponentially with more bodies.)

Collisions are easy... Just calculate the distance between the 2 objects sqrt (x^2 + y^2 + z^2) and see if that is less than the radius of the 2 objects added together.

Ahh well, must go to work, I'll explain more later.

Well, Lockheed, I appreciate the input and all, but I've already determined ahead of time that I'm not going to use gravitation; at least not in the same formulatic sense used in real life.

With regards to collisions, your equation is flawed in that checking the distances between them and accounting the radius's does not work accurately if one body or both bodies are moving fast enough to pass each other in one tick; which would cause them to have evaded collision in the code, even though they would have collided in actuality. THAT is why I'm also trying to find an equation to calculate if the coordinates of a given object cross the orbital path of another object, so that the precise and accurate collision can be detected.

And Kjartan, I haven't tested out the info in your second post yet, but I will here in a few, and thanks for the added explanations. The problem between degrees/radians would make sense to me as having thrown the results off.

just let gravity do the work for you. I've got an orbital system running on my MUD based purely off of laws of gravitation / realistic collision equations. Telnet to americanairport.us port 3000 if you'd like to have a look. You can do a 'bodyinfo 4' to look at Earth and 'bodyinfo 1' to look at the Sun, and then 'orbitinfo 4 1' to see how Earth is orbiting the Sun.

Basically I'm just calculating the gravitational influence of all the other bodies on body A, and then figuring out it's acceleration, then updating it's velocity and position based off that acceleration. Then I move onto body B and do the same. It's not taking up as much CPU as you'd think (but load will increase exponentially with more bodies.)

Wow, you are going to have so many problems with that. :-)

The first is that the number of calculations for n bodies is n(n-1). (Actually geometric, not exponential). With 2 bodies you have just two calculations. But with a thousand bodies you have nearly a million calculations per 'tick'.

The second is that you're using a discrete process to try and simulate what is actually a continuous one. Imagine two bodies heading almost directly towards each other, and very close together at the beginning of the 'tick'. The gravitational attraction at close range is very large, so the acceleration is huge, and the calculated 'next velocity' is very high.

You assume the velocity is constant for the duration of the tick, so at the beginning of the next tick, the two bodies will have passed each other and be a long way apart. But they're still modelled as travelling rapidly. So you have now have a very weak gravtitational influence slowing them down, and they keep on heading off at very high speed.

What this means is that objects consistently end up travelling faster than they should do. But the real kicker is that the effect is cumulative. Each body in the system will gradually move faster and faster and faster compared to the 'real' situation. And the faster they're moving, the more significant the effect of this error becomes. Before you know it, everything is hurtling around at relativistic velocities (or superluminal velocities if you aren't correcting for relativistic effects).

Your approach is feasible, but only if you have a really short tick length. It may be useful to have a non-uniform tick, as well - if two objects are beginning to get close to each other, use a shorter tick just for those two than you do for any object that is only affected by distant objects. Even then you will probably also will need to artificially reset the total kinetic energy of the system to its original value from time to time.


Sigma, I'm afraid your problem, as stated, is actually unsolvable. Possibly you've heard of the 'three body problem'? The end result is that you can't actually calculate precisely the trajectory of an object moving under gravity if the number of objects exerting a measurable gravitational pull is greater than 2.

You can, of course, sometimes a get a reasonable approximation by considering just two bodies. The orbit of the moon round the Earth, for example, is primarily affected by the gravity of the Earth - the influence of the Sun is relatively minor. So if you model the Earth/Moon as a single body, and calculate how its centre of mass orbits the Sun, then calculate how the Earth and Moon cause each other to orbit around their common centre of mass, you get a rough idea. Similarly, the mass of a space-ship is negligible compared to the mass of a planet, so you don't have to worry about the planet orbiting around the spaceship.

Shasharak, Sigma doesn't want to use gravity...that was all lockheed claiming it would be easier to simulate the involved physics that to completely ignore them. It sounds much more like Sigma just wants to have his planets travel orbits around certain points, which sounds easy to me.

Just imagine the orbit as a circle. Now we define a function f such that f(t,p) = the x,y coordinate pair corresponding to t hours where the orbit period is p hours (or days, weeks, months, w/e). Keeping in mind that there are 2*pi radians in a circle, our equation might look something like this:

f(t,p): x := cos(t/p*2*pi), y := sin(t/p*2*pi)

Now, it's been quite a while since I took a math course, so that might be wrong, but I see no reason it should be. On the other hand, you're not really interested in finding the x,y pair relative to 0,0 and with an orbital length of 1. So lets make the necessary modifications.

f(t,p): x := x0 + d*cos(t/p*2*pi), y := y0 + d*sin(t/p*2*pi)

x0,y0 are the CENTER of the circle, they are NOT points on the circle. Please keep that in mind. One of the nice things for such a system is that x0,y0 do not have to be constant, this equation should (If I've not been stupid) work perfectly for moons orbitting planets orbitting suns orbitting black holes. If it seems funny that a moon might be going retardedly fast relative to the black hole at times, don't forget that our earth spins at roughly 1700 km/h (and with superlow gas drainage!). That sounds pretty fast, but the solar system is travelling at 900,000 km/h relative to our galaxy (Milky Way). Our Galaxy itself is travelling over 1,000,000 km/h relative to our Local Group...so don't be surprised by fast speeds in space :) It's all relative; most people never realize they're hurtling through space at velocities that'd make bullets jealous.

(Comment added by Kelson on Fri Sep 23 6:19:20 2005)

If you're going to use it for the example I gave above, it would be necessary to update all space objects beginning with the black hole and moving outward in order of influence. Updating the moon first would just be odd...

(hierarchy: black hole, star, planet, moon)

I never addressed the collision question; here's a partial answer.

Here's how to compute whether there is a collision over a time interval. Here I assume both objects move along a linear path from their starting positions to their ending positions. Also I assume both objects are spherical. It gets more complicated if the paths aren't straight, but you can safely approximate them as straight as long as you take small enough time steps.

I will call the objects 1 and 2. They are both spheres, object 1 with radius r1 and object 2 with radius r2. For the purposes of this argument, it doesn't matter what the starting and ending time of the time interval are.

Object 1 moves from x = x1a, y = y1a at the start of the interval to x = x1b, y = y1b at the end of the interval.

Object 2 moves from x = x2a, y = y2a at the start of the interval to x = x2b, y = y2b at the end of the interval.

I will go ahead and show how I arrive at my answer, in case you are feeling adventurous and want to generate an answer for non-straight paths. You should be able to do non-straight paths with the same type of reasoning, the math is just more complicated. If you aren't interested in how I got to my answer, you can skip down to the code snippet at the end of this post.

First, I 'parametrize' the paths of the two objects over the interval. That means that I create a function of some new variable, call it s, that shows where the objects are at each point between the start and end of the interval. I will arbitrarily say that the start of the interval is at s = 0 and the end of the interval is at s = 1. I said 'a function' but actually I need two functions for each object, one for its x position and one for its y position. Recall that I said I was assuming the paths were linear, so the function should be linear, it should look like

x1( s ) = a * s + b

and since we need x( 0 ) to be x1a and x( 1 ) to be x1b, this uniquely defines a and b, and leads to

x1( s ) = x1a + s * ( x1b - x1a )

The same argument gives us

y1( s ) = y1a + s * ( y1b - y1a )

These functions x1( s ) and y1( s ) specify the path for object 1. The path for object 2 works exactly the same way, and is specified by

x2( s ) = x2a + s * ( x2b - x2a )

y2( s ) = y2a + s * ( y2b - y2a )

Next we write out the function for the distance between object 1 and object 2 as a function of s. This can be computed using the pythagorean theorem, which tells us that

dist( s ) = sqrt( ( x1(s) - x2(s) )^2 + ( y1(s) - y2(s) )^2 )

When is there a collision? There is a collision whenever dist( s ) = r1 + r2. That is, the distance between the centers of the two objects exactly equals the sum of the radii of the objects. That is

r1 + r2 = sqrt( ( x1(s) - x2(s) )^2 + ( y1(s) - y2(s) )^2 )

It's hard to work with that sqrt in there, so let's square both sides

( r1 + r2 )^2 = ( x1(s) - x2(s) )^2 + ( y1(s) - y2(s) )^2

Now I plug in the formulas for x1(s), x2(s), y1(s), and y2(s) to get the ugly equation

( r1 + r2 )^2 = ( x1a + s * ( x1b - x1a ) - x2a - s * ( x2b - x2a ) )^2
+ ( y1a + s * ( y1b - y1a ) - y2a - s * ( y2b - y2a ) )^2

To simplify the notation I am going to define some constants:
c1 = ( r1 + r2 ) * ( r1 + r2 )
c2 = x1a - x2a
c3 = x1b - x1a - x2b + x2a
c4 = y1a - y2a
c5 = y1b - y1a - y2b + y2a

Then the above equation becomes

c1 = ( c2 + s * c3 )^2 + ( c4 + s * c5 )^2

or

( c1 - c2^2 - c4^2 ) + s * ( - 2 * c2 * c3 - 2 * c4 * c5 ) + s^2 ( - c3^2 - c5^2 ) = 0

which you will note is the quadratic formula for s, a s^2 + b s + c = 0, where

a = - c3^2 - c5^2
b = - 2 c2 c3 - 2 c4 c5
c = c1 - c2^2 - c4^2

and is solved by

s = ( - b +- sqrt( b * b - 4 * a * c ) ) / ( 2 * a )

So, you would probably want to have code do all of this work. You end up with several possibilities you need to examine:

1. The two objects never collide at any s. This is true if b * b - 4 * a * c < 0.
2. The two objects collide at exactly one s. This generally isn't going to happen. It would only happen if b * b - 4 * a * c = 0.
3. The two objects collide at two s's. That happens if b * b - 4 * a * c > 0.

If we are in case 1, we're fine, there's no collision. In case 2 or 3, we need to check the s at which there are collisions and see if they lie between 0 and 1. If the collision time s is below 0, then there won't be a collision this timestep even though if you projected the linear paths of the two objects back in time, they would have collided at some point in the past. If the collision time s is greater than 1, then they won't collide during this interval although they will collide at some point in the future if they keep on their present paths.

Finally, if there are two collision s's (case 3) and they are both between 0 and 1, we would want to use the earlier one.

Here's a C code snippet that implements all of this arithmetic. I have tested it a little bit, it appears to work.

/**
 * check_collision: returns a number between 0 and 1 if the simultaneous
 * linear paths from (x1a,y1a) to (x1b,y1b) and from (x2a,y2a) to (x2b,y2b)
 * would result in a collision between objects of radius r1 and r2.
 * Returns -1 if there is no collision, returns the fraction of the step
 * (between 0 and 1) at which the collision occurs if there is a collision.
 */
double check_collision( double x1a, double y1a, double x1b, double y1b,
   double r1, double x2a, double y2a, double x2b, double y2b, double r2 ) {
   double c1, c2, c3, c4, c5, a, b, c, disc, s1, s2;

   c1 = ( r1 + r2 ) * ( r1 + r2 );
   c2 = x1a - x2a;
   c3 = x1b - x1a - x2b + x2a;
   c4 = y1a - y2a;
   c5 = y1b - y1a - y2b + y2a;
   a = - c3 * c3 - c5 * c5;
   b = - 2.0 * c2 * c3 - 2.0 * c4 * c5;
   c = c1 - c2 * c2 - c4 * c4;
   disc = b * b - 4.0 * a * c;
   if ( disc < 0 ) {
      return -1.0;
   } else if ( disc == 0 ) {
      s1 = - b / ( 2.0 * a );
      if ( ( s1 >= 0 ) && ( s1 <= 1.0 ) ) {
         return s1;
      } else {
         return -1.0;
      }
   } else {
      s1 = ( - b + sqrt( disc ) ) / ( 2.0 * a );
      s2 = ( - b - sqrt( disc ) ) / ( 2.0 * a );
      if ( ( s1 >= 0.0 ) && ( s1 <= 1.0 ) ) {
         if ( ( s2 >= 0.0 ) && ( s2 <= s1 ) ) {
            return s2;
         } else {
            return s1;
         }
      } else if ( ( s2 >= 0.0 ) && ( s2 <= 1.0 ) ) {
         return s2;
      }
   }
}

(Comment added by Kjartan on Fri Sep 23 12:07:46 2005)

Ah, I see I forgot a case in the if's at the end of the code snippet. After that last "else if"'s closing }, there should be

      } else {
         return -1.0;
      }