Calculating force vectors between N bodies returns incorrect values?

Question

I'm trying to write a function that calculates the force vectors acting on a body from N bodies using F = GMm/r^2 (given the masses and initial positions in .txt files), and stores the values in a dynamically allocated array. I think the issue stems from the calculation of the radius (distance between 2 bodies), as the values returned are incorrect. The position vectors in the initial_positions.txt file are in this format (without the headers):

pos_x pos_y pos_z
1 2 3
4 5 6
7 8 9
1 2 3
4 5 6

The masses in the masses.txt file are in this format:

1
2
3
4
5

So that the body with mass 1 has initial position (1, 2, 3), the body with mass 2 has initial position (4, 5, 6) etc.

My code to read in the files and calculate the force vectors:

#include <stdio.h>
#include <stdlib.h>

int NumberOfBodies(void) //finds the number of bodies from masses.txt file.
    {
    char character;
    char previousCharacter;
    int numberOfBodies = 1;

    FILE *file = fopen("masses.txt", "r");

    if (file == NULL)
        {
        printf("\nUnable to access the 'masses.txt' file.\n");
        exit(1);
        }
    else
        {
        while ((character = fgetc(file)) != EOF)
            {
            if (character == '\n' && previousCharacter != '\n')
                {
                numberOfBodies++;
                }
            previousCharacter = character;
            }
        }

    fclose(file);

    return numberOfBodies;
    }

double *ReadMasses(int numberOfBodies) //reads masses.
    {
    int row;
    int line;
    double *masses = malloc(sizeof(double) * numberOfBodies);

    FILE *file = fopen("masses.txt", "r");

    if (file == NULL)
        {
        printf("\nUnable to access the 'masses.txt' file.\n");
        exit(1);
        }

    for (row = 0; row < numberOfBodies; row++)
        {
        line = fscanf(file, "%lf", &masses[row]);

        if (line == EOF)
            {
            break;
            }
        }

    fclose(file);

    return masses;
    }

double **ReadInitialPositions(int numberOfBodies) //reads initial positions.
    {
    int row;
    int scan;
    double **initialPositions = malloc(sizeof(double*) * numberOfBodies);

    for (row = 0; row < numberOfBodies; row++)
        {
        initialPositions[row] = malloc(sizeof(double) * 3); //hardcoded as we only consider x, y, and z components of position.
        }

    FILE *file = fopen("initial_positions.txt", "r");

    if (file == NULL)
        {
        printf("\nUnable to access the 'initial_positions.txt' file.\n");
        exit(1);
        }

    for (row = 0; row < numberOfBodies; row++)
        {
        scan = fscanf(file, "%lf %lf %lf", &initialPositions[row][0], &initialPositions[row][1], &initialPositions[row][2]);

        if (scan == EOF)
            {
            break;
            }
        }

    fclose(file);

    return initialPositions;
    }

double **CalculateForces(int numberOfBodies, double *masses, double **initialPositions) //function  to calculate force vectors.
    {
    int row;
    int column;
    int currentBody = 0;
    double radius;
    double gravitationalConstant = 6.6743;
    double **forces = malloc(sizeof(double*) * numberOfBodies);

    for (row = 0; row < numberOfBodies; row++)
        {
        forces[row] = malloc(sizeof(double) * 3);
        }

    for (row = 0; row < numberOfBodies; row++)
        {
        for (column = 0; column < 3; column++)
            {
            if (row != currentBody)
                {
                radius = (initialPositions[row][column] - initialPositions[row][currentBody]); //I suspect the issue stems from this line.
                forces[row][column] = (gravitationalConstant * masses[row] * masses[currentBody]) / (radius * radius);
                currentBody++;
                }
            else
                {
                forces[row][column] = 0;
                currentBody++;
                }
            }
        }

    for (row = 0; row < numberOfBodies; row++)
        {
        for (column = 0; column < 3; column++)
            {
            printf(" %lf", forces[row][column]); //Prints force vectors.
            }
        printf(" \n");
        }

    return forces;
    }

int main(void)
    {
    int numberOfBodies;
    double *masses;
    double **initialPositions;

    numberOfBodies = NumberOfBodies();
    masses = ReadMasses(numberOfBodies);
    initialPositions = ReadInitialPositions(numberOfBodies);

    CalculateForces(numberOfBodies, masses, initialPositions);

    return 0;
    }

NumberOfBodies(), ReadMasses(), and ReadInitialPositions() all seem to be working as intended. Thanks in advance!:)

Answer 1

This has all the earmarks of code written without a firm grasp of the math. Try writing out the math first.

GMm/r^2 gives the scalar force between two bodies. It acts along the direction vector between the bodies. Consequently, the scalar must be split into its vector components. Splitting is just just multiplying by the direction vector scaled to unit length.

More to the point, if you are computing the force between bodies a and b, then

Rab^2 = (xb - xa)^2 + (yb - ya)^2 + (zb - za)^ 2
Fab = ma mb / Rab^2

A unit vector from a to b has coordinates

Uabx = (xb - xa) / Rab
Uaby = (yb - ya) / Rab
Uabz = (zb - za) / Rab

Of course Rab = sqrt(Rab^2) computed above.

The three components of the force acting on a due to b are

Fabx = Uabx * Fab
Faby = Uaby * Fab
Fabz = Uabz * Fab

The force acting on b due to a is Fb = -Fa.

Working this out as an algorithm:

For each pair a,b of bodies
  Find r^2 the distance between a and b.
  Find Fab the scalar force between a and b.
  Find Ua the unit vector from a toward b.
  Find vector Fa by splitting Fab into three components using Ua.
  Add Fa into the total force acting on a.
  Add -Fa into the total force acting on b. 

A quick look at your code shows it can't be very close to correct. Ex: it computes three different values of r for a pair of bodies. There can be only one distance between two points in space!

Hint: To get all pairs of integers (ignoring order) in [0..n-1] that don't include self-pairs like (1, 1), the standard pattern is loops that look like this:

for (j = 1; j < n; j++)
  for (i = 0; i < j; i++)

This will iterate over i-j pairs in the order (0, 1), (0, 2), (1, 2), (0, 3), (1, 3), (2, 3), ...