module driver

  use typedef, only: &
       int64, real64

  use errors, only: &
       IERR_OK, &
       IERR_NONCONV, &
       IERR_TIMEOUT, &
       IERR_YSCALE, &
       IERR_SMALLSTEP, &
       IERR_TERMINATE

  use parameters, only: &
       verb_driver, &
       verb_eta

  implicit none

  ! iextdefault = 0 - ExtPol
  ! iextdefault = 1 - ExtRat

  integer(kind=int64), save :: &
       iextdefault = 0_int64, &
       iflagsdefault = 0_int64

  integer(kind=int64), parameter :: &
       IFLAG_TRUNCATE = 0_int64

  integer(kind=int64), parameter :: &
       tmaxdefault  = 9223372036854775807_int64, &
       nstepdefault = 9223372036854775807_int64, &
       nsavedefault = 9223372036854775807_int64

  real(kind=real64), parameter :: &
       dtmindefault =  1.d0, &
       dtmaxdefault =  1.d99, &
       fbackupcut   = 0.1d0

  real(kind=real64), parameter :: &
       GOLDEN  = 0.5d0 * (sqrt(5.d0) + 1.d0), &
       GOLDENI = 0.5d0 * (sqrt(5.d0) - 1.d0), &
       DTCUT = 0.125d0

contains

  subroutine step_driver(f, r, t0, y0, dt, dt0_, dtmin_, dtmax_, yscale, &
       eps, nstep, nsave, tmax_, n, yy, tt, hnew, neval, nloop, mstep, ierr, iext)

    ! nstep  > 0 - output those exact points dt apart, n+1 including (t0,y0)
    ! nstep == 0 - return 1 step of dt0, (testing)
    ! integrate from t0 to t0 + dt * nstep

    use typedef, only: &
         real64, int64

    use integrate, only: &
         template, regularize, bs

    use formatting, only: &
         fl => ihformat, &
         ft => itformat

    use integrator_flags, only: &
         terminate_integration

    use states, only: &
         status, &
         STATUS_OK, STATUS_COLLIDE, STATUS_INTERACT

    implicit none

    character(len=11), parameter :: &
         name = 'step_driver'

    real(kind=real64), parameter :: &
         epsh = 1.d-12

    procedure(template) :: &
         f
    procedure(regularize), pointer :: &
         r
    integer(kind=int64), intent(in) :: &
         n, nstep, nsave, tmax_
    real(kind=real64), intent(in) :: &
         t0, dt, dt0_, eps, dtmin_, dtmax_
    real(kind=real64), intent(in), dimension(n) :: &
         y0, yscale
    integer(kind=int64), intent(in), optional :: &
         iext

    integer(kind=int64), intent(out) :: &
         ierr, neval, nloop, mstep
    real(kind=real64), intent(out) :: &
         hnew
    real(kind=real64), intent(out), dimension(nsave) :: &
         tt
    real(kind=real64), intent(out), dimension(n,nsave) :: &
         yy

    ! local valriables

    real(kind=real64), dimension(n) :: &
         y, y1, ryscale
    real(kind=real64) :: &
         t, h, h1, t1, t2, ftr, dt0, dtmin, dtmax
    logical :: &
         ltr, lr
    integer(kind=int64) :: &
         iext1, jeval, jerr, i, j, k
    integer(kind=int64) :: &
         count, count_rate, count_max, time0, time1, time2, eta, tmax

    if (nsave < 1) &
         error stop '[' // name // '] requires nsave > 0'

    if (present(iext)) then
       iext1 = iext
    else
       iext1 = iextdefault
    endif

    if (tmax_ <= 0) then
       tmax = tmaxdefault
    else
       tmax = tmax_
    endif

    neval = 0
    nloop = 0

    if (any(yscale(:)==0.d0)) then
       ierr = IERR_YSCALE
       return
    endif
    ryscale = 1.d0 / yscale

    lr = associated(r)

    t = t0
    y = y0
    if (lr) call r(y)
    t1 = t0 + dt * nstep
    t2 = t0
    mstep = 0

    ! time reversal
    ltr = (dt < 0.d0)
    if (ltr) then
       ftr = -1.d0
       if (verb_driver) then
          write(6, "(' [',A,'] Running backward in time.')") &
               name
       endif
    else
       ftr = 1.d0
    endif

    dt0   = sign(dt0_, ftr)

    if (dt  == 0.d0) error stop '[' // name // '] require non-zero time step dt'
    if (dt0 == 0.d0) error stop '[' // name // '] require non-zero inital time step dt0'

    call system_clock(count, count_rate, count_max)
    time0 = count / count_rate
    time1 = time0

    if (nstep == 0) then
       h = dt0
       call bs(f, t, y, h, eps, ryscale, iext1, n, y1, hnew, jeval, jerr)
       neval = neval + jeval
       nloop = nloop + 1
       yy(:,1) = y1
       tt(1) = t + h
       mstep = 1
       ierr = jerr
       return
    endif

    dtmin = dtmin_
    if (dtmin <= 0.d0) then
       if (dtmin == 0.d0) then
          dtmin = dt0
       else
          dtmin = dtmindefault
       endif
       if (verb_driver) then
          write(6, "(' [',A,'] Setting minimum time step to ',1p,e13.5)") &
               name, dtmin
       endif
    endif

    dtmax = dtmax_
    if (dtmax <= 0.d0) then
       if (dtmax == 0.d0) then
          dtmax = dt0
       else
          dtmax = dtmaxdefault
       endif
       if (verb_driver) then
          write(6, "(' [',A,'] Setting maximum time step to ',1p,e13.5)") &
               name, dtmax
       endif
    endif

    if (dtmin > dtmax) error stop '[' // name // '] dtmin > dtmax'

    dtmin = sign(dtmin, ftr)
    dtmax = sign(dtmax, ftr)

    terminate_integration = .false.
    tt(1) = t
    yy(:,1) = y(:)
    hnew = dt0
    h1 = 0.d0
    mstep = 1

    do i = 2, nstep + 1
       h1 = h1 + dt
       j = 0
       do while ((h1 > epsh * dt).neqv.ltr)
          j = j + 1
          if ((hnew * GOLDEN > h1).neqv.ltr) then
             h = h1
          else
             if ((hnew > h1).neqv.ltr) then
                h = h1
             else
                h = hnew
             endif
          endif
          call bs(f, t, y, h, eps, ryscale, iext1, n, y1, hnew, jeval, jerr)
          neval = neval + jeval
          if (terminate_integration) then
             if (((status == STATUS_COLLIDE).or.(status == STATUS_INTERACT)).and. &
                  ((h > dtmin * (1.d0 + epsh)).neqv.ltr)) then
                status = STATUS_OK
                terminate_integration = .False.
                hnew = max(h * DTCUT, dtmin)
                if (verb_driver) then
                   write(6, "(' [',A,'] reducing timestep due to collision to',1p,e13.5)") &
                        name, hnew
                endif
                cycle
             endif
             if (verb_driver) then
                write(6, "(' [',A,'] terminating.')") &
                     name
             endif
             ierr = IERR_TERMINATE
             goto 1000
          endif
          call system_clock(count, count_rate, count_max)
          time2 = count / count_rate
          if (time2 - time1 >= 10) then
             eta = nint((t1 - t) * (time2 - time1) / (t - t2))
             t2 = t
             if (verb_eta) then
                write(6, "(' [',A,'] ',a7,a18,1p,e23.15,e13.5,i5,a20,a8)")  &
                     name, ft(time2-time0,7), fl(i,18), t, h, jeval, fl(neval,20), ft(eta,8)
             endif
             time1 = time2
             if (time2 - time0 >= tmax) then
                if (verb_driver) then
                   write(6, "(' [',A,'] maximum execution time exceeded: ',a7)") &
                        name, ft(tmax,7)
                endif
                ierr = IERR_TIMEOUT
                goto 1000
             endif
          end if
          if ((hnew < dtmin).neqv.ltr) then
             if (verb_driver) then
                write(6, "(' [',A,'] timestep too small',1p,e13.5)") &
                     name, hnew
             endif
             ierr = IERR_SMALLSTEP
             goto 1000
          endif
          if (jerr == 0) then
             t = t + h
             if (lr) call r(y1)
             y = y1
             h1 = h1 - h
             if ((h == hnew).and.((hnew < dtmax).neqv.ltr)) then
                hnew = hnew * GOLDEN
             end if
             if ((hnew > dtmax).neqv.ltr) then
                hnew = dtmax
             endif
             cycle
          end if
          if (jerr == -1) then
             hnew = h * DTCUT
             cycle
          endif
          if (jerr == -2) then
             hnew = hnew * fbackupcut
             cycle
          endif
          if (verb_driver) then
             write(6, "(' [',A,'] integrator error ',I3)") &
                  name, jerr
          endif
          ierr = jerr
          goto 1000
       end do
       k = mod(i - 1, nsave) + 1
       mstep = i
       tt(k) = t
       yy(:,k) = y
       nloop = nloop + j
    enddo
    ierr = IERR_OK
    goto 2000

1000 continue
    nloop = nloop + j

2000 continue
    if (mstep > nsave) then
       k = mod(mstep, nsave)
       tt(:)   = cshift(tt,shift=k)
       yy(:,:) = cshift(yy,shift=k,dim=2)
    endif
    return

  end subroutine step_driver


  subroutine dynamic_driver(f, r, t0, y0, dt, dt0_, dtmin_, dtmax_, dtd_, yscale, &
       eps, nstep_, nsave_, tmax_, n, hnew, neval, nloop, ierr, iext, iflags)

    ! integrate from t0 to t0 + dt, max nstep steps
    ! store all data in module

    use typedef, only: &
         real64, int64

    use integrate, only: &
         template, regularize, bs

    use formatting, only: &
         fl => ihformat, &
         ft => itformat

    use integral_data, only: &
         initialize, append

    use integrator_flags, only: &
         terminate_integration

    use states, only: &
         status, &
         STATUS_OK, STATUS_COLLIDE, STATUS_INTERACT

    implicit none

    character(len=14), parameter :: &
         name = 'dynamic_driver'

    real(kind=real64), parameter :: &
         epsh = 1.d-13

    procedure(template) :: &
         f
    procedure(regularize), pointer :: &
         r
    integer(kind=int64), intent(in) :: &
         n, nstep_, nsave_, tmax_
    real(kind=real64), intent(in) :: &
         t0, dt, dt0_, eps, dtd_, dtmax_, dtmin_
    real(kind=real64), intent(in), dimension(n) :: &
         y0, yscale
    integer(kind=int64), intent(in), optional :: &
         iext, iflags

    real(kind=real64), intent(out) :: &
         hnew
    integer(kind=int64), intent(out) :: &
         ierr, neval, nloop

    ! local valriables

    real(kind=real64), dimension(n) :: &
         y, y1, ryscale
    real(kind=real64) :: &
         t, h, t1, t1m, t1p, dtx, t2, ftr, dtmin, dtmax, dt0, dtd
    logical :: &
         ltr, ltruncate, lr
    integer(kind=int64) :: &
         iext1, jeval, jerr, tmax
    integer(kind=int64) :: &
         i, nstep, nsave, iflags1
    integer(kind=int64) :: &
         count, count_rate, count_max, time0, time1, time2, eta

    neval = 0
    nloop = 0

    if (any(yscale(:)==0.d0)) then
       ierr = IERR_YSCALE
       return
    endif
    ryscale = 1.d0 / yscale

    if (present(iext)) then
       iext1 = iext
    else
       iext1 = iextdefault
    endif

    if (present(iflags)) then
       iflags1 = iflags
    else
       iflags1 = iflagsdefault
    endif
    ltruncate = .not.btest(iflags1, IFLAG_TRUNCATE)
    lr = associated(r)

    if (tmax_ <= 0) then
       tmax = tmaxdefault
    else
       tmax = tmax_
    endif

    if (nsave_ <= 0) then
       nsave = nsavedefault
    else
       nsave = nsave_
    endif

    call initialize(n, nsave)

    if (nstep_ <= 0) then
       nstep = nstepdefault
    else
       nstep = nstep_
    endif

    t1 = t0 + dt
    t1m = t1 * (1.d0 - epsh)
    t1p = t1 * (1.d0 + epsh)
    t2 = t0

    ! time reversal
    ltr = (dt < 0.d0)
    if (ltr) then
       ftr = -1.d0
       if (verb_driver) then
          write(6, "(' [',A,'] Running backward in time.')") &
               name
       endif
    else
       ftr = 1.d0
    endif

    dtmin = dtmin_
    if (dtmin <= 0.d0) then
       if (dtmin == 0.d0) then
          dtmin = dt0
       else
          dtmin = dtmindefault
       endif
       if (verb_driver) then
          write(6, "(' [',A,'] Setting minimum time step to ',1p,e13.5)") &
               name, dtmin
       endif
    endif

    dtmax = dtmax_
    if (dtmax <= 0.d0) then
       if (dtmax == 0.d0) then
          dtmax = dt0
       else
          dtmax = dtmaxdefault
       endif
       if (verb_driver) then
          write(6, "(' [',A,'] Setting maximum time step to ',1p,e13.5)") &
               name, dtmax
       endif
    endif

    if (dtmin > dtmax) error stop '[' // name // '] dtmin > dtmax'

    dt0   = sign(dt0_ , ftr)
    dtd   = sign(dtd_ , ftr)
    dtmin = sign(dtmin, ftr)
    dtmax = sign(dtmax, ftr)

    if (dt  == 0.d0) error stop '[' // name // '] require non-zero time step dt'
    if (dt0 == 0.d0) error stop '[' // name // '] require non-zero inital time step dt0'

    call system_clock(count, count_rate, count_max)
    time0 = count / count_rate
    time1 = time0

    terminate_integration = .false.
    ierr = IERR_OK
    dtx = 0.d0
    t = t0
    y = y0
    if (lr) call r(y)
    call append(t, y)
    hnew = dt0
    do i = 2, nstep + 1
       h = hnew
       if (ltruncate.and.((t + h > t1).neqv.ltr)) then
          h = t1p - t
       endif
       call bs(f, t, y, h, eps, ryscale, iext1, n, y1, hnew, jeval, jerr)
       neval = neval + jeval
       if (terminate_integration) then
          if (((status == STATUS_COLLIDE).or.(status == STATUS_INTERACT)).and. &
               ((h > dtmin * (1.d0 + epsh)).neqv.ltr)) then
             status = STATUS_OK
             terminate_integration = .False.
             hnew = max(h * DTCUT, dtmin)
             if (verb_driver) then
                write(6, "(' [',A,'] reducing timestep due to collision to',1p,e13.5)") &
                     name, hnew
             endif
             cycle
          endif
          ierr = IERR_TERMINATE
          if (verb_driver) then
             write(6, "(' [',A,'] terminating.')") &
                  name
          endif
          exit
       endif
       call system_clock(count, count_rate, count_max)
       time2 = count / count_rate
       if (time2 - time1 >= 10) then
          eta = nint((t1 - t) * (time2 - time1) / (t - t2))
          t2 = t
          if (verb_eta) then
             write(6, "(' [',A,'] ',a7,a18,1p,e23.15,e13.5,i5,a20,a8)") &
                  name, ft(time2-time0,7), fl(i,18), t, h, jeval, fl(neval,20), ft(eta,8)
          endif
          time1 = time2
          if (time2 - time0 >= tmax) then
             if (verb_driver) then
                write(6, "(' [',A,'] maximum execution time exceeded: ',a7)") &
                     name, ft(tmax,7)
             endif
             ierr = IERR_TIMEOUT
             exit
          endif
       end if
       if ((hnew < dtmin).neqv.ltr) then
          ierr = IERR_SMALLSTEP
          if (verb_driver) then
             write(6, "(' [',A,'] timestep too small',1p,e13.5)") &
                  name, hnew
          endif
          exit
       endif
       if (jerr == 0) then
          t = t + h
          if (lr) call r(y1)
          y = y1
          dtx = dtx + h
          if ((dtx > dtd).neqv.ltr) then
             call append(t, y)
             dtx = 0.d0
          endif
          if ((t > t1m).neqv.ltr) exit
          if ((h == hnew).and.((hnew < dtmax).neqv.ltr)) then
             hnew = hnew * GOLDEN
          end if
          if ((hnew > dtmax).neqv.ltr) then
             hnew = dtmax
          endif
          cycle
       end if
       if (jerr == -1) then
          hnew = hnew * DTCUT
          cycle
       endif
       if (jerr == -2) then
          hnew = hnew * fbackupcut
          cycle
       endif
       ierr = jerr
       exit
    end do
    if ((dtx > 0.d0).neqv.ltr) call append(t, y)
    nloop = i - 1
    return
  end subroutine dynamic_driver

end module driver