c-----------------------------------------------------------
      program nbody
c-----------------------------------------------------------
c     This program (adapted from deut.f) uses LSODA to
c     integrate the equations of motion for n particles
c     interacting gravitationally.  The maximum number of
c     particles is set with nbdmax.
c-----------------------------------------------------------
c Assignments
      implicit none
      character*5 cdnm
      parameter (cdnm='nbody')
      integer iargc,indlnb,i4arg
      real*8 r8arg
      real*8 r8_never
      parameter (r8_never=-1.0d-60)
      integer nbdmax,neqmax
      parameter (nbdmax=5)
      parameter (neqmax=6*nbdmax)
      integer maxout
      parameter (maxout=10000)
      real*8 y0(neqmax)
      real*8 vtout(maxout),vyout(maxout,neqmax),work(maxout)
      real*8 tfin,tdel
      real*8 v1,v2,v3,v4,v5,v6,v7
      real*8 x_com,y_com,z_com,E_tot
      real*8 KE_tot,PE_tot,P_tot,J_tot
      real*8 m_tot,xi,yi,zi,xv,yv,zv,xj,yj,zj,r2
      real*8 x_mom,y_mom,z_mom
      real*8 x_ang,y_ang,z_ang
      real*8 m1,m2,m3,md,x1,x2,x3,xd,vd,rd,dd
      integer ir
      integer ntout,itout,ntout_succ
      integer nbd,ibd,jbd,neq,ieq,iof,jof
      logical ltrace,vtrace,errortrace,lsodatrace
      parameter (ltrace=.false.)
      parameter (vtrace=.false.)
      parameter (lsodatrace=.false.)
      character*32 argi
c-----------------------------------------------------------
c LSODA assignments
      external fcn,jac
      real*8 y(neqmax)
      real*8 tbgn,tend
      integer itol
      real*8 rtol,atol
      integer itask,istate,iopt
      integer lrw
      parameter (lrw=22+neqmax*16)
      real*8 rwork(lrw)
      integer liw
      parameter (liw=20+neqmax)
      integer iwork(liw)
      integer jt
      real*8 tol
      real*8 default_tol
      parameter (default_tol=1.0d-6)
      include 'fcn.inc'
      G = 1.0d0
c     G = 6.672d-6
c-----------------------------------------------------------
c Parse input (final t, delta t, optional tolerance,trace)
      if(iargc().lt.2) goto 99
      tfin = r8arg(1,r8_never)
      tdel = r8arg(2,r8_never)
      tol = r8arg(3,default_tol)
      errortrace = .false.
      if(iargc().ge.3) call getarg(3,argi)
      if(argi.eq.'trace') errortrace = .true.
      if(iargc().ge.4) call getarg(4,argi)
      if(argi.eq.'trace') errortrace = .true.
      if(tfin.eq.r8_never.or.tdel.eq.r8_never) goto 99
      if(ltrace) write(0,*) tfin,tdel,tol,errortrace
      if(vtrace) then
         v1 = sqrt(1.0d0*0.2d0/2.0d0**2)
         v2 = sqrt(9.0d0*1.8d0/2.0d0**2)
c        write(0,*) 'Two-body orbit @-.2,1.8 ',v1,v2
         m1 = 9.0d0
         m2 = 1.0d0
         m3 = 0.001d0
         md = m2 + m3
         rd = tfin
         dd = tdel
         x1 = -md*rd/(m1+md)
         xd = x1 + rd
         v1 = -sqrt(-G*md*x1/rd**2)
         vd = +sqrt(+G*m1*xd/rd**2)
         x2 = m3*dd/md
         x3 = x2 - dd
         v2 = -sqrt(+G*m3*x2/dd**2)
         v3 = +sqrt(-G*m2*x3/dd**2)
         x2 = xd + x2
         x3 = xd + x3
         v2 = vd + v2
         v3 = vd + v3
         write(0,*) '9.000 ',x1,' 0 0 0 ',v1,' 0'
         write(0,*) '1.000 ',x2,' 0 0 0 ',v2,' 0'
         write(0,*) '0.001 ',x3,' 0 0 0 ',v3,' 0'
      endif
c-----------------------------------------------------------
c Determine integration time steps
      if(tdel.le.0.0d0.or.tdel.gt.tfin) then
         write(0,*) cdnm//': delta t range (0<dt<=tfinal)'
         stop
      endif
      ntout = 1 + tfin/tdel
      if(ntout.gt.maxout) then
         write(0,*) cdnm//': delta t too small'
         stop
      endif
      vtout(1) = 0.0d0
      do itout=2,ntout
         vtout(itout) = vtout(itout-1) + tdel
      enddo
c-----------------------------------------------------------
c Read masses and initial position/velocity from stdin
      nbd = 0
5     continue
      read(5,*,iostat=ir,end=6) v1,v2,v3,v4,v5,v6,v7
      if(ir.eq.0) then
         if(nbd.eq.nbdmax) then
            write(0,*) cdnm//': max bodies input ',nbdmax
            goto 6
         else 
            nbd = nbd + 1
            iof = (nbd-1)*6
            M(nbd) = v1
            y0(iof+1) = v2
            y0(iof+2) = v5
            y0(iof+3) = v3
            y0(iof+4) = v6
            y0(iof+5) = v4
            y0(iof+6) = v7
         endif
      endif
      goto 5
6     continue
      if(nbd.le.0) then
         write(0,*) cdnm//': no input read from stdin'
         stop
      else
         neq = 6*nbd
      endif
c-----------------------------------------------------------
c Set LSODA parameters and initialize solution
      itol = 1
      rtol = tol
      atol = tol
      itask = 1 
      iopt = 0
      jt = 2
      do ieq=1,neq
         y(ieq) = y0(ieq)
         vyout(1,ieq) = y0(ieq)
      enddo
c-----------------------------------------------------------
c Perform integration
      do itout=2,ntout
         istate = 1
         tbgn = vtout(itout-1)
         tend = vtout(itout)
         call lsoda(fcn,neq,y,tbgn,tend,
     &              itol,rtol,atol,itask,    
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)
         if(lsodatrace) then
            write(0,10) cdnm, itout, ntout, vtout(itout), 
     &      vtout(itout+1)
10          format(' ',a,': Step ',i4,'  t = ',1pe10.3,
     &      ' .. ',1pe10.3)
            write(0,*) cdnm//': lsoda returns ', istate
         endif
         if(istate.lt.0) then
            write(0,20) cdnm, istate, itout, ntout, 
     &                    vtout(itout-1), vtout(itout)
20          format(/' ',a,': Error return ',i2,
     &              ' from integrator LSODA.'/
     &              '       At output time ',i5,' of ',i5/
     &              '       Interval ',1pe11.3,' .. ',
     &             1pe11.3/)
            ntout_succ = itout - 1
            goto 30
         endif
         do ieq=1,neq
            vyout(itout,ieq) = y(ieq)
         enddo
      enddo
      ntout_succ = ntout
30    continue
c-----------------------------------------------------------
c Write tracing output to standard error
      if(errortrace) then
         m_tot = 0.0d0
         do ibd=1,nbd
            m_tot = m_tot + M(ibd)
         enddo
         do itout=1,ntout_succ
            x_com = 0.0d0
            y_com = 0.0d0
            z_com = 0.0d0
            E_tot = 0.0d0
            KE_tot = 0.0d0
            PE_tot = 0.0d0
            P_tot = 0.0d0
            J_tot = 0.0d0
            x_mom = 0.0d0
            y_mom = 0.0d0
            z_mom = 0.0d0
            x_ang = 0.0d0
            y_ang = 0.0d0
            z_ang = 0.0d0
            do ibd=1,nbd
               iof = (ibd-1)*6
               xi = vyout(itout,iof+1)
               yi = vyout(itout,iof+3)
               zi = vyout(itout,iof+5)
               xv = vyout(itout,iof+2)
               yv = vyout(itout,iof+4)
               zv = vyout(itout,iof+6)
               x_com = x_com + M(ibd)*xi
               y_com = y_com + M(ibd)*yi
               z_com = z_com + M(ibd)*zi
               KE_tot = KE_tot + 
     &                  0.5d0*M(ibd)*(xv**2+yv**2+zv**2)
               x_mom = x_mom + M(ibd)*xv
               y_mom = y_mom + M(ibd)*yv
               z_mom = z_mom + M(ibd)*zv
               if(ibd.gt.1) then
                  do jbd=1,ibd-1
                     jof = (jbd-1)*6
                     xj = vyout(itout,jof+1)
                     yj = vyout(itout,jof+3)
                     zj = vyout(itout,jof+5)
                     r2 = (xi-xj)**2+(yi-yj)**2+(zi-zj)**2
                     r2 = sqrt(r2)
                     if(r2.gt.0.0d0) then
                        PE_tot = PE_tot - G*M(ibd)*M(jbd)/r2
                     endif
                  enddo
               endif
            enddo
            if(m_tot.gt.0.0d0) then
               x_com = x_com/m_tot
               y_com = y_com/m_tot
               z_com = z_com/m_tot
            endif
            do ibd=1,nbd
               iof = (ibd-1)*6
               xi = vyout(itout,iof+1) - x_com
               yi = vyout(itout,iof+3) - y_com
               zi = vyout(itout,iof+5) - z_com
               xv = vyout(itout,iof+2)*M(ibd)
               yv = vyout(itout,iof+4)*M(ibd)
               zv = vyout(itout,iof+6)*M(ibd)
               x_ang = x_ang + yi*zv - zi*yv
               y_ang = y_ang + zi*xv - xi*zv
               z_ang = z_ang + xi*yv - yi*xv
            enddo
            E_tot = KE_tot + PE_tot
            P_tot = sqrt(x_mom**2+y_mom**2+z_mom**2)
            J_tot = sqrt(x_ang**2+y_ang**2+z_ang**2)
            write(0,'(1P,9E26.16)') vtout(itout),
c           write(0,'(1P,9E9.1)') vtout(itout),
     &         x_com,y_com,z_com,E_tot,
     &         KE_tot,PE_tot,P_tot,J_tot
         enddo
      endif
c-----------------------------------------------------------
c Write standard output and exit
      write(*,*) nbd
      do ibd=1,nbd
         write(*,'(1P,E25.16)') M(ibd)
      enddo
      do itout=1,ntout_succ
         write(*,'(1P,E25.16)') vtout(itout)
         do ibd=1,nbd
            iof = (ibd-1)*6
            write(*,'(1P,3E25.16)') 
     &      (vyout(itout,iof+ieq),ieq=1,5,2)
         enddo
      enddo
      stop
c-----------------------------------------------------------
c Usage exit
99    continue
      write(0,*) 'usage: '//cdnm//' <t final> <dt out> ',
     &           '[<tol> <trace>]'
      write(0,*) ' '
      write(0,*) '       Input read from standard input'
      stop
      end