What's wrong with my math in this function to update the position of a planet near a star?

Question

Initially the code seems to work as the planet curves toward the star, but then as it should either get pulled into the star or make an orbit, it just gets pushed away in the opposite direction. What am I doing wrong? It seems like the change in behavior happens when the sign of distance changes.

xPos += t * velocity/10000 * cos(direction / 180 * 3.14);
yPos += t * velocity/10000 * sin(direction / 180 * 3.14);

double gravity = 1000;
double starX = 1920/2;
double starY = 1080/2;

double deltaX = xPos - starX;
double deltaY = yPos - starY;
double distance = sqrt(pow(deltaX, 2) + pow(deltaY, 2));
int modifier = 1;
if (xPos > starX){
    modifier = -1;
}
double angle = atan(deltaY / deltaX) * 180 / 3.14;

xPos += t * gravity / pow(distance,2) * modifier * cos(angle / 180 * 3.14);
yPos += t * gravity / pow(distance,2) * modifier * sin(angle / 180 * 3.14);

xPos and yPos is the current position of the planet. t is the time since the last update. velocity and direction are the initial velocity and angle that the planet begins at. This force never changes since we are in space. gravity is the gravitational constant. starX and starY is the position of the star and distance is the distance between the star and the planet. angle is the angle between the star and the planet (I have a feeling that this is what is causing the unwanted behavior.)

I apologize if this is the wrong place for this question but I figured it was the right place since this is more of a math question than a programming one.

Answer 1

The following should work at least approximately. The granularity of time may cause the planet not to follow a closed orbit. To improve precision you should use a smaller $\Delta t$ and/or learn something about the numerical analysis of handling this type of differential equations more accurately.

Initialize variables (run tests to find combos that work for DeltaT):

DeltaT= ; (=time between frames)
xPos= ; (= x-coordinate of the planet at the beginning)
yPos= ; (= y-coordinate of the planet at the beginning)
xSpeed= ;  (= x-component of the planet's speed at the beginning)
ySpeed= ;  (= y-component of the planet's speed at the beginning)

Then a Do - Until loop in pseudocode (comments on lines starting with a %-sign):

Do
% update the planet's position
xPos=xPos+DeltaT*xSpeed;
yPos=yPos+DeltaT*ySpeed;
% calculate the current relative position to the star and the acceleration due to gravity
DeltaX=xPos-xStar;
DeltaY=yPos-yStar;
distance=sqrt(pow(DeltaX,2)+pow(DeltaY,2));
GravityMagnitude=GravityConst*StarMass/pow(distance,2);
% update the planet's speed
xSpeed=xSpeed-DeltaT*GravityMagnitude*DeltaX/distance;
ySpeed=ySpeed-DeltaT*GravityMagnitude*DeltaY/distance;
Until EndOfSimulation;

The acceleration caused by gravity is a vector. It's magnitude follows the inverse square law, but we need to break that down to $x/y$-components as well. This is exactly what the multipliers DeltaX/distance (and same for $y$) achieve.

Once you are happy with this it's time for some additions! Introduce several planets/stars. All moving, all being affected by the gravity of all the others. Trust me, even with three celestial bodies (of comparable masses) strange things will happen. Enjoy!