Help Gravity with 3D coordinations

[Sword7]

Hello folks,

I am new to this physics forum since I searched for that through Google and found this.

I am figuring out about gravity calculation with 3D coordinations (or vector). I wrote my small program to simulate gravity for flying around the planet Earth. I tried that but a spacecraft quickly went away from Earth instead of orbit. I searched for vector calculation for gravity through the Internet but no luck or not enough information. Does anyone have any complete information about 3D coordination/vector calculation for gravity?

I have two books called "Fundamental of Astrodynamics" and "Orbit Determination for Microcomputers". I still am figuring them out because I am new to vector calculations.

My program for 3D gravity simulation are:

double G = 6.67428e-11;
double M = 5.9742e+24;
double m = 100;
double px, py, pz;
double ax, ay, az;
double Fg, a, r;

px = 7000.0 * 1000;
py = 7000.0 * 1000;
pz = 7000.0 * 1000;

ax = 0.0;
ay = 14.0 * 1000;
az = 0.0;

for(idx = 0; idx < 10000; idx++) {
r = sqrt((px*px) + (py*py) + (pz*pz));
Fg = (G * M * m) / (r * r);
a = Fg / m;

ax += a; ay += a; az += a;
px -= ax; py -= ay; pz -= az;

printf("Spacecraft = %g,%g,%g: Fg = %g N (%g m/s)\n",
px, py, pz, Fg, a);
}

Thanks!
Sword7

 

[Nabeshin]

r = sqrt((px*px) + (py*py) + (pz*pz));
Fg = (G * M * m) / (r * r);
a = Fg / m;

ax += a; ay += a; az += a;

You calculate the full radius vector, which you use in your computation of the force vector, which is the full force vector (the one parallel to the radius vector). You cannot use this acceleration in the x y and z directions -- you have to decompose the accelerations (or forces) into components first before you calculate the new positions.

 

[Sword7]

Hello folks,

Ok, I read my books again and learned about unit vector for gravity formula. I changed my code and tried again. Note: center of Earth is (0, 0, 0).

Fg = -( GMn / r^2 ) * (r / |r|);

:
:
r = sqrt((px*px) + (py*py) + (pz*pz));
Fgx = -(G * M * m * px) / (r*r*r);
Fgy = -(G * M * m * py) / (r*r*r);
Fgz = -(G * M * m * pz) / (r*r*r);

ax = Fgx / m;
ay = Fgy / m;
az = Fgz / m;

px += ax; py += ay; pz += az;
:
:

Results:

Spacecraft accelerated around Earth then quickly went away. At near center, its peak was near speed of light! Then it just started slow down but never return back. I think that spacecraft fell down to surface too fast instead of staying in orbit. However, there is another problem. When one of px/py/pz became negative, one of ax/ay/az jumped. Does anyone know any solution? That's why I want to write my own space simulator.

Thanks again,
Sword7

 

[Nabeshin]

I think you're still calculating the component forces wrong, what you have is this:
$$F_{gx}=-\frac{GMm}{r^3}p_{x}$$
The force of gravity is given by:
$$F_{g}=-\frac{GMm}{r^3}\vec{r}$$
Or, equivilently,
$$F_{g}=-\frac{GMm}{r^2}\hat{r}$$

Your formula isn't directly in either of these forms, so I think what you should use is something like this:
$$F_{gx}=-\frac{GMm}{p_{x}^2}\hat{p_{x}}$$

Which should correctly yield the force in the x dimension (the unit vector isn't really necessary in your calculations though).

 

[Sword7]

Hello,

Hmm. I noticed that nice math fonts that you tried to write. How do I write equations by using math formula/symbols fonts?

Ok, I was looking for unit vector information through my books and Internet but not much information like r-> or r^. How convert position vectors to unit vector? Does anyone have any good information about unit vector because I am very new to vector math?

Thanks again,
Sword7

 

[HallsofIvy]

You can use "LaTex". Some people have software that converts "WYSIWYG" to LaTex form but it is not difficult to do it "by hand" (and I think it produces a cleaner file). Look at
https://www.physicsforums.com/showthread.php?t=8997

A "unit vector" is simply a vector of length 1. To convert any (non-zero) vector to a unit vector in the same direction, divide each component by the length of the vector.

 

[D H]

The force of gravity is given by:
$$F_{g}=-\frac{GMm}{r^3}\vec{r}$$

That is exactly what he is doing now.

Your formula isn't directly in either of these forms, so I think what you should use is something like this:
$$F_{gx}=-\frac{GMm}{p_{x}^2}\hat{p_{x}}$$

Don't do that!

======================================================

Sword7, acceleration is the second derivative of position. Acceleration changes velocity, not position. Velocity changes position. The problem is that you are not representing velocity. You need to represent both the position and velocity of the object.

There are many ways to propagate state over time. You are using the simplest such technique, the Euler method. Euler integration is never going to be particularly accurate. That said, there are several variants even on the Euler method for a second order differential equation. Basic Euler uses

$$\aligned {\boldsymbol a}(t_{n}) &= -\,\frac {GM}{||{\boldsymbol x}(t_{n})||^3}{\boldsymbol x}(t_{n}) \\ {\boldsymbol x}(t_{n+1}) &= {\boldsymbol x}(t_{n}) + {\boldsymbol v}(t_{n})\Delta t \\ {\boldsymbol v}(t_{n+1}) &= {\boldsymbol v}(t_{n}) + {\boldsymbol a}(t_{n})\Delta t \endaligned$$

In other words, calculate the acceleration based on the position, update the position vector using the current velocity vector, and update the velocity vector using the acceleration just calculated.

A slight variation called the symplectic Euler method or semi-implicit Euler method will give much better results than basic Euler:

$$\aligned {\boldsymbol a}(t_{n}) &= -\,\frac {GM}{||{\boldsymbol x}(t_{n})||^3}{\boldsymbol x}(t_{n}) \\ {\boldsymbol v}(t_{n+1}) &= {\boldsymbol v}(t_{n}) + {\boldsymbol a}(t_{n})\Delta t \\ {\boldsymbol x}(t_{n+1}) &= {\boldsymbol x}(t_{n}) + {\boldsymbol v}(t_{n+1})\Delta t \endaligned$$

In other words, calculate the acceleration based on the position, update the velocity vector using the acceleration just calculated, and update the position vector using the updated velocity vector.

There are much, much better techniques than either of these. Ask if you are interested.

 

[Nabeshin]

That is exactly what he is doing now.

Is it? It looks to me like he's multiplying by the component rather than the full radius vector. The equation I provided was for the force along the radius vector, not yet in components. The two are certainly not equivilent, but it may be the case that what he's doing turns out to be correct and I'm simply not seeing it.

Don't do that!

...Why not?

 

[D H]

Is it? It looks to me like he's multiplying by the component rather than the full radius vector.

As you noted, a correct expression for the force due to gravity is

$$F_{g}=-\frac{GMm}{r^3}\vec{r}$$[/QUOTE]

Denote the components of the position vector as $p_x, p_y, p_z[/i$. Then [itex]r^2 = {p_x}^2+{p_y}^2+{p_z}^2[/tex]. The x,y,z components of the force vector are thus

$$\aligned F_x &=-\frac{GMm}{r^3}p_x \\ F_y &=-\frac{GMm}{r^3}p_y \\ F_z &=-\frac{GMm}{r^3}p_z \endaligned$$

This is exactly what Sword7 did.

Your formula isn't directly in either of these forms, so I think what you should use is something like this:
$$F_{gx}=-\frac{GMm}{p_{x}^2}\hat{p_{x}}$$

Don't do that!

...Why not?

Because it's wrong.

Consider the case $p_x=p_y=0, p_z=r\ne0$. You're equation would yield infinite values for the x- and y- components of the force vector.

 

[Nabeshin]

As you noted, a correct expression for the force due to gravity is

$$F_{g}=-\frac{GMm}{r^3}\vec{r}$$

Denote the components of the position vector as $p_x, p_y, p_z[/i$. Then [itex]r^2 = {p_x}^2+{p_y}^2+{p_z}^2[/tex]. The x,y,z components of the force vector are thus

$$\aligned F_x &=-\frac{GMm}{r^3}p_x \\ F_y &=-\frac{GMm}{r^3}p_y \\ F_z &=-\frac{GMm}{r^3}p_z \endaligned$$

This is exactly what Sword7 did.Because it's wrong.

Consider the case $p_x=p_y=0, p_z=r\ne0$. You're equation would yield infinite values for the x- and y- components of the force vector.

Ok I see now.

Spacecraft accelerated around Earth then quickly went away. At near center, its peak was near speed of light! Then it just started slow down but never return back. I think that spacecraft fell down to surface too fast instead of staying in orbit. However, there is another problem. When one of px/py/pz became negative, one of ax/ay/az jumped. Does anyone know any solution? That's why I want to write my own space simulator.

This sounds to me like your step size is too large, and so the when the p value jumped suddenly from -500 to 500 you see a large jump in one of the a's. As far as your object going near the speed of light, what are your initial conditions? Because that just sounds wrong. (Do you have a condition that stops the program if r<Re?)

 

[D H]

Read the second part of post #7. He is not representing velocity.

 

[Nabeshin]

Read the second part of post #7. He is not representing velocity.

Wow how did I miss that. Sorry.

 

[D H]

Wow how did I miss that. Sorry.

'Cause its hidden. In the first post, he was representing velocity, but he called it ax, ay, az. The problem there was that he has updating the velocity incorrectly by using acceleration as a scale. In Sword7's next post (post #3), he corrected the calculation of acceleration as a vector rather than a scalar, but messed things up by getting rid of velocity.

 

[Sword7]

Hello again,

Good news! Huston, we finally got my spacecraft in orbit...

I finally resolved a problem with my gravity calculation after I fixed my code. I corrected initial velocity setting from 14km/s to 7km/s and finally got spacecraft in orbit instead of runaway. That was elliptic orbit. I now recongized that 14km/s is too fast for Earth orbit.

I recongized that I miswrote my equations above and now mean that 'Ax += (Fgx / m); ...', not 'Ax += a; ...'. I enhanced my code by adding Vx, Vy, and Vz variable for current velocity updating.

Also, I added crash detection (r < Re). I set all initial velocity to zero and ran it. It fell down to Earth surface like dropping a ball and my code printed "Crash!" and exit. It worked so well.

Everything worked so well. I now am working on n-body equations that are affecting orbit path and 3D animations right now. Later I will add orbital determination and atmosphere techniques for air-resistance (re-entry, launches, and orbit decay).

DH, I saw one of your postings said that much better techniques than that technique from my astrodynamics book. Yes, I am interested in that.

Thank for your help.
Sword7

 

[D H]

Here's a start: http://en.wikipedia.org/wiki/Numerical_ordinary_differential_equations.