# Orbital dynamics math review for a computer simulation

ElectronicStudent123
Homework Statement:
I have been tasked with creating a computer simulation that (algebraically) models the behavior of an asteroid as it approaches the Earth. I have a drawing area of 300 Mm by 200 Mm, which I scale down to be 1200 by 800 on screen and 12 by 8 on paper. The asteroid begins its journey in Q3 at point (-150Mm x, -100Mm y) and moves towards the Earth which its situated at the origin.

There are two variables associated with the asteroid: initial angle and initial speed. The angle can be between 10 degrees to 80 degrees and the speed can be between 1 km/s to 20 km/s. For programming purposes, I use scanf() to input the initial values, but the final product will use potentiometers to set the values.

The Earth is designated as point P(x,y) and the asteroid is designated as point C(x,y). The mass of the Earth is #define Mp 5.972e24 and the mass of the asteroid is #define Mc 26.99e9.

The purpose of this exercise is to better understand computing concepts like multitasking, graphics, multithreading, semaphores, IPCs, etc.
Relevant Equations:
Vcx = Vc * cosA
Vcy = Vc * sinA
F = (G * Mp * Mc)/(R^2)
R^2 = (Px - Cx)^2 + (Py - Cy)^2
B = arctan((Py - Cy)/(Px - Cx)
ax = (F/Mc) * cosB
ay = (F/Mc) * sinB
Vcx = V0 + (ax * t)
Vcy = V0 + (ay * t)
Cx = Cx-old + (Vcx * t)
Cy = Cy-old + (Vcy * t)
The program begins by prompting the user to input the initial angle and speed. Once these values have been input, the program calculates the initial velocity components for the asteroid in two dimensions.

C:
printf("Input an initial angle from 10 to 80 degrees\n");
scanf("%f", &theta);
rewind(stdin);
printf("Input an speed from 1000 m/s to 20,000 m/s\n");
scanf("%f", &V0);
rewind(stdin);
printf("You have entered: %.1f degrees and %.1f km/s\n", theta, V0);
printf("Computing initial x and y components...\n");
Vcx = V0 * (cos(theta));
Vcy = V0 * (sin(theta));
printf("Vcx = %.1f m/s and Vcy = %.1f m/s\n\n", Vcx, Vcy);

Once the initial velocity components have been calculated, a pipe for passing information between threads is created, along with two threads. The calculations are handled in the parent thread and the graphics are handled in the child thread. If any potential helpers would like me to post the entire code I can, but for now I am just going to post the calculations thread and its relevant code.

C:
#define G 6.67408e-11 //universal graviational constant
#define Mp 5.972e24 //mass of the earth
#define Mc 26.99e9 //mass of apophis
#define ts 1200
#define pi 3.14159
#define todeg 180/pi
#define scalar 250000

const double Px = 0;
const double Py = 0;
double Cx = -150e6;
double Cy = -100e6;
float ax = 0;
float ay = 0;
float a = 0;
float v = 0;
float Fx, Fy;
float Q1xy[2];
float Data[4]; //for use with DisplayData

while (endflag != NULL)
{
endflag = fopen("666", "r");
double Rsq, F;
float B = 0;
float cosA = 0;
float sinA = 0;
float A1 = 0;
float A2 = 0;
float A3 = 0;
float A4 = 0;

if(Cx > 0 && Cy > 0) //Q1
{
printf("***YOU ARE IN Q1***\n");
printf("Cx-now = %g Mm and Cy-now = %g Mm\n", Cx, Cy);
Rsq = ((Px - Cx) * (Px - Cx)) + ((Py - Cy) * (Py - Cy));
printf("Rsq = %g m^2\n", Rsq);
F = (G * Mp * Mc) / Rsq;
printf("Force = %g N\n", F);
B = BCalc(Cx, Cy); //BCalc is a function which calculates the angle and forces it positive if it is negative
A1 = B;
cosA = cos(A1) * todeg;
printf("cosA = %.1f deg\n", cosA);
sinA = sin(A1) * todeg;
printf("sinA = %.1f deg\n", sinA);
ax = (F / Mc) * cos(A1);
ay = (F / Mc) * sin(A1);
if(Cx < Px && Cy < Py)
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
else
{
Vcx = Vcx - (ax * ts);
Vcy = Vcy - (ay * ts);
}

}
else if(Cx < 0 && Cy > 0) //Q2
{
printf("***YOU ARE IN Q2***\n");
printf("Cx-now = %g Mm and Cy-now = %g Mm\n", Cx, Cy);
Rsq = ((Px - Cx) * (Px - Cx)) + ((Py - Cy) * (Py - Cy));
printf("Rsq = %g m^2\n", Rsq);
F = (G * Mp * Mc) / Rsq;
printf("Force = %g N\n", F);
//B = atan((Py - Cy) / (Px - Cx));
B = BCalc(Cx, Cy);
A2 = B + (pi / 2);
cosA = cos(A2) * todeg;
printf("cosA = %.1f deg\n", cosA);
sinA = sin(A2) * todeg;
printf("sinA = %.1f deg\n", sinA);
ax = (F / Mc) * cos(A2);
ay = (F / Mc) * sin(A2);
if(Cx < Px && Cy < Py)
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
else
{
Vcx = Vcx - (ax * ts);
Vcy = Vcy - (ay * ts);
}

}
else if(Cx < 0 && Cy < 0) //Q3 (asteroid starts here)
{
printf("***YOU ARE IN Q3***\n");
printf("Cx-now = %g Mm and Cy-now = %g Mm\n", Cx, Cy);
Rsq = ((Px - Cx) * (Px - Cx)) + ((Py - Cy) * (Py - Cy));
printf("Rsq = %g m^2\n", Rsq);
F = (G * Mp * Mc) / Rsq;
printf("Force = %g N\n", F);
//B = atan((Py - Cy) / (Px - Cx));
B = BCalc(Cx, Cy);
A3 = B + pi;
cosA = cos(A3) * todeg;
printf("cosA = %.1f deg\n", cosA);
sinA = sin(A3) * todeg;
printf("sinA = %.1f deg\n", sinA);
ax = (F / Mc) * cos(A3);
ay = (F / Mc) * sin(A3);
if(Cx < Px && Cy < Py)
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
else
{
Vcx = Vcx - (ax * ts);
Vcy = Vcy - (ay * ts);
}
}
else if(Cx > 0 && Cy < 0) //Q4
{
printf("***YOU ARE IN Q4***\n");
printf("Cx-now = %g Mm and Cy-now = %g Mm\n", Cx, Cy);
Rsq = ((Px - Cx) * (Px - Cx)) + ((Py - Cy) * (Py - Cy));
printf("Rsq = %g m^2\n", Rsq);
F = (G * Mp * Mc) / Rsq;
printf("Force = %g N\n", F);
//B = atan((Py - Cy) / (Px - Cx));
B = BCalc(Cx, Cy);
A4 = B + ((3 * pi) / 2);
cosA = cos(A4) * todeg;
printf("cosA = %.1f deg\n", cosA);
sinA = sin(A4) * todeg;
printf("sinA = %.1f deg\n", sinA);
ax = (F / Mc) * cos(A4);
ay = (F / Mc) * sin(A4);
if(Cx < Px && Cy < Py)
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
else
{
Vcx = Vcx - (ax * ts);
Vcy = Vcy - (ay * ts);
}
}
a = sqrt((ax * ax) + (ay * ay));
v = sqrt((Vcx * Vcx) + (Vcy * Vcy));
Cx += (Vcx * ts);
Cy += (Vcy * ts);
printf("acceleration = %.3f m/s^2\n", a);
printf("velocity = %.3f m/s\n", v);
printf("ax = %.3f m/s^2 and ay = %.3f m/s^2\n", ax, ay);
printf("Vcx = %.3f m/s and Vcy = %.3f m/s\n", Vcx, Vcy);
printf("Cx-new = %g Mm and Cy-new = %g Mm\n", Cx, Cy);
Fx = DownSizer(Cx);
Fy = DownSizer(Cy);
printf("Downsizing Cx and Cy yields: %.1f and %.1f\n", Fx, Fy);
Data[0] = Fx;
Data[1] = Fy;
Data[2] = Vcx;
Data[3] = Vcy;
char message[20];
memset(message, '\0', sizeof(message));
sprintf(message, "%.1f %.1f", Fx, Fy);
printf("Message = %s\n\n", message);
write(id[1], message, strlen(message)); //write message to pipe
sleep(1);
startgfx = fopen("1", "w"); //create flag to tell gfx to start the display
sleep(1);

}

float BCalc(double Cx, double Cy)
{
float tempB = 0;
tempB = atan((Py - Cy) / (Px - Cx));
if(tempB < 0)
{
tempB *= -1;
return tempB;
}
else return tempB;
}

The calculations thread enters a while() loop that is terminated by the removal of an empty file in the directory in which the code is located. This empty file is removed if a) the asteroid leaves the drawing area b) the asteroid collides with the planet c) the asteroid enters a stable orbit. I have not been able to achieve a stable orbit yet, which is why I am here.

The calculations loop begins by determining what quadrant the asteroid is in by using a series of conditional statements.

C:
if(Cx > 0 && Cy > 0) //Q1
else if(Cx < 0 && Cy > 0) //Q2
else if(Cx < 0 && Cy < 0) //Q3 (asteroid starts here)
else if(Cx > 0 && Cy < 0) //Q4

Regardless of which quadrant the asteroid is in, the square of the distance between the asteroid and the Earth are calculated, as well as the force between them. I have a function called BCalc() which takes the current position of Cx and Cy as its arguments. This function calculates the angle between the asteroid and the Earth, and if the angle is negative, forces it to be positive (I am not sure if this is necessary).

C:
float BCalc(double Cx, double Cy)
{
float tempB = 0;
tempB = atan((Py - Cy) / (Px - Cx));
if(tempB < 0)
{
tempB *= -1;
return tempB;
}

Depending on which quadrant the asteroid is in, I declare a new angle, either A1, A2, A3 or A4 for use in the remaining calculations.

C:
A1 = B;
A2 = B + (pi / 2);
A3 = B + pi;
A4 = B + ((3 * pi) / 2);

My program then proceeds to calculate the acceleration for the x-axis and the y-axis based on the cos and sin of either A1, A2, A3, A4. With these accelerations, I then update the x and y velocity components based on two conditional statements.

C:
if(Cx < Px && Cy < Py)
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
else
{
Vcx = Vcx - (ax * ts);
Vcy = Vcy - (ay * ts);
}

These conditionals are the same regardless of whatever quadrant the asteroid is in. Perhaps the problem with my code is here?

My drawing area is 300 Mm by 200 Mm, which I scale down to 1200 by 800. From a calculations perspective, I treat this as if there are four quadrants. Unfortunately, the graphics library that we have to use, initializes a drawing area that only behaves as if there is one quadrant (Q1) so I have to "shift" the Earth and the position of the asteroid relative to the Earth to Q1 for the purpose of displaying the graphics. This may be a source of error, but I believe that I have this part correct. The Earth is initially at the origin, and it is shifted to the middle of the drawing area (600, 400). Since the position of the asteroid will always be relative to the Earth, I take the x and y coordinates of the asteroid and add 600 to the x coordinate and 400 to the y coordinate to get the equivalent Q1 coordinate for use in the graphics thread.

Attached is an output of the program where we can see the path of the asteroid being altered by the Earths gravity (25 degrees initial angle with 8000 m/s initial speed). The second attachment is of my four quadrants for the calculations thread and the one quadrant for the graphics thread.
If you have any other questions, please let me know. I am making this thread in hopes that someone with more knowledge about this stuff can review my math and help find any errors that exist. Thanks for reading.

#### Attachments

• 1.jpg
24.8 KB · Views: 90
• 2.jpg
45.7 KB · Views: 74
Last edited by a moderator:

ElectronicStudent123
I have made a couple changes to the code, which seem to have had a positive impact, but a stable orbit still eludes me.

I have replaced all instances of:
C:
if(Cx < Px && Cy < Py)
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
else
{
Vcx = Vcx - (ax * ts);
Vcy = Vcy - (ay * ts);
}
with a single instance of:
C:
if (Cx > Px || Cy > Py)
{
Vcx -= (ax * ts);
Vcy -= (ay * ts);
}
else
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}
With the exception of Q1 where A1 = B, I have replaced the following:
C:
A2 = B + (pi / 2);
A3 = B + pi;
A4 = B + ((3 * pi) / 2);

with:
C:
A2 = pi - B;
A3 = B - pi;
A4 = (2 * pi) - B;

These changes I believe make sense. Attached is an image of the path the asteroid takes with 30 degrees as the initial angle and 8000 m/s as the initial speed. If I leave the angle as 30 degrees and increase the initial speed, similar effects happen, but instead of looping around, the path is altered and the asteroid keeps going past the planet. The "tightness" of the curl as the asteroid wraps around planet is definitely being affected by the speed, which I am fairly certain makes physical sense.

#### Attachments

• 3.jpg
16.5 KB · Views: 81
Homework Helper
Gold Member
If I were to create something like this, I would start by trying to simulate a situation that I know the answer to just to make sure that I got all the parts working as they should and all bugs are removed. Why don't you start your asteroid at some known distance ##r## on the y-axis with initial speed ##v=\sqrt{\dfrac{GM_E}{r}}## directed along the x-axis.
If all works according to design Then
the asteroid should follow a circular path
Else
start troubleshooting knowing what the path should look like
End if

ElectronicStudent123
If I were to create something like this, I would start by trying to simulate a situation that I know the answer to just to make sure that I got all the parts working as they should and all bugs are removed. Why don't you start your asteroid at some known distance ##r## on the y-axis with initial speed ##v=\sqrt{\dfrac{GM_E}{r}}## directed along the x-axis.
If all works according to design Then
the asteroid should follow a circular path
Else
start troubleshooting knowing what the path should look like
End if
Hello and thanks for the reply. I think that early on in this project the idea of starting with a known angle and speed that would yield a stable orbit crossed my mind, but I quickly discarded the idea when I realized that I do not know how to calculate a stable orbit.

By setting Cy = 25e6 and Cx = 0, the asteroid starts just north of the Earth. Using your equation with r = 25e6, I find an initial speed for the x-component to be 3994 m/s.

Assuming that the angle is 0 here... I know that Vcx = Vc * cos(0) which means that Vcx = Vc. The initial speed input of my program uses the total magnitude of both the Vx and Vy components. In this case for an angle of 0, Vcy = Vc * sin(0) which equals 0. So the total magnitude is just Vcx in this case? I think that this makes sense but it has been a while since I have done anything like this.

In any event, if Cx = 0 and Cy = 25e6, and with an initial angle of 0 and an initial speed of 3994, the asteroid makes a perfect half circle from (0, 25e6) to (0, -25e6). Once the asteroid crosses back into Q3, it all falls apart and I lose my potentially stable orbit. I have included a picture of the output for your suggestion to this post (0-3994.png).

I have included another picture, called 30-8000.png, which shows the asteroid beginning from the bottom left of the screen as I initially had it. The asteroid passes from Q3 into Q4 into Q1 and into Q2, in what looks like a pretty healthy circle. Once it crosses back into Q3, I once again lose my potentially stable orbit.

What your suggestion and my initial attempt both have in common is that they both fall apart in Q3. I have not solved this problem yet, but I believe that the issue has to do with this bit of code:
C:
if (Cx > Px || Cy > Py)
{
Vcx -= (ax * ts);
Vcy -= (ay * ts);
}
else
{
Vcx += (ax * ts);
Vcy += (ay * ts);
}

The if statement is true for Q4, Q1 and Q2, but once the asteroid enters Q3, the else statement takes over which causes an increase in both the x and y components of the velocity.

The code posted above is valid for the beginning of a run, where the asteroid starts from the bottom left and picks up speed as it approaches the planet and wraps around from either the left or the right. The issue here, I believe, is that once the asteroid wraps back around to Q3, the code forces the asteroid back from the direction in which it started from.

My solution to this problem is still a work in progress, but I would like to thank you for your suggestion. Knowing that I lose my potentially stable orbit as the asteroid passes from Q2 into Q3, and seeing the asteroid lose its potentially stable orbit as it passes from Q4 into Q3 as part of your suggestion really makes me suspect that the issue is with the transition back into Q3.

#### Attachments

• 0-3994.png
115.1 KB · Views: 72
• 30-8000.png
82.2 KB · Views: 75
Last edited:
Homework Helper
Gold Member
I have not solved this problem yet, but I believe that the issue has to do with this bit of code:
It probably does. How do you update ax and ay? For a circular orbit, assuming counterclockwise motion,
In Q1 ax <0, ay<0
In Q2 ax >0, ay<0
In Q3 ax >0, ay>0
In Q4 ax<0, ay>0

Your if statement doesn't seem to be ensuring this.

Homework Helper
Gold Member
2022 Award
With the asteroid at coordinates (x,y) relative to Earth as origin,
##k=\frac{Gm_E}{(x^2+y^2)^\frac 32}##
##\dot x-=kx\Delta t##
Etc.

A key problem with simulation is accumulated error. You can minimise this in various ways.
One is to keep checking the angular momentum and energy and making small adjustments to try to keep those correct. Exactly how you should make the adjustments I'm not sure.
Another is, when updating the position variables, take the current velocity as somewhere between the previously calculated velocity and the newly calculated velocity.

ElectronicStudent123
It probably does. How do you update ax and ay? For a circular orbit, assuming counterclockwise motion,
In Q1 ax <0, ay<0
In Q2 ax >0, ay<0
In Q3 ax >0, ay>0
In Q4 ax<0, ay>0

Your if statement doesn't seem to be ensuring this.

I update the accelerations inside of each conditional statement depending on which quadrant I am in.

C:
A1 = B;
A2 = pi - B;
A3 = B - pi;
A4 = (2 * pi) - B;

Depending on which quadrant I am in, I add or subtract 180 degrees or 360 degrees from the calculated angle B to adjust for my quadrant. Then I update the acceleration with: ax = (F / Mc) * cos(A) and ay = (F / Mc) * sin(A) where A can be A1, A2, A3 or A4. I am starting to question whether anything that I have done makes sense.

With the asteroid at coordinates (x,y) relative to Earth as origin,
##k=\frac{Gm_E}{(x^2+y^2)^\frac 32}##
##\dot x-=kx\Delta t##
Etc.

A key problem with simulation is accumulated error. You can minimise this in various ways.
One is to keep checking the angular momentum and energy and making small adjustments to try to keep those correct. Exactly how you should make the adjustments I'm not sure.
Another is, when updating the position variables, take the current velocity as somewhere between the previously calculated velocity and the newly calculated velocity.

Can you tell me more about the equations that you have posted? I have never seen them before. I will try to implement these suggestions.

When calculating the square of the distance between the asteroid and the Earth with R^2 = (Px - Cx)^2 + (Py - Cy)^2 does the size of Px and Cx matter? I don't think that it does because you end up squaring them anyway and forcing them positive. The same seems to be true for calculating the initial angle with B = arctan((Py - Cy)/(Px - Cx)), with the exception of the squares.

Homework Helper
Gold Member
Can you explain why you need Px and Cx? It seems to me that only one is needed, ##x##, the x-component of the position vector from the force center (here the center of the Earth) and the asteroid. Likewise for ##y##. The formula suggested by @haruspex in #6 is based on this idea and that's the meaning of ##x## and ##y## therein. Note that these are algebraic quantities, meaning they can be positive or negative. In Q1 both are positive, in Q2 ##x<0## and ##y>0## and so on. Using ##x## and ##y## in this fashion takes care of the quadrant problem automatically without any if statements and also allows finding the angle by using the appropriate arctangent.

ElectronicStudent123
Can you explain why you need Px and Cx? It seems to me that only one is needed, ##x##, the x-component of the position vector from the force center (here the center of the Earth) and the asteroid. Likewise for ##y##. The formula suggested by @haruspex in #6 is based on this idea and that's the meaning of ##x## and ##y## therein. Note that these are algebraic quantities, meaning they can be positive or negative. In Q1 both are positive, in Q2 ##x<0## and ##y>0## and so on. Using ##x## and ##y## in this fashion takes care of the quadrant problem automatically without any if statements and also allows finding the angle by using the appropriate arctangent.

I was under the impression that I need Px and Cx, as well as Py and Cy, to calculate the distance between the asteroid and the Earth. The distance between the asteroid and the Earth are needed to calculate the force exerted on the asteroid by the Earth? And with this force, which changes with each iteration of the loop, I can update the acceleration components? Is this wrong?

What is k in the equation that haruspex mentioned?

Homework Helper
Gold Member
2022 Award
What is k in the equation that haruspex mentioned?
It is what it says in my post. You can calculate it once at each iteration then use it in updating both x and y velocities.

Last edited:
Homework Helper
Gold Member
You still have not explained what Cx, Px and Cy, Py stand for. The gravitational attraction force in Cartesian ##x## and ##y## coordinates is given by $$\vec F=-GM_Em\frac{x~\hat x+y~\hat y}{(x^2+y^2)^{3/2}}.$$Here ##x## and ##y## are, as mentioned earlier, the coordinates of the asteroid relative to the Earth which is assumed to be at the origin. The acceleration is obtained by dividing out the mass of the asteroid. Splitting it into Cartesian components,$$a_x= -GM_E\frac{x}{(x^2+y^2)^{3/2}};~~a_y= -GM_E\frac{y}{(x^2+y^2)^{3/2}}$$The ##k## as defined by @haruspex is the common factor ##\frac{GM_E}{(x^2+y^2)^{3/2}}## that appears in both ##a_x## and ##a_y##. Note that when the acceleration is expressed in this fashion, it is directed towards the force center and gets updated after ##x## and ##y## are updated.

Homework Helper
Gold Member
2022 Award
Can you tell me more about the equations that you have posted?

With the asteroid at coordinates (x,y) relative to Earth as origin,
I.e., x=Cx-Px, y=Cy-Py:
The force of attraction is ##\frac{Gm_Em_A}{(x^2+y^2)}##
The acceleration has magnitude ##\frac{Gm_E}{(x^2+y^2)}##
To get the magnitude of the component of this in the x direction we multiply by ##\frac x{(x^2+y^2)^\frac 12}##. But we know the direction will be opposite to x, so we get the signed component as ##-\frac {xGm_E}{(x^2+y^2)^\frac 32}##
In the y direction, the same except y in the numerator instead of x, so it is convenient to calculate
##k=\frac{Gm_E}{(x^2+y^2)^\frac 32}##
##\dot x-=kx\Delta t##
##\dot y-=ky\Delta t##.