      implicit NONE

      integer NEQ,MAX
      PARAMETER (NEQ=7,MAX=100000)

      real*8 y(NEQ)

      integer j,ijob,itp,ntp
      character*80 fbin
      character*100 fhyd
      character*5 ext
      character*1 extph

      integer i,itpbuf,iplbuf,ipl(MAX),nenc
      real*8 tbuf,abuf,qbuf,rbuf
      real*8 time(MAX),a(MAX),e(MAX)

      real*8 au,year,dpi,rhill(10),mass(10)

      integer idum
      real*8 xi,ran3

      integer eflag,iph,ialpha,ifailfl
      real*8 diambin,mbin,abin,ebin,incbin,qbin
      real*8 capombin,omegabin,capmbin,ri,phii
      real*8 t,tf,eps,h1,hmin,hdid,dt

      integer ia,idiam,na,ndiam
      real*8 abin0

      integer ihyd(MAX),nhyd
      real*8 thyd(MAX),whyd(MAX),totw

      real*8 telim

      real*8 gm,gmp,ahyp,ehyp
      common /par1/ gm,gmp,ahyp,ehyp

c --------------------------------------------------------------

c     elimination info
      telim=6089164.5d0

c     constants
      au=149597870.7d0
      year=365.25d0*86400.d0
      dpi=acos(-1.d0)

c     Total number of selected particles
      ntp=1000000

c     two loops
      ndiam=1
      na=100
      do idiam=1,ndiam
         diambin=2.0       !+1.0*(float(idiam)-1.0)/float(ndiam-1)
         diambin=10.d0**diambin
         
         do ia=91,100
            abin0=2.0+3.0*(float(ia)-1.0)/float(na-1)     
            abin0=10.0**abin0
            abin0=abin0/au

c     Define binary parameters
      mbin=4.d0/3.d0*dpi*(diambin*1.d3/2.d0)**3*1000.d0*6.67408d-11 
     &     * year**2/(1000.d0*au)**3 

c      write(*,*)'mbin=',mbin

      gm=mbin+mbin
      
      do j=1,ntp                ! loop over particles  

         eflag=0                ! expecting to survive

         if(abin0*au.le.diambin) then
            eflag=2
            abin=abin0
            ebin=0.d0
            goto 301
         end if

c     Define the binary orbit
      abin=abin0   
      ebin=0.d0
      incbin=acos(2.d0*ran3(idum)-1.d0)
      capombin=ran3(idum)*2.d0*dpi
      omegabin=ran3(idum)*2.d0*dpi
      capmbin=ran3(idum)*2.d0*dpi

c     loop over phases, only one phase for Trojans
      do iph=1,1

c     get the particle job and id
         ijob=(j-1)/500+1
         itp=j-(ijob-1)*500 

c         write(*,*)j,itp,ijob   ! check on this
         
c     Set up extension for the input file
         write(ext,78) ijob
 78      format(i5.5)

c     Read the weights for each specific particle
      fhyd = 'caseI_trojan/hydra_trojan.dat_'//ext
      i=1
      open(1,file=fhyd)
 400  continue
      read(1,*,err=401,end=401)thyd(i),ihyd(i),whyd(i)
      if(ihyd(i).eq.itp) i=i+1
      goto 400
 401  continue
      nhyd=i-1   ! this is the number of outputs in the hydra file for this specific particle
      close(1)

      if(nhyd.eq.0) goto 666    ! ignore this particle if it is not in hydra file

c     compute the total weight for this particle
      totw=0.d0
      do i=1,nhyd
         totw=totw+whyd(i)
      end do

c      write(*,*)ijob,itp,nhyd,totw      

c     Read the binary file for specific particle
      fbin = 'caseI_trojan/binary.dat_'//ext  
      i=1
      open(1,file=fbin)
 200  continue
      read(1,*,err=201,end=201)tbuf,itpbuf,iplbuf,abuf,qbuf,rbuf
      if(itpbuf.eq.itp.and.tbuf.lt.thyd(nhyd)) then    
         time(i)=tbuf
         ipl(i)=iplbuf
         a(i)=abuf
         e(i)=1.d0+qbuf/abuf
c         write(*,*)itpbuf,iplbuf,tbuf,qbuf
         i=i+1
      end if
      goto 200
 201  continue
      nenc=i-1
      close(1)

c      goto 666

c     Big loop over all encounters
      do i=1,nenc

         if(time(i).lt.telim) then
            rhill(2)=5.21*(3.77d-2/3.d0/39.5d0)**0.333333d0
            rhill(3)=9.55*(1.13d-2/3.d0/39.5d0)**0.333333d0
            rhill(4)=13.0*(2.00d-3/3.d0/39.5d0)**0.333333d0
            rhill(5)=19.22*(2.00d-3/3.d0/39.5d0)**0.333333d0
            rhill(6)=30.11*(2.00d-3/3.d0/39.5d0)**0.333333d0
            mass(2)=3.769225110089709d-002
            mass(3)=1.128589982365781d-002
            mass(4)=0.200000000000000d-003
            mass(5)=0.200000000000000d-003
            mass(6)=0.200000000000000d-003
         else
            rhill(2)=5.21*(3.77d-2/3.d0/39.5d0)**0.333333d0
            rhill(3)=9.55*(1.13d-2/3.d0/39.5d0)**0.333333d0
            rhill(4)=19.22*(2.00d-3/3.d0/39.5d0)**0.333333d0
            rhill(5)=30.11*(2.00d-3/3.d0/39.5d0)**0.333333d0
            rhill(6)= 30.11*(2.00d-3/3.d0/39.5d0)**0.333333d0           
            mass(2)=3.769225110089709d-002
            mass(3)=1.128589982365781d-002
            mass(4)=0.200000000000000d-003
            mass(5)=0.200000000000000d-003
            mass(6)=0.200000000000000d-003
         end if

      gmp=mass(ipl(i))
      ialpha=-1
      call orbel_el2xv(gm,ialpha,abin,ebin,incbin,
     &     capombin,omegabin,capmbin,
     &     y(1),y(2),y(3),y(4),y(5),y(6))

      if(e(i).gt.1.d0) then 
         ri=2.d0*rhill(ipl(i))  
         phii=acos((a(i)*(e(i)*e(i)-1.d0)/ri-1.d0)/e(i))
         y(7)=-phii
      else
         ri=2.d0*rhill(ipl(i))  
         if(ri.gt.2.d0*a(i)) then
            phii=dpi
            y(7)=-phii
         else
            phii=acos((a(i)*(1.d0-e(i)*e(i))/ri-1.d0)/e(i))
            y(7)=-phii
         end if
      end if

      ahyp=a(i)
      ehyp=e(i)

c     Integrate over encounters
      t=0.d0
      dt=1.d-2
      hmin=1.d-10
      hdid=1.d-4
      eps=1.d-8

 10   continue
       tf=t+dt
       h1=hdid

c       write(*,*)t

c .... here is the propagator (t -> tf) ............................
       call odeint(y,neq,t,tf,eps,h1,hmin,hdid,ifailfl)
       if (ifailfl.ne.0) then
          eflag=3
          goto 301
       end if
c .... here is the propagator ......................................

       call orbel_xv2el(y(1),y(2),y(3),y(4),y(5),y(6),gm,ialpha,
     &           abin,ebin,incbin,capombin,omegabin,capmbin)
c       write(*,*)tf,abin,ebin,incbin,capombin,omegabin,capmbin

c       write(*,*)tf,y(7),abin,ebin

       t=tf

      call orbel_xv2el(y(1),y(2),y(3),y(4),y(5),y(6),gm,ialpha,
     &           abin,ebin,incbin,capombin,omegabin,capmbin)

      if(ialpha.ge.0) then
         eflag=1
         goto 301               ! binary dissociated
      end if

      qbin=abin*(1.d0-ebin)*au
      if(qbin.le.diambin) then
         eflag=2
         goto 301               ! collision
      end if

      if (y(7).lt.phii) goto 10

c-----..............................................................

      incbin=acos(2.d0*ran3(idum)-1.d0)
      capombin=ran3(idum)*2.d0*dpi
      omegabin=ran3(idum)*2.d0*dpi
      capmbin=ran3(idum)*2.d0*dpi

c      write(*,*)abin*au,ebin,a(i)*(e(i)-1.d0)

      end do  ! over encounters i

      end do  ! over phases iph

c     Output
 301  continue
      write(*,8848)diambin,abin0*au,ijob,itp,abin*au,ebin,eflag,totw
 8848 format(f6.1,1x,f9.1,1x,i4,1x,i3,1x,g11.4,1x,g11.4,1x,i1,1x,f9.1)


 666  continue

      end do  ! over particles j

      end do ! over abin0 ia
      end do ! over diambin idiam

c
      stop
      end

c ====================================================================
      SUBROUTINE bs_derivs(nv,x,y,dy)
      implicit none
      integer nv,NMAX
      PARAMETER (NMAX=7)
      real*8 x,y(NMAX),dy(NMAX)

      real*8 r,rr,xx,yy,r1x,r1y,r1z,r2x,r2y,r2z,r1,r2

      real*8 gm,gmp,ahyp,ehyp
      common /par1/ gm,gmp,ahyp,ehyp

      integer i

c ---------------------------------------------------------------
      
      dy(1)=y(4)
      dy(2)=y(5)
      dy(3)=y(6)
      r=sqrt(y(1)*y(1)+y(2)*y(2)+y(3)*y(3))
       if(ehyp.gt.1.d0) then
         rr=ahyp*(ehyp*ehyp-1.d0)/(1.d0+ehyp*cos(y(7)))
      else
         rr=ahyp*(1.d0-ehyp*ehyp)/(1.d0+ehyp*cos(y(7)))
      end if
      xx=rr*cos(y(7))
      yy=rr*sin(y(7))
      r1x=xx-0.5d0*y(1)   ! ASSUMED THAT M1=M2 !!!!!!!!!!!!!!!!!
      r1y=yy-0.5d0*y(2)
      r1z=-0.5*y(3)
      r2x=xx+0.5d0*y(1)   ! ASSUMED THAT M1=M2 !!!!!!!!!!!!!!!!!
      r2y=yy+0.5d0*y(2)
      r2z=0.5*y(3)
      r1=sqrt(r1x*r1x+r1y*r1y+r1z*r1z)
      r2=sqrt(r2x*r2x+r2y*r2y+r2z*r2z)

      dy(4)=-gm*y(1)/(r*r*r) - gmp*(r2x/(r2*r2*r2)-r1x/(r1*r1*r1)) 
      dy(5)=-gm*y(2)/(r*r*r) - gmp*(r2y/(r2*r2*r2)-r1y/(r1*r1*r1)) 
      dy(6)=-gm*y(3)/(r*r*r) - gmp*(r2z/(r2*r2*r2)-r1z/(r1*r1*r1)) 

      if(ehyp.gt.1.d0) then
         dy(7)=sqrt(gmp)/((ahyp*(ehyp*ehyp-1.d0))**1.5d0) *
     &        (1.d0+ehyp*cos(y(7)))**2
      else
         dy(7)=sqrt(gmp)/((ahyp*(1.d0-ehyp*ehyp))**1.5d0) *
     &        (1.d0+ehyp*cos(y(7)))**2
      end if

c      write(*,*)'time=',x
c      write(*,*)(y(i),i=1,7)
c      write(*,*)(dy(i),i=1,7)

c
      return
      end
c ====================================================================


c ====================================================================
      SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,hdid,ifailfl)
c This is the driver that make the BS scheme go through
c the required step x1 --> x2 in time; outputs the resulting
c state vector in ystart
      implicit none
      INTEGER i,nvar,NMAX,ifailfl
      REAL*8 eps,h1,hmin,x1,x2,TINY
      PARAMETER (NMAX=7,TINY=1.d-30)
      REAL*8 h,hdid,hnext,x,dydx(NMAX),y(NMAX),yscal(NMAX)
      REAL*8 ystart(NMAX)
c -------------------------------------------------------------------
      x=x1
      h=sign(h1,x2-x1)
      do i=1,nvar
       y(i)=ystart(i)
      enddo
c .... we start here...
10    continue
       call bs_derivs(nvar,x,y,dydx)
       do i=1,nvar
        yscal(i)=dabs(y(i))+dabs(h*dydx(i))+TINY
       enddo
       if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x
c .... here the BS step with size h is called ........................
       call bsstep(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,ifailfl)
       if (ifailfl.ne.0) return
c .... here the BS step with size h is called ........................
       if((x-x2)*(x2-x1).ge.0.d0)then
        do i=1,nvar
         ystart(i)=y(i)
        enddo
        return
       endif
       if(dabs(hnext).lt.hmin) then
c        write(*,1111)
c        write(*,*)hdid,hnext,hmin
        ifailfl=1
        return
       endif
       h=hnext
      goto 10
1111  format('stepsize smaller than minimum in odeint')
      end
c ====================================================================
      SUBROUTINE bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,ifailfl)
      implicit none
      INTEGER nv,NMAX,KMAXX,IMAX,ifailfl
      REAL*8 eps,hdid,hnext,htry,x,SAFE1,SAFE2,
     *REDMAX,REDMIN,TINY,SCALMX
      PARAMETER (NMAX=7,KMAXX=8,IMAX=KMAXX+1,
     *SAFE1=0.25d0,SAFE2=0.7d0,REDMAX=1.d-5,REDMIN=0.7d0,TINY=1.d-30,
     *SCALMX=0.1d0)
CU    USES mmid,pzextr
      INTEGER i,iq,k,kk,km,kmax,kopt,nseq(IMAX)
      REAL*8 eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew,
     *a(IMAX),alf(KMAXX,KMAXX),err(KMAXX),yerr(NMAX),ysav(NMAX),
     *yseq(NMAX),dydx(NMAX),y(NMAX),yscal(NMAX)
      LOGICAL first,reduct
      SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew
      DATA first/.true./,epsold/-1.d0/
      DATA nseq /2,4,6,8,10,12,14,16,18/
c ---------------------------------------------------------------
      ifailfl=0 ! check on fail

      if (eps.ne.epsold) then
c ... this is basically done only once and then saved
       hnext=-1.d29
       xnew=-1.d29
       eps1=SAFE1*eps
       a(1)=nseq(1)+1
       do k=1,KMAXX
        a(k+1)=a(k)+nseq(k+1)
       enddo
       do iq=2,KMAXX
        do k=1,iq-1
         alf(k,iq)=eps1**((a(k+1)-a(iq+1))/((a(iq+1)-a(1)+1.)*(2*k+1)))
        enddo
       enddo
       epsold=eps
       do kopt=2,KMAXX-1
        if(a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt)) goto 1
       enddo
1      kmax=kopt
      endif
c ...
      h=htry
      do i=1,nv
       ysav(i)=y(i)
      enddo
      if (h.ne.hnext.or.x.ne.xnew) then
       first=.true.
       kopt=kmax
      endif
      reduct=.false.
2     do 17 k=1,kmax
       xnew=x+h
       if (xnew.eq.x) then
ccc         write(*,1111)
         ifailfl=1
         return
       endif    
       call mmid(ysav,dydx,nv,x,h,nseq(k),yseq)  ! -> mmid()
       xest=(h/nseq(k))**2
       call pzextr(k,xest,yseq,y,yerr,nv)        ! -> pzextr()
       if(k.ne.1)then
        errmax=TINY
        do i=1,nv
         errmax=max(errmax,abs(yerr(i)/yscal(i)))
        enddo
        errmax=errmax/eps
        km=k-1
        err(km)=(errmax/SAFE1)**(1./(2*km+1))
       endif
       if(k.ne.1.and.(k.ge.kopt-1.or.first))then
        if(errmax.lt.1.)goto 4
        if(k.eq.kmax.or.k.eq.kopt+1)then
         red=SAFE2/err(km)
         goto 3
        else if(k.eq.kopt)then
         if(alf(kopt-1,kopt).lt.err(km))then
          red=1.d0/err(km)
          goto 3
         endif
        else if(kopt.eq.kmax)then
         if(alf(km,kmax-1).lt.err(km))then
          red=alf(km,kmax-1)*SAFE2/err(km)
          goto 3
         endif
        else if(alf(km,kopt).lt.err(km))then
         red=alf(km,kopt-1)/err(km)
         goto 3
        endif
       endif
17    continue
3     red=min(red,REDMIN)
      red=max(red,REDMAX)
      h=h*red
      reduct=.true.
      goto 2
4     x=xnew
      hdid=h
      first=.false.
      wrkmin=1.d35
      do kk=1,km
       fact=max(err(kk),SCALMX)
       work=fact*a(kk+1)
       if(work.lt.wrkmin)then
        scale=fact
        wrkmin=work
        kopt=kk+1
       endif
      enddo
      hnext=h/scale
      if(kopt.ge.k.and.kopt.ne.kmax.and..not.reduct)then
       fact=max(scale/alf(kopt-1,kopt),SCALMX)
       if(a(kopt+1)*fact.le.wrkmin)then
        hnext=h/fact
        kopt=kopt+1
       endif
      endif

1111  format('step size underflow in bsstep')
      
      return
      END
c ====================================================================
      SUBROUTINE mmid(y,dydx,nvar,xs,htot,nstep,yout)
      implicit none
      INTEGER nstep,nvar,NMAX
      REAL*8 htot,xs
      PARAMETER (NMAX=7)
      INTEGER i,n
      REAL*8 h,h2,swap,x
      REAL*8 ym(NMAX),yn(NMAX),dydx(NMAX),y(NMAX),yout(NMAX)
c ---------------------------------------------------------------
      h=htot/nstep
      do i=1,nvar
       ym(i)=y(i)
       yn(i)=y(i)+h*dydx(i)
      enddo
      x=xs+h
      call bs_derivs(nvar,x,yn,yout)
      h2=2.d0*h
      do n=2,nstep
       do i=1,nvar
        swap=ym(i)+h2*yout(i)
        ym(i)=yn(i)
        yn(i)=swap
       enddo
       x=x+h
       call bs_derivs(nvar,x,yn,yout)
      enddo
      do i=1,nvar
       yout(i)=0.5d0*(ym(i)+yn(i)+h*yout(i))
      enddo
      return
      END
c ====================================================================
      SUBROUTINE pzextr(iest,xest,yest,yz,dy,nv)
      implicit none
      INTEGER iest,nv,IMAX,NMAX
      REAL*8 xest
      PARAMETER (IMAX=13,NMAX=7)
      INTEGER j,k1
      REAL*8 delta,f1,f2,q,d(NMAX),qcol(NMAX,IMAX),x(IMAX)
      REAL*8 dy(NMAX),yest(NMAX),yz(NMAX)
      SAVE qcol,x
c ---------------------------------------------------------------
      x(iest)=xest
      do j=1,nv
       dy(j)=yest(j)
       yz(j)=yest(j)
      enddo
      if(iest.eq.1) then
       do j=1,nv
        qcol(j,1)=yest(j)
       enddo
      else
       do j=1,nv
        d(j)=yest(j)
       enddo
       do k1=1,iest-1
        delta=1.d0/(x(iest-k1)-xest)
        f1=xest*delta
        f2=x(iest-k1)*delta
        do j=1,nv
         q=qcol(j,k1)
         qcol(j,k1)=dy(j)
         delta=d(j)-q
         dy(j)=f1*delta
         d(j)=f2*delta
         yz(j)=yz(j)+dy(j)
        enddo
       enddo
       do j=1,nv
        qcol(j,iest)=dy(j)
       enddo
      endif
      return
      END

c ===============================================================
      function ran3(idum)
      implicit real*8(a-h,o-z)
      parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.d-9)
      dimension ma(55)
      data iff /0/
      save
c --------------------------------------------------------------------
      if((idum.ge.0).and.(iff.ne.0)) goto 20
        iff=1
        mj=mseed-iabs(idum)
        mj=mod(mj,mbig)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.mz)mk=mk+mbig
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
20    continue
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.mz)mj=mj+mbig
      ma(inext)=mj
      ran3=mj*fac
      return
      end
