×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

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

 Forum Search FAQs Links MVPs

## 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!

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

I'm not sure if I understood you right, because you write once about [F_1x, F_1y, F_1z, F_2x, F_2y, F_2z, ...] and then about F_ij=... Does this mean that you have n*3 matrix F_ij where i=1,..,n and j=1,2,3 (or x,y,z) ?
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 M and the Vector v could 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

subroutine set_force(force_array,...)
real, dimesnion(3*n) :: force_array
real :: f
...
do i=1,n
do j=1,3
! compute Force
f = (G * mass(i) * mass(j) * (r(i,i) - r(i,j))/((abs(r(i,i) - r(i,j)))**3))
! store Force into array
force_array(3*(i-1) + j) = f
end do
end do
end subroutine set_force 

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

real function get_force(i,j,...)
...
! return force from array
get_force = force_array(3*(i-1) + j)
end function get_force 

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

(OP)
Okay great, thanks for your help!

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

#### Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!