1
\$\begingroup\$

I don't know if this question get her place in physics.stackoverflow forum. I prefer post here, with code etc.

So, I want simulate 2D movements (no rotation) of an asteroid in space affected by a planet's gravity, the problem is how integrate that ? Any better integration for performance ?

Here is the code in C:

struct vec2
{
  float x;
  float y;
};

struct body
{
  double mass; //kg
  struct vec2 prev_pos; //for semi-implicit euler method
  struct vec2 pos; //position
  struct vec2 v; //velocity
};

void move_asteroid(float dt) {
  //semi-implicit euler integration #1 -moving the object
  asteroid.pos.x += asteroid.v.x * dt;
  asteroid.pos.y += asteroid.v.y * dt;
  
  //START
  float x = planet.pos.x - asteroid.pos.x;
  float y = planet.pos.y - asteroid.pos.y;
  float d = x * x + y * y; //distance^2 between asteroid and center of planet
  
  //acceleration
  float ah = G * ((planet.mass * asteroid.mass) / d); //gravity
 
  //get components of gravity
  float rad = atan2(y, x);
  float ax = ah * cos(rad);
  float ay = ah * sin(rad);

  
  //semi-implicit euler integration #2 -update velocity
  asteroid.v.x += ax * dt;
  asteroid.v.y += ay * dt;
}

this is for 'fun' and sadly my knowledge is quite not good, I try to find with my vocabulary (not native english), if you have directions, ideas for better implementation (like study RK4 for correctness etc), it's welcome

\$\endgroup\$
8
  • 1
    \$\begingroup\$ You can replace your trigonometry with a vector normalization instead (ax = x * (ah / sqrt(d))) but otherwise this looks reasonable. Is there anything you're unsatisfied with about how this is performing in your game? \$\endgroup\$
    – DMGregory
    Jan 31 at 15:21
  • \$\begingroup\$ Actually this working well, (just need handle case where d=0). I started using a wrong method (movements looks good but values was totally wrong). So I want to be sure it's the right way todo. Vector normalization ? That was the thing I tried to do but with big misconception (like this ugly thing: ax = ah * (x*x/d)). This is too early but sqrt() isn't bottleneck ? simulation will take lot of object/planets. I keep it in a side, thank you! \$\endgroup\$
    – Mendes
    Jan 31 at 15:48
  • 1
    \$\begingroup\$ Re: bottlenecks, compare the costs of one square root and one division versus three trig functions: an arctangent, a cosine, and a sine. \$\endgroup\$
    – DMGregory
    Jan 31 at 15:51
  • 1
    \$\begingroup\$ Are you interested in actually simulating it step-by-step by integration, or are you open to analytical solutions? \$\endgroup\$
    – Sacha
    Jan 31 at 17:13
  • 2
    \$\begingroup\$ The "ah" you currently compute is the force, not the acceleration. To get the latter, remove the asteroid.mass (a more massive asteroid must not fall faster, the increased force is counterbalanced by the increased inertia, hence why the asteroid mass is irrelevant to get the gravitationnal acceleration) \$\endgroup\$
    – Keelhaul
    Feb 1 at 8:28
2
\$\begingroup\$

Following correction from @Keelhaul's comment and performance improvement from @DMGregory's comment.
Here is the right code:

struct vec2
{
  float x;
  float y;
};

struct body
{
  double mass; //kg
  struct vec2 prev_pos; //for semi-implicit euler method
  struct vec2 pos; //position
  struct vec2 v; //velocity
};

void move_asteroid(float dt) {
  //semi-implicit euler integration #1 -moving the object
  asteroid.pos.x += asteroid.v.x * dt;
  asteroid.pos.y += asteroid.v.y * dt;
  
  //START
  float x = planet.pos.x - asteroid.pos.x;
  float y = planet.pos.y - asteroid.pos.y;
  float d = x * x + y * y; //distance^2 between asteroid and center of planet
  
  //acceleration
  /*
    F = G * ((objectA.mass * objectB.mass)/d^2) //gravitational force
    a_objectA = F/objectA.mass //acceleration applied to objectA
  */
  
  float ah = G * (planet.mass / d); //acceleration applied to asteroid
 
  //get x,y components of acceleration
  float norm = ah / sqrt(d);
  float ax = x * norm;
  float ay = y * norm;

  
  //semi-implicit euler integration #2 -update velocity
  asteroid.v.x += ax * dt;
  asteroid.v.y += ay * dt;
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.