Little orbit sim

[george7378]

Hi everyone,

I made a little simulator recently, just to test out some stuff in a book I was lent. It was done in C++, and it shows how an orbit progresses over time based on a few inputs. It's very basic (it came out at less than 100 lines of code in the end) but I'm quite proud of it. There are also many more features I'm planning to add (such as launching into an orbit first, and the ability to change orbits mid-sim) but I thought I'd share it for others to have a look at:

http://gkastro.0sites.org/OrbitSim.exe

Maybe it will develop into something great

 

[Wishbone]

If you have access to any graphing component, don't hesitate to use it. Besides obvious users' pleasure, there's the added benefit of being able to spot for bugs and inconsistencies visually.

 

[george7378]

That's a good idea - I think it's probably beyond my experience at the moment though. The book also tells you how to get x and y coordinates for the satellite at any time, so that would allow it to plot the shape of the orbit too.

 

[Wishbone]

If you have arrays of the coordinates handy, the quick and dirty way to display a static graph is using DISLIN: http://www.mps.mpg.de/dislin/

 

[Spacethingy]

Cool. With a Windows Forms that could be animated...
:thumbup:

 

[Wishbone]

I did not touch upon .NET because for some it is a touchy question. But yes, there're bunches of components like NPlot or packs from Microsoft Research etc.

 

[luki1997a]

Can you explain how does it work? Any algorythm or something?

 

[george7378]

Here's the source code with comments explaining it:

#include <iostream>
#include <cmath>
#include <windows.h>
#include <conio.h>

#define Rearth             6378.14
#define GMearth            398600.4415
#define Pi                 3.141597

using namespace std;

int main (){

	int runprogram = 1;
	double perigee, apogee, rapogee, rperigee;
	double time = 0;
	double timerate = 1;
	char key;

	cout << "Enter desired apogee altitude in KM: " << endl;
	cin >> apogee;
	cout << "Enter desired perigee altitude in KM: " << endl;
	cin >> perigee;

	rapogee = apogee + Rearth;   //Radius of Earth + apogee
	rperigee = perigee + Rearth;   //Radius of Earth + perigee
    
	double semimajor = (rapogee + rperigee) / 2;         //Semi-major axis
	double eccentricity = (rapogee - rperigee) / (rapogee + rperigee);     //Eccentricity of orbit
	double meanmotion = sqrt((GMearth / (semimajor*semimajor*semimajor)));    //Mean motion (vector used to calculate how much of an orbit has been completed during a time interval)

	while (runprogram == 1){          //Loops the program

		if (_kbhit())
		{key = _getch();
		if (key == 114)
			timerate = timerate - 1;
		else if (key == 116)
			timerate = timerate + 1;             //Time acceleration controls

		if (timerate < 1)
			timerate = 1;}                  //Time acceleration can't fall below 1

		time = time + timerate;              //Progression of time as program loops

		double meananomaly = meanmotion * time;              //Mean anomaly (fraction of orbit completed)

		double eccentricanomaly0 = meananomaly;             //Iteration to work out Eccentric Anomaly angle

        double Margin = 1.0E-8; 

              while (true)
			  {
              double eccentricanomaly = eccentricanomaly0 - (eccentricanomaly0-eccentricity*sin(eccentricanomaly0)-meananomaly)/(1.0-eccentricity*cos(eccentricanomaly0));

              if (abs(eccentricanomaly - eccentricanomaly0) < Margin)
              {
              eccentricanomaly0 = eccentricanomaly;

              break;
              }

              eccentricanomaly0 = eccentricanomaly;            //End of iteration
              }

	    double trueanomaly = acos((cos(eccentricanomaly0) - eccentricity) / (1 - (eccentricity * cos(eccentricanomaly0))));             //True anomaly (angle between perigee and current point in orbit)

		double radius = semimajor * (1 - (eccentricity * cos(eccentricanomaly0)));
		double altitude = radius - 6378.14;        //Altitude calculation

		double flightpath = atan((eccentricity * sin(trueanomaly)) / (1 + (eccentricity * cos(trueanomaly))));
		double velocity = (sqrt(GMearth * semimajor * (1 - (eccentricity*eccentricity)))) / (radius * cos(flightpath));                //Velocity calculation
		int period = 2 * Pi * (sqrt((semimajor*semimajor*semimajor)/GMearth));          //Orbit period
		double completeorbits = time / period;        //Complete orbits

		cout << "To slow time, press 'r'\n" << "To speed up time, press 't'\n\n" << endl;
		cout << "Current vectors:\n" << endl;
		cout << "Altitude: " << altitude << " KM" << endl;
		cout << "Velocity: " << velocity << " KM/Sec\n" << endl;
		cout << "Orbital period: " << period << " Sec" << endl;
		cout << "Complete orbits: " << completeorbits << endl;           //Outputs

		system ("cls");}

return 0;
}

It basically just loops that set of equations and outputs the result at every step.

---------- Post added at 15:03 ---------- Previous post was at 14:23 ----------

I wonder if Orbiter has the same equations at its core?

 

[RisingFury]

Your Pi is incorrect...

3.141593, not 3.14159

Rounded, of course...

---------- Post added at 16:25 ---------- Previous post was at 16:23 ----------

I wonder if Orbiter has the same equations at its core?

No, it doesn't.

Instead of calculating the velocity and assuming the object will fly in an ellipse, it sums all the forces acting on the body to get acceleration. Then, knowing acceleration, velocity and current position, it propagates the state into next frame.

You'll find a more detailed explanation of how that's done in
\Orbiter\Doc\Technotes\dynamics.pdf

 

[luki1997a]

:hello:
I've created a gui program based on your source:

Many thanks for the source:thumbup:, so I should ask you: Can I release it?

 

[Jarvitä]

Here's the source code with comments explaining it:

Why did you write a C program in C++?

 

[Wishbone]

Why did you write a C program in C++?

One li'l step at a time . Learning doesn't come in an instant...

 

[luki1997a]

So I can release it right?

 

[george7378]

:hello:
I've created a gui program based on your source:

Many thanks for the source:thumbup:, so I should ask you: Can I release it?

Very sorry for the delay! Yes, it's OK for you to release it. If you used any of the source directly, a little mention would be nice, but it doesn't really matter - it's your program. Thanks for making my code more presentable! I would love to learn to program graphics like that - did you use C++?

Why did you write a C program in C++?

Not really sure what you mean by that! Have I used a mixture of the two languages? I learned C++ from an Internet tutorial anyway, so I apologise for any silly mistakes! I always use that syntax and structure - it works OK (most of the time).

 

[n72.75]

I don't see where he used depreciated syntax.

It looks fine to me.

 

[george7378]

I was also thinking - it would be quite cool to add an atmosphere simulation to it - maybe I could try and use the fluid dynamics drag equation, and define different densities at different altitudes, and use F=ma to change the velocity based on the drag.

 

[luki1997a]

Thanks george, I'll release it soon.:thumbup:

About coding langauge:
It's MS Visual Basic 2010 Express
It's very easy, easier than any other langauge(I think), because you just drag'n'drop buttons, textboxes etc. to your Form. It uses .NET Framework 4.

Atmosphere simulation:
There's a little blue line which indicates atmosphere up to 200 kms.

You can ask me about everything you want

 

[Turbinator]

Is this Orbiter 0.0001 ?

Neutrinos... neutrinos everywhere.

 

[luki1997a]

:lol:
'Simple Orbit Simulator':thumbup:

 

[RisingFury]

I was also thinking - it would be quite cool to add an atmosphere simulation to it - maybe I could try and use the fluid dynamics drag equation, and define different densities at different altitudes, and use F=ma to change the velocity based on the drag.

For a compressible fluid such as gas, the simplest estimate of density with height is the Rho = Rho0 * e^-(z / z0) formula. Rho0 is density at ground level and z0 controls how quickly the density drops with altitude. This holds well for low altitudes, but not very well for high altitudes in areas with low density. Also, it doesn't take into account temperature changes at different altitudes.

As for the force of drag...
F = 1/2 * c * Rho * S * v^2 equation works well in the subsonic regime up to about 250 m/s and then starts to break down. Supersonic drag is a lot more complex, but if you want to implement it, there are sources out there, if you look around.