!-------------------------------------------------------------------- 
!                    Run  N-body simulations
!                    using integration scheme with variable time-step
!
!                    Logic: - read parameters, coordinates, velocities, 
!                             masses from file 'particles.dat' 
!                           - run for RunTime
!                           - save the data for continuation of the run 
!-------------------------------------------------------------------- 
PROGRAM Move2
Use DataStructures
REAL*8    Ekin,Epot,Etot,center(3)

                          ! open files
      Open(1,file='DATA/particles.dat',form='unformatted')
      Open(2,file='DATA/trajectories.dat')
      Open(3,file='DATA/energies.dat',position='append')
      Open(8,file='DATA/analysis.dat',position='append')

      CALL ReadMoment('DATA/particlesA.dat')

      write(*,*) ' Number of particles=',Np,' Time=',time
      write(*,*) ' Time step          =',dtstep,' Step=',iStep
      write(*,*) ' Force softening    =',eps
      write(*,*) ' Run Time           =',RunTime

                          ! set parameters
      Nsteps     = INT(RunTime/dtstep)  ! number of steps for this run
      PrintTime  = time +dPrintTime         ! next output time
      Ntraj      = min(5,Np)                          ! save this number of trajectories
 
      indx       = 0 
      E0          =-0.4802930                       ! initial energy
      Do i=1,Nsteps                                  ! main loop
         Call SelectStep                            !   assign value of time-step to each particle
         Do
            Call MoveParticles 
             If(iMultiple>0)Then
                Call TestRestart(iRestart)  ! if particles move too fast, restart current step
                If(iRestart == 0)exit
             EndIf
         EndDo
         time = time + dtstep
         iStep= iStep+ 1
         iOrbit = INT(time/Period)

         If(time.ge.PrintTime)Then                                     ! make analysis when time is right
            Call Energies(Ekin,Epot) ; Etot = Ekin+Epot 
            write(*,20) &
            time,iStep,Ekin,Epot,(Etot/E0-1)*100.,Ekin/Etot,MaxLev,Naccel/float(Np),(iLevel(j),j=MaxLev,0,-1)
            write(3,20) &
            time,iStep,Ekin,Epot,Etot,Ekin/Etot,MaxLev,Naccel/float(Np),(iLevel(j),j=MaxLev,0,-1)
            Call Analyze(indx,center)
            Call SaveMoment('DATA/particlesA.dat')
            Call WriteSnapshot(center)
            Call WriteTrajectories(Ntraj,center)     ! this dumps very 100th trajectory
            PrintTime  = PrintTime + dPrintTime
        EndIf 
      EndDo
stop 
      Call SaveMoment('DATA/particlesA.dat')
 20   format(f11.3,i11,4g15.7,i4,f10.1,2i6,20i4)

      Stop
      End Program MOVE2
!-------------------------------------------------------------------- 
!               Find Acceleations for particles on given Level
!               Coordinates of Np particles are xa,ya,za
!                They must be at the same moment of time
!               Accelerations are placed to {gx,gy,gz}
      SUBROUTINE GetAccelerations(Level)
!-------------------------------------------------------------------- 
 use DataStructures
      REAL*8   rr,ax,ay,az,x,y,z
!Naccel =Naccel +iLevel(Level)  ! count how many time we evaluate accelerations
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,rr,ax,ay,az,x,y,z)
       Do i=1,Np       ! sum contributions for particles on given level
         If(LevPart(i).eq.Level)Then
            ax = 0.; ay= 0.; az = 0.
            x =xa(i); y=ya(i); z=za(i)
            Do j=1,Np
               rr = aMass(j)/(sqrt( (x-xa(j))**2+(y-ya(j))**2+ (z-za(j))**2+ eps2 ) )**3
               ax =ax  -rr*(x-xa(j))
               ay =ay  -rr*(y-ya(j))
               az =az  -rr*(z-za(j))
            EndDo
            gx(i) =ax; gy(i) = ay; gz(i) = az 

          EndIf
      EndDo 
      Return
      End SUBROUTINE GetAccelerations
!-------------------------------------------------------------------- 
!               Find Acceleations for all particles 
!               If  indx == 0 => {xc,yc,zc} -> {gx,gy,gz}
!               otherwise: {xa,ya,za} -> {gxa,gya,gza}
      SUBROUTINE GetAccel(indx)
!-------------------------------------------------------------------- 
 use DataStructures
      REAL*8   rr,ax,ay,az,x,y,z

!Naccel =Naccel +Np  ! count how many time we evaluate accelerations
If(indx == 0)Then     !      Use the main set of coords and accelerations
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,rr,ax,ay,az,x,y,z) 
       Do i=1,Np       ! sum contributions for particles on given level
            ax = 0.; ay= 0.; az = 0.
            x =xc(i); y=yc(i); z=zc(i)
            Do j=1,Np
               rr = aMass(j)/(sqrt( (x-xc(j))**2+(y-yc(j))**2+ (z-zc(j))**2+ eps2 ) )**3
               ax =ax  -rr*(x-xc(j))
               ay =ay  -rr*(y-yc(j))
               az =az  -rr*(z-zc(j))
            EndDo
            gx(i) =ax; gy(i) = ay; gz(i) = az 
      EndDo 
Else    !      Use the auxiliary set of coords and accelerations
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,rr,ax,ay,az,x,y,z) 
       Do i=1,Np       ! sum contributions for particles on given level
            ax = 0.; ay= 0.; az = 0.
            x =xa(i); y=ya(i); z=za(i)
            Do j=1,Np
               rr = aMass(j)/(sqrt( (x-xa(j))**2+(y-ya(j))**2+ (z-za(j))**2+ eps2 ) )**3
               ax =ax  -rr*(x-xa(j))
               ay =ay  -rr*(y-ya(j))
               az =az  -rr*(z-za(j))
            EndDo
            gxa(i) =ax; gya(i) = ay; gza(i) = az 
      EndDo 
EndIf
      End SUBROUTINE GetAccel
!-------------------------------------------------------------------- 
!          Move particles on all levels for one large 
!          time-step dt. Number of current levels is
!          0 - MaxLev.
!          Zero level is has the smallest time-step dmin
!          MaxLev has time-step dt     
      SUBROUTINE MoveParticles
!-------------------------------------------------------------------- 
use DataStructures
Real*8 :: dmin,dtime,t0,t,dd

If(iMultiple == 0 .or. MaxLev == 0)Then          ! No mutistepping
dt0 = dtstep/2.
!$OMP PARALLEL DO DEFAULT(SHARED)                                    
          Do i=1,Np
               xc(i) =xa(i)+vxa(i)*dt0
               yc(i) =ya(i)+vya(i)*dt0
               zc(i) =za(i)+vza(i)*dt0
               vx(i) =vxa(i);  vy(i) =vya(i); vz(i)   =vza(i)
          EndDo
          Return
EndIf

                                   ! Multistepping -----------------------------
!$OMP PARALLEL DO DEFAULT(SHARED)  
Do i=1,Np                           ! initialize arrays
   timeP(i) =0.
   iLevPart(i) =0.
EndDo

Nsteps =2**MaxLev                  ! number of small steps to make
dmin  = dtstep/2.**MaxLev          ! smallest time-step
t0 = 0.
   do i=0,(Nsteps-1)                        ! loop over all possible small steps
      t = t0+i*dmin                               ! time relative to the begining of step
      Lev1 =0
      Do Level=MaxLev,1,-1            ! find which level to advance
            dd =i/float(2**Level)-0.5
            if(abs(dd-INT(dd))<1.d-5)Then
                Lev1 =Level
                dtime = dmin*2**Level
                exit
            EndIf
       EndDo
       If(Lev1.ne.0)Call MoveLevel(t-dtime/2.,dtime,Level)! if found, move Lev1
                              Call MoveLevel(t,dmin,0)                         ! always move zero level
   enddo
return
end SUBROUTINE MoveParticles
!----------------------------------------------------------------------------------------- 
!           Change time-steps and levels of particles.
!            MaxLev = current maximum level of stepping
!               level= 0            <-- smallest step
!               level=MaxLev <-- largest   step = dtstep
!            
      SUBROUTINE SelectStep
!---------------------------------------------------------------------------------------- 
 use DataStructures
                      ! construct coordinates {xa,ya,za} at t1 =t0 +dt/2
dt0 = dtstep/2.
d2  = dtstep**2/2.

!$OMP PARALLEL DO DEFAULT(SHARED)                                      
          Do i=1,Np
                   xa(i) =xc(i)+vx(i)*dt0
                   ya(i) =yc(i)+vy(i)*dt0
                   za(i) =zc(i)+vz(i)*dt0
          EndDo
   Call GetAccel(1)         ! get accelerations at t0+dt/2
                                          ! find velocities at moment t1 
!$OMP PARALLEL DO DEFAULT(SHARED)                                     
          Do i=1,Np
                   vxa(i) = vx(i) +gxa(i)*dtstep        ! kick
                   vya(i) = vy(i) +gya(i)*dtstep
                   vza(i) = vz(i) +gza(i)*dtstep
          EndDo
If(iMultiple == 0) Return                                 ! exit if single step is used
                             !-------------------------------------------------------------------------------------
                                                          ! find levels of particles
   dmin    =  dtstep/2.**MaxLev   ! current smalest time-step     
   dvmin   = 1.d+10
   dvmax   = 0.
   LevMin  =  100
   LevMax  = -100                   
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(i,v,va,g,dt,dx,dv,k) &
!$OMP REDUCTION(MIN:LevMin,dvmin) &                     
!$OMP REDUCTION(MAX:LevMax,dvmax)                      
          Do i=1,Np
                v = sqrt(vx(i)**2 +vy(i)**2 +vz(i)**2)    ! velocity   0
                g = sqrt((gxa(i))**2 +(gya(i))**2 +(gza(i))**2)    ! velocity   0
               va= sqrt(vxa(i)**2+vya(i)**2+vza(i)**2)  ! velocity   1
                dt = dmin*2**(LevPart(i))
                v= max(v,va)
               dx =v*dt
               dv =g*dt

              If(dx>1.5*aStep.or.dv>1.5*Alph*v)Then
                          do k=1,5                                            ! k is the number of levels to add for refinement
                               if(dv< Alph*(k+1.5)*v.and.dx<aStep*(k+1.5))exit
                          enddo
                   LevPart(i) =Max(LevPart(i)-k,MaxLev-iMultiple)  ! make time-step smaller

               Else if(dx<0.75*aStep.and.dv<0.75*Alph*v)Then
                   LevPart(i) =min(LevPart(i)+1,MaxLev)                  ! increase time-step
               EndIf

               LevMin =min(LevPart(i),LevMin)
               LevMax =max(LevPart(i),LevMax)
          EndDo
MaxLev = MaxLev-LevMin                         ! update current max level
!write(*,*) ' MaxLev now=',MaxLev,' min=',LevMin
!write(*,*) '                 dv max/min=',dvmax,dvmin

!$OMP PARALLEL DO DEFAULT(SHARED)                                      
          Do i=1,Np
             LevPart(i) =LevPart(i) -LevMin
          EndDo
          iLevel    = 0                                ! count particles on different levels
          Do i=1,Np
             iLevel(LevPart(i)) =iLevel(LevPart(i)) +1
          EndDo

!$OMP PARALLEL DO DEFAULT(SHARED)                                      
           Do  i=1,Np                               ! save {x} and {V} for re-starting the time-step
                   vxa(i) = vx(i)  
                   vya(i) = vy(i) 
                   vza(i) = vz(i) 
                   gxa(i) = xc(i)  
                   gya(i) = yc(i) 
                   gza(i) = zc(i)
            EndDo
 
end SUBROUTINE SelectStep

!----------------------------------------------------------------------------------------- 
!                     Decide whether restarting is needed
!           
!                     If re-starting: restore coords and velocities
!                                             and update refinements levels
!            
      SUBROUTINE TestRestart(iRestart)
!---------------------------------------------------------------------------------------- 
 use DataStructures
                      ! construct coordinates {xa,ya,za} at t1 =t0 +dt/2
iRestart =0
i1 = 0; i2 =0; i3 =0; i4 =0
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(i) &
!$OMP  REDUCTION(+:i1,i2,i3,i4)                   
           Do  i=1,Np              
               iLevPart(i) = MIN(iLevPart(i),iMultiple+LevPart(i)-MaxLev)
               If(iLevPart(i)==1)i1=i1+1
               If(iLevPart(i)==2)i2=i2+1
               If(iLevPart(i)==3)i3=i3+1
               If(iLevPart(i)==4)i4=i4+1
            EndDo

 If(i3+i4==0)return                            ! exit if re-start not needed
                             !--------------------------------------------------------------------
          write(*,'(20x,a,6i5)') 'Restarting: number of particles adjustments=',i1,i2,i3,i4
                                                          ! find levels of particles
iRestart = i3+i4
LevMin = 0
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(i) &
!$OMP REDUCTION(MIN:LevMin)                      
          Do i=1,Np
                   LevPart(i) =Max(LevPart(i)-iLevPart(i),MaxLev-iMultiple)  ! make time-step smaller
                   LevMin =min(LevPart(i),LevMin)
          EndDo
MaxLev = MaxLev-LevMin                         ! update current max level

!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(i) 
          Do i=1,Np
             LevPart(i) =LevPart(i) -LevMin
          EndDo
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(i) SCHEDULE(DYNAMIC,100)                     
           Do  i=1,Np                               ! restore {x} and {V} for re-starting the time-step
                   vx(i) = vxa(i)  
                   vy(i) = vya(i) 
                   vz(i) = vza(i) 
                   xc(i)   = gxa(i)  
                   yc(i)   = gya(i) 
                   zc(i)   = gza(i)
            EndDo
 
end SUBROUTINE TestRestart

!----------------------------------------------------------------------------------------- 
!          Move particles on current level  'Level'.
!          time-step is 'dtime'. 
!          'time' is the beginning of time-step for this level 
!             1) advance all particles to '1/2' moment of level 'Level'
!             2) find accelerations for particles on Level
!             3) advance particles on Level
      SUBROUTINE MoveLevel(t,dtime,Level)
!---------------------------------------------------------------------------------------- 
 use DataStructures
Real*8 :: dmin,dtime,t,time0,thalf,tend,dt0

thalf = t +dtime/2.
tend  = t +dtime
                      ! construct coordinates {xa,ya,za} for given moment time
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(dt0)                      
          Do i=1,Np
               If(LevPart(i).lt.Level)Then ! copy  xc -> xa
                   xa(i) =xc(i); ya(i) = yc(i); za(i) =zc(i)
               Else                                       ! advance coordinates
                   dt0 = thalf-timeP(i)
                   xa(i) =xc(i)+vx(i)*dt0       ! drift
                   ya(i) =yc(i)+vy(i)*dt0
                   za(i) =zc(i)+vz(i)*dt0
               endif
          EndDo
   Call GetAccelerations(Level)
dt =dtime/2.
aa = (dtime/Alph)**2    ! 

!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(i,v0x,v0y,v0z,g2,v2,iAdd)                      
          Do i=1,Np
               If(LevPart(i).eq.Level)Then 
                                                         ! advance coordinates
                   v0x = vx(i);  v0y = vy(i);  v0z = vz(i)
                   vx(i) = v0x  +gx(i)*dtime        ! kick
                   vy(i) = v0y +gy(i)*dtime
                   vz(i) = v0z +gz(i)*dtime

                   xc(i) =xc(i)+(vx(i)+v0x)*dt   ! drift
                   yc(i) =yc(i)+(vy(i)+v0y)*dt
                   zc(i) =zc(i)+(vz(i)+v0z)*dt
                   g2 = gx(i)**2+gy(i)**2 +gz(i)**2   ! acceleration
                   v2 = vx(i)**2+vy(i)**2 +vz(i)**2                   
                   timeP(i) = tend                      ! change time
                   iAdd = 0                                  ! test rate of change of velocity
                   If(g2*aa > 2.25*v2)Then
                          iAdd =1
                          If(g2*aa > 6.25*v2)Then
                              iAdd =2
                              If(g2*aa > 12.25*v2)Then
                                 iAdd =3
                                 If(g2*aa > 20.25*v2)Then
                                    iAdd =4
                                 EndIf
                              EndIf
                           EndIf
                    EndIf
                    iLevPart(i) =MAX(iLevPart(i),iAdd)          ! update map of refinement
               endif
          EndDo
end SUBROUTINE MoveLevel
!-------------------------------------------------------------------- 
!               write Ntraj trajectories into file 2
      SUBROUTINE WriteTrajectories(Ntraj,center)
!-------------------------------------------------------------------- 
 use DataStructures
Real*8 center(3)

      write(2,10)time,iStep,(xc(i)-center(1),yc(i)-center(2),zc(i)-center(3),i=1,Ntraj)

 10   format(g12.5,i8,20(2x,3f9.4))
      Return
End SUBROUTINE WriteTrajectories