## Fortran 90 N-body problem code not working, help please!

## Fortran 90 N-body problem code not working, help please!

(OP)

Hey guys,

as an assignment, we are coding up the N-body problem for the solar system. I'm running into some problems however.

We need to code up the forces from one body to another, F_ij, the force on i due to j. We have been told do not code it in matrix form, as for problems with a million bodies, this becomes a 3 * a million * a million matrix, which represents formiddable memory challenges. So we have to do it as a rank 1 array, i.e. F = [F_1x, F_1y, F_1z, F_2x, F_2y, F_2z, ...] where for example, F_1x is the total force on body 1 in the x-direction.

The force has the form F_ij = G * m_i * m_j * (r_i - r_j)/((abs(r_i - r_j)^3)) - just Newton's force law. Another condition to be satisfied is F_ij = - F_ji. Obviously F_ii = 0.

I'm running into a big problem here, as I have no idea how to do it. What I have done so far hasn't worked. My latest attempt is below. nDim = number of dimensions, so 3, and in this case we're only doing the sun and the planets up to Jupiter, so 6 bodies (nBody). r is a rank-2 array with the positions of each of the bodies, with dimension r(nDim,nBody), mass is a rank-1 array of the masses of the bodies, F is the force array that I am trying to solve for. I have also defined another variable, Force, which is an intermediate step which I don't believe I need but I'm really stuck. nVector is just nDim * nBody * 2, that was defined for a different purpose. Never mind it, shouldn't make a difference.

subroutine forceij(nDim,nBody,nVector,mass,r,F)

real(DP), dimension(:,:), intent(in) :: r(:,:)

real(DP), dimension(:), intent(in) :: mass

integer, intent(in) :: nBody, nDim, nVector

real(DP), dimension(:), intent(out), allocatable :: F(:)

integer :: i, j, k ! counting variables

real(DP), dimension(nDim,nBody) :: Force

allocate(F(nVector/2))

do i = 1, nBody-1

do j = i+1, nBody

Force(i,j) = (G * mass(i) * mass(j) * (r(i,i) - r(i,j))/((abs(r(i,i) - r(i,j)))**3))

F(k) = sum(Force(:,j))

end do

end do

end subroutine forceij

Any help/suggestions/tips and whatnot is appreciated, thanks!

as an assignment, we are coding up the N-body problem for the solar system. I'm running into some problems however.

We need to code up the forces from one body to another, F_ij, the force on i due to j. We have been told do not code it in matrix form, as for problems with a million bodies, this becomes a 3 * a million * a million matrix, which represents formiddable memory challenges. So we have to do it as a rank 1 array, i.e. F = [F_1x, F_1y, F_1z, F_2x, F_2y, F_2z, ...] where for example, F_1x is the total force on body 1 in the x-direction.

The force has the form F_ij = G * m_i * m_j * (r_i - r_j)/((abs(r_i - r_j)^3)) - just Newton's force law. Another condition to be satisfied is F_ij = - F_ji. Obviously F_ii = 0.

I'm running into a big problem here, as I have no idea how to do it. What I have done so far hasn't worked. My latest attempt is below. nDim = number of dimensions, so 3, and in this case we're only doing the sun and the planets up to Jupiter, so 6 bodies (nBody). r is a rank-2 array with the positions of each of the bodies, with dimension r(nDim,nBody), mass is a rank-1 array of the masses of the bodies, F is the force array that I am trying to solve for. I have also defined another variable, Force, which is an intermediate step which I don't believe I need but I'm really stuck. nVector is just nDim * nBody * 2, that was defined for a different purpose. Never mind it, shouldn't make a difference.

subroutine forceij(nDim,nBody,nVector,mass,r,F)

real(DP), dimension(:,:), intent(in) :: r(:,:)

real(DP), dimension(:), intent(in) :: mass

integer, intent(in) :: nBody, nDim, nVector

real(DP), dimension(:), intent(out), allocatable :: F(:)

integer :: i, j, k ! counting variables

real(DP), dimension(nDim,nBody) :: Force

allocate(F(nVector/2))

do i = 1, nBody-1

do j = i+1, nBody

Force(i,j) = (G * mass(i) * mass(j) * (r(i,i) - r(i,j))/((abs(r(i,i) - r(i,j)))**3))

F(k) = sum(Force(:,j))

end do

end do

end subroutine forceij

Any help/suggestions/tips and whatnot is appreciated, thanks!

## RE: Fortran 90 N-body problem code not working, help please!

If so, then you want to store this matrix as an 1-dimensional array?

In this case the relation between the elements of the Matrix

Mand the Vectorvcould be:M(i,j) = v(3*(i-1) + j)so M(1,1)= v(1), M(1,2) = v(2), .., M(3,2)=v(8), M(3,3)= v(9), ..etc

In this case I would first compute all forces into the array, e.g:

## CODE

and for further computation with forces I would create a function which for given values of i and j returns the value of the force from the array

## CODE

## RE: Fortran 90 N-body problem code not working, help please!