c This file implements a N-body problem in 2D. The body with the
c mass m_i interacts a body m_j by the Newton's gravity law:
c
c                         m_i m_j
c    grav. force = gamma ------  , r=((x_i-x_j)**2+(y_i-y_j)**2)**1/2
c                          r**2
c
c Hence we have the equations of motion:
c
c     r''_i = sum (grav. force (j), j=1,..,N, j<>i) / m_i  (**)
c
c
c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   F 2 N D
c-----!-----------------------------------------------------------------
c
c Compute the right hand side of the ODE (**) :
c
      subroutine f2nd(n,t,u,fw)
      implicit none
      integer n
      real*8 t,u(n),fw(n)

      integer mmax
      parameter (mmax=10000)
      real*8 fx(mmax),fy(mmax),fpx(mmax),fpy(mmax),mm(mmax)
      common /fdata/fx,fy,fpx,fpy,mm

      integer i,j,nh
      real*8 r,x,y,sx,sy
      real*8 xj,yj
	  real*8 gamma

	  data gamma /6.6720d-11/

      nh=n/2

      do i=1,nh
         sx=0d0
         sy=0d0
         x=u(2*i-1)
         y=u(2*i)
         do j=1,nh
            if (j.ne.i) then
               xj=u(2*j-1)
               yj=u(2*j)
               r=((x-xj)**2+(y-yj)**2)**1.5d0
c               if (r.lt.1d-3) r=1d-3
               sx=sx+gamma*mm(j)*(xj-x)/r
               sy=sy+gamma*mm(j)*(yj-y)/r
            end if
         end do
         fw(2*i-1)=sx
         fw(2*i)=sy
      end do
      return
      end


c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   P R O B
c-----!-----------------------------------------------------------------
c
c Read the initial values from a file.
c 
c
      subroutine prob(name,n,t0,te)
      implicit none
      character*12 name
      integer n,mmax,i,nh

      real*8 t0,te

      parameter (mmax=10000)
      real*8 fx(mmax),fy(mmax),fpx(mmax),fpy(mmax),mm(mmax)
      common /fdata/fx,fy,fpx,fpy,mm

      character*12 cnamestr
      common /cname/ cnamestr

      open (11,file='moon.dat',err=100,status='old')
      read (11,*)
      read (11,*)
      read (11,*) name,nh,t0,te


      n=2*nh

      cnamestr=name

      do i=1,nh
         read(11,*,err=200) mm(i),fx(i),fy(i),fpx(i),fpy(i)
      end do

      close (11)
      return

 100  continue
      stop "star::prob, can not open file"

 200  continue
      stop "error reading star data!"
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   I N I T
c-----!-----------------------------------------------------------------
c
c Copy inital values form input arrays to up and u.
c
      subroutine init(n,u,up)
      implicit none
      integer n,nh
      real*8 u(*),up(*),i

      integer mmax
      parameter (mmax=10000)
      real*8 fx(mmax),fy(mmax),fpx(mmax),fpy(mmax),mm(mmax)
      common /fdata/fx,fy,fpx,fpy,mm

      nh=n/2

      do i=1,nh
         u(2*i-1)=fx(i)
         u(2*i)=fy(i)
         up(2*i-1)=fpx(i)
         up(2*i)=fpy(i)
      end do

      return
      end


c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   S O L U T
c-----!-----------------------------------------------------------------
c
c Read the reference solution from a file.
c
      subroutine solut(n,t,u)
      implicit none

      integer n,i,nh
      real*8 t,u(*)

      character*12 cnamestr
      common /cname/cnamestr

      open (11,file='moon_solut',err=200,status='old')
      read (11,*)
      read (11,*)


      do i=1,n
         read (11,*) u(i)
      end do
      close (11)
      return

 200  continue
      stop "can not read reference solution"

      return
      end