Modeling stellar orbits

Started by
94 comments, last by taby 13 years, 5 months ago
Alright, here's a new version. The feature set is exactly the same, actually lesser since I don't include the controls at the moment. I'm more interested to know if it runs for you JoeCooper.

A little more cleanup and then it's back to implementing more physics-y stufff. I'll save my next post for that unless something comes up!

Oh, I reversed the controls for zooming/changing spawn mass. Now you hold control and scroll the mouse wheel to up the mass.

Edit: I suppose I could create a new thread for this question, but seeing as it is related to modeling stellar orbits I guess I can just put it here:

@h4tt3n

How did you calculate the Periapsis and Apoapsis of your orbits? I wanna do some extra stuff depending on the eccentricity of the orbit in real time, but I can't calculate it without those values correct?

Edit2: Okay, so I found out how to calculate the eccentricity of my orbit based on the orbital state vectors. However I'm getting some strange eccentricities. If the stars aren't actually in orbit and just go shooting off, the eccentricity is greater than 1 like it should be (although always just barely, as in 1.0004) however when an orbit is sustained, the eccentricity is always 0.99xxx. I was under the impression that a perfectly circular orbit should have eccentricity of 0, yet I cannot get even close to this result no matter how circular my orbit is.

I'm gonna sleep on it and then maybe try to work out the bugs myself, if not I'll probably post my code here to see if any of you guys can help!

[Edited by - ChugginWindex on December 2, 2010 10:29:03 PM]
Advertisement
It works!

It was missing some runtime DLLs but I downloaded them. MSVCsomethingsomething and MSVCsomethingelse.
Oh whoops! Yeah those are the two that I made sure were with the last version every time I uploaded it. Just forgot this time...

I fixed the link in my previous post, it points to one that has those now.
Quote:Original post by h4tt3n
@Chuggin

I've been working with this stuff for years. Don't ever let it discourage you that others have come further than you. I don't mind posting the source, but it won't help you unless you know a little about orbital mechanics. Google that and orbital elements, and soon you'll have figured it out yourself. It's much more rewarding that way.


Mike, thank you for posting that encouragement!
Graham Rhodes Moderator, Math & Physics forum @ gamedev.net
Okay, so I've got two major problems in calculating the ellipses for my orbits:

The first is that my eccentricity, semi-major axis and all other orbital elements are being calculated absolutely correctly if the the orbiting star has a mass of 1.0. Anything above this gives me huge negatives for eccentricity (around -20). I think it has to do with how I am calculating my orbital angular momentum. I'm posting the code that does all the work in the hopes that someone might spot something I'm missing. (NOTE: I can't seem to find readable information on how to actually do all of this math, so I actually found a website that generates the orbital elements from the orbital state vectors and viewed the page source, haha)

OrbitalElements Body::getOrbitalElements(Body& centralBody) {	OrbitalElements elems;	elems.MajorMass = centralBody.getMass();	elems.MinorMass = m_mass;	elems.RelativePosition = m_position.distanceFrom(centralBody.getPosition());	elems.Velocity = m_velocity;	elems.V2 = m_velocity.distanceSquared(Vector2f(0.0f, 0.0f)); //the velocity squared, for ease	elems.LinearMomentum = Vector2f(m_velocity.X*m_mass, m_velocity.Y*m_mass);	elems.AngularMomentum = elems.RelativePosition.X*elems.LinearMomentum.Y - elems.RelativePosition.Y*elems.LinearMomentum.X;        /* This commented out version of calculating AngularMomentum was one that wikipedia suggested, however when I did it this way my eccentricities were always incorrect */	//elems.AngularMomentum = (elems.RelativePosition.X*elems.Velocity.Y + elems.RelativePosition.Y*elems.Velocity.X);	elems.Mu = m_gravityConstant * (elems.MajorMass + elems.MinorMass); //standard gravitational parameter	float p = elems.AngularMomentum * elems.AngularMomentum / elems.Mu; //umm, I dont actually kow what "p" is, but it seems to be used in the code that I looked at for inspiration, and it all works when MinorMass==1.0!	float r2 = elems.RelativePosition.X * elems.RelativePosition.X + elems.RelativePosition.Y * elems.RelativePosition.Y; //the cartesian position vector	elems.R = sqrt(r2); 	elems.SemiMajorAxis = 1.0f/(2.0f / elems.R - elems.V2 / elems.Mu);	elems.Eccentricity2 = 1.0f - p/elems.SemiMajorAxis;	return elems;}


Problem #2 is with the ellipses themselves. Like I said all my math comes out correct for satellites of mass 1.0. I've been trying to draw the orbital ellipses in this case, and I'm having a hard time orienting my ellipse so that the central body is always at one of the foci. I honestly have no idea on how to do this, I thought I could work it out using the normal of the relative position between the satellite and the central star, but that doesn't actually solve the problem.

I'll post the code I'm using at the moment to draw the ellipse as well:

OrbitalElements elems = (*iter)->getOrbitalElements(centralBody);			float ec2 = elems.Eccentricity2;			if(ec2 < 1.0f) {				//we've got something to draw!				float semiminorAxis = elems.SemiMajorAxis*sqrt(1.0f - ec2);				Vector2f normalized = elems.RelativePosition;				normalized.X /= elems.R;				normalized.Y /= elems.R;				float xs = pos.X + normalized.X*elems.SemiMajorAxis;				float ys = pos.Y + normalized.Y*elems.SemiMajorAxis;				glBegin(GL_LINE_LOOP);					for(float i = 0.0f; i < 3.141592654f*2.0f; i+= (3.141592654f/40.0f)) {						glVertex2f(xs + cos(i)*elems.SemiMajorAxis,ys + sin(i)*semiminorAxis);					}				glEnd();			}
Hi there, a few comments on the code:

I think your elems.Velocity should be relative to the central body, otherwise it will only work if the central body stands completely still relative to the reference frame.

elems.V2 = m_velocity.distanceSquared(Vector2f(0.0f, 0.0f)); This strikes me as an odd way of getting velocity squared. Isn't there a magnitudesquared function?

Besides, neither your orbital element function nor your ellipse drawing function seem to handle rotated orbits, so it won't be able to draw an orbit unless the semimajor axis is horizontal.

As for problem #2, you need to calculate the distance between the focii, which is simply eccentricity*semimajor axis (iirc). Then simply put the ellipse center half that distance away from the central body. Again, rotated orbits not accounted for.

Oh, would you mind posting a link to the website you mentioned? Just want to check it out.

Cheers,
Mike
Sure, the link to the site I was using the page source for was:

http://orbitsimulator.com/formulas/OrbitalElements.html

My "strange" way of calculating Velocity squared is because I wrote my own Vector class. All that function does is calculate the delta between two vectors and then do all of calculating magnitude EXCEPT for the actual square root at the end of pythagorean's theorem. By passing the "from vector" as (0, 0) it's the same as saying elems.Velocity.X*elems.Velocity.X + elems.Velocity.Y*elems.Velocity.Y...I just didn't feel like writing all that in the code so I made use of an existing function to do it...

You're right, I'm not using velocity relative to the central body, I had no idea I had to although it seems obvious now that you mention it.

I hate to ask so much cuz I really do like learning on my own, but I've been trying to get this working for a while now. Would you mind giving me some insight into how I would go about "handling rotated orbits". Like I said my ellipse is perfectly fine for a mass of 1.0, but it doesn't rotate.

Thanks!
My first tries at this didn't work very well either, so take it cool :-)

As for actually drawing a rotated ellipse, there's probably loads of code samples already out there, but if you keep beeing stuck I don't mind giving you my home-made and very fast routine later on.

As for calculating a rotated trajectory, I'll give you just one hint: eccentricity-vector

Good huntin'

Cheers,
Mike
Does that mean that I can't just obtain the relative velocity by subtracting the central body's velocity from the satellite's velocity at every instant in time?

Edit: Okay, I repositioned the ellipse so that it now draws like this:

float ecc = sqrt(ec2);				float semiminorAxis = elems.SemiMajorAxis*sqrt(1.0f - ec2);				Vector2f focusPos = m_centralBody->getPosition();				Vector2f normalized = elems.RelativePosition;				normalized.X /= elems.R;				normalized.Y /= elems.R;				float xs = focusPos.X - normalized.X*elems.SemiMajorAxis*ecc/2;				float ys = focusPos.Y - normalized.Y*elems.SemiMajorAxis*ecc/2;				glBegin(GL_LINE_LOOP);					for(float i = 0.0f; i < 3.141592654f*2.0f; i+= (3.141592654f/40.0f)) {						glVertex2f(xs + cos(i)*elems.SemiMajorAxis,ys + sin(i)*semiminorAxis);					}				glEnd();


However like you said, my rotations are still far off because I'm not doing anything with them at all. I'm going to start looking around for examples but so far I haven't actually found anything. There's tons of info on drawing ellipses with horizontal and vertical semimajor axis, but nothing that really does the trick for rotated ones. I'm going to keep looking...

Also, I'm still struggling to see the connection between the eccentricity vector and the rotational trajectory. Wikipedia doesn't seem to have any info on that either.

[Edited by - ChugginWindex on December 3, 2010 2:13:17 PM]
Yes you can. In code it should look something like:

elems.Velocity = m_velocity - centralBody.velocity;

This topic is closed to new replies.

Advertisement