My Gravitational Pull Algorithm behaves Oddly in Certain Situations

Question

While trying to create my own physics engine (don't try persuading me not to), I decided to create a class for each pixel, called Particle, this system has an x and a y, and a x and y velocity, as shown below. Unfortunately, the code for calculateGravitationalVelocity doesn't abide by the laws of physics in certain situations. For example, if the x of the particle and the x of the other particle is the same, the particle will fall towards the object realistically, but when the particle gets too close, it pings off towards the positive x. I am only going to include the class source code, but I can include the source code of the other file, though it is partly written in SFML

Particle.cpp:

#include <iostream>
#include <string>
#include <math.h>

class Particle
{
    private:
        //Coords:
        double x, y;
        //Velocities:
        double xVelocity = 0;
        double yVelocity = 0;
        //Material:
        std::string material = "Generic";
        //Mass:
        double mass = 0;
    public:
        //Coords:
        void setCoords(double, double);
        float getCoords(char);
        //Velocities:
        void giveVelocity(char, float);
        void setVelocity(char, float);
        float getVelocity(char);
        //Gravitational Velocity:
        void calculateGravitationalVelocity(Particle);
        //Material:
        void setMaterial(std::string);
        std::string getMaterial();
        //Mass:
        void setMass(double);
        double getMass();
        //Update:
        void update();
};

//Coords:
void Particle::setCoords(double newX, double newY)
{
    x = newX;
    y = newY;
}
float Particle::getCoords(char axis)
{
    if (axis == 'x')
    {
        //return floor(x);
        return x;
    }
    else if (axis == 'y')
    {
        //return floor(y);
        return y;
    }
}

//Velocities:
void Particle::giveVelocity(char axis, float addedVelocity)
{
    if (axis == 'x') {xVelocity = xVelocity + addedVelocity;}
    else if (axis == 'y') {yVelocity = yVelocity + addedVelocity;}
}
void Particle::setVelocity(char axis, float newVelocity)
{
    if (axis == 'x') {xVelocity = newVelocity;}
    else if (axis == 'y') {yVelocity = newVelocity;}
}
float Particle::getVelocity(char axis)
{
    if (axis == 'x') {return xVelocity;}//floor(xVelocity);}
    else if (axis == 'y') {return xVelocity;}//floor(yVelocity);}
}

//Gravitational Velocity (Where the problems probably are):
void Particle::calculateGravitationalVelocity(Particle distantParticle)
{
    //Physics constants:
    const double pi = 3.14159265359; //Pi
    const double G = 0.00000000006673; //Gravitational Constant (or Big G)

    //Big Triangle Trigonometry:
    //Get coords of moving particle:
    double x1 = x;
    double y1 = y;
    //Get coords of particle with gravity:
    double x2 = distantParticle.getCoords('x');
    double y2 = distantParticle.getCoords('y');
    if (x1 != x2)
    {
        //Work out the angle:
        double A = atan((y2 - y1) / (x2 - x1)) * 180 / pi;
        //Remove the minus sign:
        A = fabs(A);

        //Small Triangle Trigonometry:
        //Work out the hypotenuse of the big triangle:
        double hyp = sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
        //Work out the gravitational field (hyp of small triangle):
        long double gravitationalField = G * (distantParticle.getMass() / pow(hyp, 2));
        //For testing purposes:
        //std::cout << "X: " << (cos(A) * gravitationalField) / 1000 << std::endl;
        //std::cout << "Y: " << (sin(A) * gravitationalField) / 1000 << std::endl;
        //Work out the X velocity:
        xVelocity = xVelocity + (cos(A) * gravitationalField) / 1000;
        //Work out the Y velocity:
        yVelocity = yVelocity + (sin(A) * gravitationalField) / 1000;
    }
    else
    {
        //Work out the hypotenuse of the big triangle:
        double hyp = sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
        //Work out the gravitational field (hyp of small triangle):
        long double gravitationalField = G * (distantParticle.getMass() / pow(hyp, 2));

        yVelocity = yVelocity + gravitationalField / 1000;
    }
}

//Material:
void Particle::setMaterial(std::string newMaterialType)
{
    material = newMaterialType;
}
std::string Particle::getMaterial()
{
    return material;
}

//Mass:
void Particle::setMass(double newMass)
{
    mass = newMass;
}
double Particle::getMass()
{
    return mass;
}

//Update:
void Particle::update()
{
    x = x + xVelocity;
    y = y + yVelocity;
}

I am sorry for the very open question, and it probably goes against the rules somewhere, but I couldn't find it. The code for working out mostly uses a two triangles to make a x and y velocity. Here is an image of what I was hoping the code would do as a triangle (sorry it doesn't look great, but I like using a whiteboard):

Answer 1

You don't need to perform any trigonometric calculation.

    ...

    //Get coords of particle with gravity:
    double x2 = distantParticle.getCoords('x');
    double y2 = distantParticle.getCoords('y');

    // Get difference vector
    double rx = x1 - x2;
    double ry = y1 - y2;

    // square of distance
    double r2 = rx * rx + ry * ry;
    // distance
    double r = sqrt (r2);

    if (r != 0) {
        // normalize difference vector
        double ux = rx / r;
        double uy = ry / r;

        // acceleration of gravity
        double a = - G * distantParticle.getMass() / r2;

        xVelocity += a * ux / 1000;
        yVelocity += a * uy / 1000;
    }
}