# gravity problem (probably precision stuff?)

Page 1 of 1

## 8 Replies - 5493 Views - Last Post: 28 March 2008 - 08:07 AMRate Topic: //<![CDATA[ rating = new ipb.rating( 'topic_rate_', { url: 'http://www.dreamincode.net/forums/index.php?app=forums&module=ajax&section=topics&do=rateTopic&t=44872&amp;s=4cdeafcb19fb9a355406326e6c524908&md5check=' + ipb.vars['secure_hash'], cur_rating: 0, rated: 0, allow_rate: 0, multi_rate: 1, show_rate_text: true } ); //]]>

### #1 taure

Reputation: 0
• Posts: 16
• Joined: 09-September 07

# gravity problem (probably precision stuff?)

Posted 02 March 2008 - 03:38 AM

Hello there, sorry that I interrupt you but I've a problem with a gravity simulation app. Supposing we've a fixed sun and a planet going round, I've inveted the following code:

```struct {
double m;
double x, y;
} sun;

struct {
double m;
double x, y;
double vx, vy;
double ax, ay;
} planet;

const double G=10;
double r;
```

and the app loop:

```	case WM_CREATE:
planet.m=1;
planet.x=35;
planet.y=135;
planet.vx=0;
planet.vy=1;

sun.m=10;
sun.x=135;
sun.y=135;

SetTimer(hwnd, 1, 1, NULL);
break;
case WM_TIMER:
hdc=GetDC(hwnd);
SetPixel(hdc, planet.x, planet.y, 1024);
ReleaseDC(hwnd, hdc);

planet.x+=planet.vx;
planet.y+=planet.vy;
planet.vx+=planet.ax;
planet.vy+=planet.ay;
r=sqrt((sun.x-planet.x)*(sun.x-planet.x)+(sun.y-planet.y)*(sun.y-planet.y));
planet.ax=G*sun.m*(sun.x-planet.x)/pow(r,3);
planet.ay=G*sun.m*(sun.y-planet.y)/pow(r,3);
break;

```

etc.

It though resulted in the planet moving round a spiral not an elliptic curve. Could anybody please help me fixing it?

This post has been edited by taure: 02 March 2008 - 08:05 AM

Is This A Good Question/Topic? 0

## Replies To: gravity problem (probably precision stuff?)

### #2 Trogdor

Reputation: 15
• Posts: 627
• Joined: 06-October 06

## Re: gravity problem (probably precision stuff?)

Posted 02 March 2008 - 05:14 AM

i am not sure without consulting the books, but here is my long-shot at it:
is this simplified newtonian physics you put into that formula ?
i understand that you expect the result to be an elliptic and instead you get a spiral, right?
Then you are probably right: you either made a typo in the formulae (though i cant find any mistakes at a quick glance) or there is a rounding problem, though i find it strange that it can be so big.
btw, you did not initialize planet.ax and planet.ay

Can anyone with a better (fresher) knowledge of classic mechanical laws take a look?

### #3 taure

Reputation: 0
• Posts: 16
• Joined: 09-September 07

## Re: gravity problem (probably precision stuff?)

Posted 02 March 2008 - 08:27 AM

I think the formula looks so odd because of long names of the structs. I wasn't though sure if I did right by simplifying it. I mean, the x and y components of acceleration should have been calculated by multiplying the whole "a" by appropriate trigonometric functions. I didn't like the idea of adding another variable for storing the angle. Instead, I multiplied Newton's classic "a=G*M/(r^2)" by (rx/r) - the ratio between vertical distance and the whole radius so I finally got "a=G*M*rx/(r^3)".

I'm afraid that setting ax and ay to 0 at the beginning didn't work.

### #4 Trogdor

Reputation: 15
• Posts: 627
• Joined: 06-October 06

## Re: gravity problem (probably precision stuff?)

Posted 02 March 2008 - 02:03 PM

ax and ay where probably automatically set to 0. though not pretty it would probably not have made a big difference.
Your math looks good to me.
It must be a rounding problem.
But as far as i know there should not be such large (even noticeable) deviations in floating point math.

What compiler are you using, and on what platform?
What is the precision of doubles in that compiler (should be in the manual somewhere)

One thing i would try is to get as much of the pow and sqrt out as possible.
since r=sqrt((sun.x-planet.x)*(sun.x-planet.x)+(sun.y-planet.y)*(sun.y-planet.y)); and you use it later in pow(r,3) you must be able to take a litle shortcut there along the lines of"
```r_sqr = (sun.x-planet.x)*(sun.x-planet.x)+(sun.y-planet.y)*(sun.y-planet.y);
r = sqrt(r_sqr);
and later use
r_sqr * r

```

but i think that if there is a rounding problem somewhere here, that is not going to fix it (unless its a really bad implementation of pow() ofcourse)

### #5 NickDMax

Reputation: 2254
• Posts: 9,245
• Joined: 18-February 07

## Re: gravity problem (probably precision stuff?)

Posted 02 March 2008 - 05:03 PM

I did this one once.

Well, the math looks ok I suppose (I didn't get indepth) but you realize that finding a nice orbital path is not as simple as picking random numbers and letting things go.

You have to choose an appropriate velocity for the planet. Too much and it will fly right by the sun, too little and it will fall into the sun. Try playing around with your velocity vector. if you continue to spiral inward, then increase it a little, if you start to fly off then back down a little (you should be able to bracket the solution pretty quickly once you have found both a fall and a fly away).

Also keep in mind that many times when it flys off it is just in a large orbit.

### #6 Trogdor

Reputation: 15
• Posts: 627
• Joined: 06-October 06

## Re: gravity problem (probably precision stuff?)

Posted 03 March 2008 - 07:35 AM

That is not true.
Unless you start with a velocity of 0, in which case it flys directly into the sun (and in this case straight through) it should always become an elliptical trajectory.

### #7 taure

Reputation: 0
• Posts: 16
• Joined: 09-September 07

## Re: gravity problem (probably precision stuff?)

Posted 03 March 2008 - 11:08 AM

The compiler is MinGW from DevC++. As Trogdor suggested I've got rid of pow() at all. The startup values should then result in a circular orbit of 100px radius, but still I get an unwinding spiral like this: a gif screenshot

### #8 Trogdor

Reputation: 15
• Posts: 627
• Joined: 06-October 06

## Re: gravity problem (probably precision stuff?)

Posted 04 March 2008 - 05:20 AM

Interesting and frustrating.
The deviation is so big and so regular!

Can it be that somewhere, perhaps in SetPixel, the doubles are cast to int? just rounded or truncated? and not by value but by reference?

If you could print out a list of values of all the steps for say the first 3 iterations you might be able to track down both rounding errors and implicit casting.

### #33 Ghost Walker

Reputation: 5
• Posts: 1
• Joined: 28-March 08

## Re: gravity problem (probably precision stuff?)

Posted 28 March 2008 - 08:07 AM

Hello,

this is a bit old topic, but I wanted to show what your problem was. The order of your computation is wrong. You have to compute first ax and ay, then velocity and then the position:

```	  r=sqrt((sun.x-planet.x)*(sun.x-planet.x)+(sun.y-planet.y)*(sun.y-planet.y));
planet.ax=G*sun.m*(sun.x-planet.x)/pow(r,3);
planet.ay=G*sun.m*(sun.y-planet.y)/pow(r,3);
planet.vx+=planet.ax;
planet.vy+=planet.ay;
planet.x+=planet.vx;
planet.y+=planet.vy;

```

That's all