N体间力向量计算结果异常，请求问题排查与解决

我正在编写一个函数，利用公式F=GMm/r²计算N个天体间的作用力向量，从.txt文件读取质量和初始位置数据，并将结果存储在动态分配的数组中。目前返回结果不正确，怀疑问题出在两个天体间的半径（距离）计算环节。

文件格式说明 #

initial_positions.txt（无表头） #

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

masses.txt #

1
2
3
4
5

说明：质量为1的天体初始位置为(1,2,3)，质量为2的天体初始位置为(4,5,6)，以此类推。

代码实现 #

#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()和ReadInitialPositions()这三个函数运行正常，恳请帮忙排查问题！

内容的提问来源于stack exchange，提问作者user11080418