module polyadic

  use parameters, only: &
       verb_interact

  use typedef, only: &
       int32, real64

  implicit none

  save

  integer(kind=int32), parameter :: &
       ndim  = 3, &
       nquat = 4, &
       ncorn = 2

  integer(kind=int32) :: &
       nstar, nstar1, nstar2, &
       nedge, ncon, norb, &
       nvorb, nvstar, nqstar

  integer(kind=int32), dimension(:, :), allocatable :: &
       ies, ise, isc, ipe, &
       iij
  integer(kind=int32), dimension(:, :, :), allocatable :: &
       iijk


contains

  subroutine setup_poly(n)

    implicit none

    integer(kind=int32), intent(IN) :: &
         n

    integer(kind=int32) :: &
         i, j, k, k1, k2

    ! stars belonging to each edge (ncorn=2)
    !    ies = reshape((/1,2, 1,3, 2,3, 1,4, 2,4, 3,4/), (/ncorn, nedge/))
    ! edge                1    2    3    4    5    6

    ! edges belonging to each star (ncon=nstar-1)
    !    ise = reshape((/1,2,4, 1,3,5, 2,3,6, 4,5,6/), (/ncon, nstar/))
    ! star                 1      2      3      4

    ! corners belonging to each star (ncon=nstar-1)
    !    isc = reshape((/1,1,1, 2,1,1, 2,2,1, 2,2,2/), (/ncon, nstar/))
    ! star                 1      2      3      4

    ! edges belonging to each star pair combination (symmetric)
    ! integer(kind=int32), dimension(nstar, nstar), parameter :: &
    !      ipe = reshape((/0, 1, 2, 4, &
    !                      1, 0, 3, 5, &
    !                      2, 3, 0, 6, &
    !                      4, 5, 6, 0 /), (/nstar, nstar/))

    ! indices for td3 combinations
    ! integer(kind=int32), dimension(nstar, nstar1), parameter :: &
    !      iij  = reshape((/2,1,1,1, 3,3,2,2, 4,4,4,3/), (/nstar, nstar1/))
    ! integer(kind=int32), dimension(nstar, nstar1, nstar2), parameter :: &
    !      iijk = reshape((/3,3,2,2, 2,1,1,1, 2,1,1,1, &
    !                       4,4,4,3, 4,4,4,3, 3,3,2,2 /), (/nstar, nstar1, nstar2/))

    nstar = n
    nstar1 = nstar - 1
    nstar2 = nstar - 2
    nedge = (nstar1 * nstar) / 2
    ncon = nstar1
    norb  = nstar - 1

    nvorb  = ndim * norb
    nvstar = ndim * nstar
    nqstar = nquat * nstar


    if (allocated(ies)) deallocate(ies)
    allocate(ies(ncorn, nedge))

    if (allocated(ise)) deallocate(ise)
    allocate(ise(ncon, nstar))

    if (allocated(isc)) deallocate(isc)
    allocate(isc(ncon, nstar))

    if (allocated(ipe)) deallocate(ipe)
    allocate(ipe(nstar, nstar))

    if (allocated(iij)) deallocate(iij)
    allocate(iij(nstar, nstar1))

    if (allocated(iijk)) deallocate(iijk)
    allocate(iijk(nstar, nstar1, nstar2))

    k1 = 0
    k2 = 1
    do i=1, nedge
       k1 = k1 + 1
       if (k1 == k2) then
          k2 = k2 + 1
          k1 = 1
       endif
       ies(1, i) = k1
       ies(2, i) = k2
       ise(k2 - 1, k1) = i
       ise(k1    , k2) = i
       isc(k2 - 1, k1) = 1
       isc(k1    , k1) = 2
       ipe(k1, k2) = i
       ipe(k2, k1) = i
    enddo
    do i=1, nstar
       ipe(i,i) = 0
       do j=1, nstar1
          if (j < i) then
             k1 = j
          else
             k1 = j + 1
          endif
          iij(i,j) = k1
          do k=1, nstar2
             k2 = k
             if (k2 ==  i) k2 = k2 + 1
             if (k2 == k1) then
                k2 = k2 + 1
                if (k2 == i) k2 = k2 + 1
             endif
             iijk(i,j,k) = k2
          enddo
       enddo
    enddo

  end subroutine setup_poly


  ! =======================================================================
  ! Direct code, full
  ! =======================================================================

  function directpf(trel, y) result(yp)

    use, intrinsic :: &
         ieee_arithmetic, only: &
         ieee_value, ieee_quiet_nan

    use typedef, only: &
         int32, real64

    use physconst, only: &
         GRAV, CLIGHT

    use stardata, only: &
         IMASS, IRADI, IAPSI, IINER, ITAUQ, &
         toffset, stars, &
         star_huntpol

    use vector, only: &
         cross_product

    use deriv_data, only: &
         local_store, local_data

    use parameters, only: &
         cutoff, has_cutoff, &
         interact, has_interact, &
         epslim, omegalim

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

    use integrator_flags, only: &
         terminate_integration

    use flags_data, only: &
         use_sd, use_td,use_tf, use_pn, &
         use_td3, &
         use_rxv, use_qp, &
         use_ori, use_orq, use_ors, &
         use_pqm, use_pqq, use_cd

    use orientation, only: &
         rotation_w2A, &
         rotation_f2A, &
         rotation_dodt, &
         rotation_dqdt, &
         rotate_diag

    use tree, only: &
         order, hierarchy

    implicit none

    real(kind=real64), intent(in) :: &
         trel
    real(kind=real64), intent(in), dimension(:) :: &
         y

    real(kind=real64), dimension(size(y, 1)) :: &
         yp

    real(kind=real64), parameter :: &
         CLIGHTm1 = 1.d0 / CLIGHT, &
         CLIGHTm2 = CLIGHTm1**2, &
         CLIGHTm4 = CLIGHTm2**2, &
         CLIGHTm5 = CLIGHTm4 * CLIGHTm1, &
         CLIGHTm6 = CLIGHTm5 * CLIGHTm1, &
         Gcm2 = GRAV * CLIGHTm2, &
         Gcm6 = Gcm2 * CLIGHTm4, &
         ONETHIRD = 1.d0 / 3.d0, &
         G3 = GRAV * 3.d0

    ! local variables

    integer(kind=int32) :: &
         i, j, k, l, ir, ij, ik, j1, k1
    real(kind=real64) :: &
         t
    real(kind=real64), dimension(stars(1)%m, nstar) :: &
         star_v
    real(kind=real64), dimension(ndim, norb) :: &
         r_j, v_j, acc_j
    real(kind=real64), dimension(nquat, nstar) :: &
         orq_v, roq_v
    real(kind=real64), dimension(ndim, ndim, nstar) :: &
         a_v
    real(kind=real64), dimension(ndim, nstar) :: &
         i_v, j_v, w_v, torque_v, ori_v, rot_v, &
         q_v
    real(kind=real64), dimension(nstar) :: &
         m_v, s_v, k_v, t_v, mi_v, s5k_v, &
         wn2_v, m2_v
    logical, dimension(nstar) :: &
         lpqm_v
    real(kind=real64), dimension(ncorn, nedge) :: &
         qrr_c
    real(kind=real64), dimension(ndim, ncorn, nedge) :: &
         qr_c

    ! coordinate setup

    real(kind=real64), dimension(ndim, nedge) :: &
         r_e, v_e
    real(kind=real64) :: &
         m_12, mi_12, m_123, mi_123, m_1234
    real(kind=real64), dimension(2) :: &
         mi_d
    real(kind=real64), dimension(2,2) :: &
         mc_y
    real(kind=real64), dimension(ndim, 2,2) :: &
         rc_y, vc_y

    ! forces

    real(kind=real64), dimension(nstar1, nstar) :: &
         gjk_x, mjk_x
    real(kind=real64), dimension(ndim, nedge) :: &
         f_e
    real(kind=real64), dimension(nedge) :: &
         rn_e, rn2_e, rnm1_e, rnm2_e, rnm3_e, rnm4_e, rnm5_e, rnm7_e, rnm8_e, &
         m_e, mp_e, mm2_e, fgr_e, ggr_e
    real(kind=real64) :: &
         f, g, h, &
         eta, phi, vn2, vr, &
         Gm, eps, an, nm1
    real(kind=real64), dimension(ndim) :: &
         v, w, rik, rij, qr
    real(kind=real64), dimension(ncorn, nedge) :: &
         wdr_c
    real(kind=real64), dimension(ndim, ncorn, nedge) :: &
         wxr_c

    ! for recording

    real(kind=real64), dimension(ndim, ncorn, nedge) :: &
         fsd_c, ftd_c, ftf_c, fpq_c, tpq_c
    real(kind=real64), dimension(ndim, nstar, nstar1, nstar2) :: &
         ftd3_d
    real(kind=real64), dimension(ndim, nedge) :: &
         fpn_e, fng_e


    ! [[[star1, star2], star3], star4] or [[star1, star2], [star3, star4]]

    ! FUTURE - note on notations for hierarchies
    ! [[[,],],] +3,-1,-1,-1
    ! [[,],[,]] +2,-1,+1,-2
    ! not supported (but redundant)
    ! [[,[,]],] +2,+1,-2,-1
    ! [,[,[,]]] +1,+1,+1,-3

    ! VARIABLE NAMES
    ! _v - values for stars as last index (nstar)
    ! _j - values for relative star coordnates in Jacobi coordinates
    !? _x - helper arrays of arbitrary dimensions
    ! _e - array of edges
    ! _y - helper arrays for coordinate construction

    if (.not.use_rxv) ERROR STOP 'require rxv'

    t = trel + toffset

    do i=1,nstar
       star_v(:,i) = star_huntpol(t, i)
    end do

    m_v(:)   = star_v(IMASS, :)
    s_v(:)   = star_v(IRADI, :)
    i_v(:,:) = star_v(IINER, :)
    k_v(:)   = star_v(IAPSI, :)
    t_v(:)   = star_v(ITAUQ, :)

    ! use q < 0 as flag for containing -1/Q, tau for q > 0

    r_j(:,:) = reshape(y(1            :  ndim*norb),         (/ndim, norb /))
    v_j(:,:) = reshape(y(1+  ndim*norb:2*ndim*norb),         (/ndim, norb /))
    j_v(:,:) = reshape(y(1+2*ndim*norb:ndim*(2*norb+nstar)), (/ndim, nstar/))
    if (use_ori) then
       if (use_orq) then
          orq_v(:,:) = reshape(y(1+2*nvorb+nvstar:2*nvorb+nvstar+nqstar), &
               (/nquat, nstar/))
          roq_v(:,:) = 0.d0
       else
          ori_v(:,:) = reshape(y(1+2*nvorb+nvstar:2*(nvorb+nvstar)), &
               (/ndim, nstar/))
          rot_v(:,:) = 0.d0
       endif
    end if

    mi_v(:) = 1.d0 / m_v(:)

    s5k_v(:) = s_v(:)**5 * k_v(:)
    m2_v(:) = m_v(:)**2

    mp_e(:) = m_v(ies(1,:)) * m_v(ies(2,:))
    m_e(:)  = m_v(ies(1,:)) + m_v(ies(2,:))

    mm2_e(:) = 1.d0 / m_e(:)**2

    ! convert from Jacobi coordinates to real relative vectors
    if (use_qp) then

       mi_d(1) = 1.d0 / m_e(1)
       mi_d(2) = 1.d0 / m_e(6)

       mc_y(1,1) =  m_v(2) * mi_d(1)
       mc_y(2,1) = -m_v(1) * mi_d(1)
       mc_y(1,2) = -m_v(4) * mi_d(2)
       mc_y(2,2) =  m_v(3) * mi_d(2)

       ! do i=1,2
       !    do j=1,2
       !       k = 1 - i + 2 * j ! use map
       !       mc_y(i,j) = float((3 - 2 * j) * (3 - 2 * i), kind=real64) * mi_d(j)
       !    enddo
       ! enddo

       do i=1,2
          do j=1,2
             rc_y(:,i,j) = mc_y(i,j) * r_j(:,j)
             vc_y(:,i,j) = mc_y(i,j) * v_j(:,j)
          enddo
       enddo

       r_e(:,1) = r_j(:,1)
       r_e(:,6) = r_j(:,2)

       r_e(:,2) = r_j(:,3) + rc_y(:,1,1) + rc_y(:,1,2)
       r_e(:,3) = r_j(:,3) + rc_y(:,1,1) + rc_y(:,2,2)

       r_e(:,4) = r_j(:,3) + rc_y(:,2,1) + rc_y(:,1,2)
       r_e(:,5) = r_j(:,3) + rc_y(:,2,1) + rc_y(:,2,2)

       v_e(:,1) = v_j(:,1)
       v_e(:,6) = v_j(:,2)

       v_e(:,2) = v_j(:,3) + vc_y(:,1,1) + vc_y(:,1,2)
       v_e(:,3) = v_j(:,3) + vc_y(:,1,1) + vc_y(:,2,2)

       v_e(:,4) = v_j(:,3) + vc_y(:,2,1) + vc_y(:,1,2)
       v_e(:,5) = v_j(:,3) + vc_y(:,2,1) + vc_y(:,2,2)

       ! do i=1,2
       !    do j=1,2
       !       k = 2 * i + j - 1 ! use map
       !       r_e(:,k) = r_j(:,3) + rc_y(:,i,1) + rc_y(:,j,2)
       !       v_e(:,k) = v_j(:,3) + vc_y(:,i,1) + vc_y(:,j,2)
       !    enddo
       ! enddo

    else

       m_12 = m_e(1)
       mi_12 = 1.d0 / m_12

       m_123 = m_12 + m_v(3)
       mi_123 = 1.d0 / m_123

       m_1234 = m_123 + m_v(4)

       mc_y(1,1) = + m_v(2) * mi_12
       mc_y(1,2) = - m_v(1) * mi_12

       mc_y(2,1) = + m_v(3) * mi_123
       mc_y(2,2) = - m_12   * mi_123

       do i=1,2
          rc_y(:,1,i) = mc_y(1,i) * r_j(:,1)
          rc_y(:,2,i) = mc_y(2,i) * r_j(:,2) + r_j(:,3)

          vc_y(:,1,i) = mc_y(1,i) * v_j(:,1)
          vc_y(:,2,i) = mc_y(2,i) * v_j(:,2) + v_j(:,3)
       enddo

       r_e(:,1) = r_j(:,1)
       r_e(:,2) = r_j(:,2)    + rc_y(:,1,1)
       r_e(:,4) = r_j(:,2)    + rc_y(:,1,2)
       r_e(:,3) = rc_y(:,2,1) + rc_y(:,1,1)
       r_e(:,5) = rc_y(:,2,1) + rc_y(:,1,2)
       r_e(:,6) = rc_y(:,2,2)

       v_e(:,1) = v_j(:,1)
       v_e(:,2) = v_j(:,2)    + vc_y(:,1,1)
       v_e(:,4) = v_j(:,2)    + vc_y(:,1,2)
       v_e(:,3) = vc_y(:,2,1) + vc_y(:,1,1)
       v_e(:,5) = vc_y(:,2,1) + vc_y(:,1,2)
       v_e(:,6) = vc_y(:,2,2)

    endif

    do l=1,nedge
       rn_e(l) = norm2(r_e(:,l))
    enddo

    if (use_ori) then
       do concurrent (i=1:nstar)

          ! alternate approach if use i_x == i_y == i_z
          ! if (abs(i_v(1,i)-i_v(2,i)) + abs(i_v(2,i)-i_v(3,i)) > 1.d-12 * i_v(3,i)) then

          ! here we assume i_x == i_y == 0 for spherical bodies
          ! and i_z == 0 if just point mass w/o spin

          if (i_v(3,i) == 0.d0) then
             w_v(:,i) = 0.d0
             lpqm_v(i) = .false.
             cycle
          endif

          if (i_v(1,i) > 0.d0) then
             if (use_orq) then
                a_v(:,:,i) = rotation_f2A(orq_v(:,i))
             else
                a_v(:,:,i) = rotation_w2A(ori_v(:,i))
             endif
             w_v(:,i) = matmul(a_v(:,:,i), matmul(j_v(:,i), a_v(:,:,i)) / i_v(:,i))
             lpqm_v(i) = .true.
          else
             w_v(:,i) = j_v(:,i) / i_v(3,i)
             lpqm_v(i) = .false.
          endif

          if (.not.(lpqm_v(i) .or. use_ors)) &
               cycle

          if (use_orq) then
             roq_v(:,i) =  rotation_dqdt(orq_v(:,i), w_v(:,i))
          else
             rot_v(:,i) =  rotation_dodt(ori_v(:,i), w_v(:,i))
          endif

          ! TODO - modify rot(:,:) for GR/SR time delay

          if (.not.(lpqm_v(i) .and. (use_pqm .or. use_pqq))) &
               cycle

          q_v(:,i) = ONETHIRD * sum(i_v(:,i)) - i_v(:,i)
          do k=1, ncon
             l = ise(k,i)
             j = isc(k,i)
             qr_c(:,j,l) = matmul(a_v(:,:,i), q_v(:,i) * matmul(r_e(:,l), a_v(:,:,i)))
             qrr_c(j,l) = sum(r_e(:,l) * qr_c(:,j,l))
          enddo
       enddo
    else
       do concurrent (i=1:nstar)
          if (i_v(3,i) > 0.d0) then
             w_v(:,i) = j_v(:,i) / i_v(3,i)
          else
             w_v(:,i) = 0.d0
          endif
          lpqm_v(i) = .false.
       end do
    end if
    wn2_v(:) = sum(w_v(:,:)**2, dim=1)

    ! fail on star interaction
    if (has_interact) then
       do l=1,nedge
          if (rn_e(l) < interact(ies(1,l),ies(2,l))) then
             if (verb_interact) then
                write(6, "(' [direct4f] object interaction ',i1,' + ',i1)") ies(:,l)
             endif
             yp(:) = ieee_value(0.d0,ieee_quiet_nan)
             status = STATUS_INTERACT
             status_stars = ies(:,l)
             terminate_integration = .true.
             return
          endif
       end do
    endif

    ! fail on star collisions
    do l=1,nedge
       if (rn_e(l) < s_v(ies(1,l)) + s_v(ies(2,l))) then
          if (verb_interact) then
             write(6, "(' [direct4f] object collision ',i1,' + ',i1)") ies(:,l)
          endif
          yp(:) = ieee_value(0.d0,ieee_quiet_nan)
          status = STATUS_COLLIDE
          status_stars = ies(:,l)
          terminate_integration = .true.
          return
       endif
    end do

    ! fail on star escape
    if (has_cutoff) then
       do j=1, norb
          if (norm2(r_j(:,j)) > cutoff(j)) then
             if (verb_interact) then
                write(6, "(' [direct4f] orbit escape ',i1)") j
             endif
             yp(:) = ieee_value(0.d0,ieee_quiet_nan)
             status = STATUS_ESCAPE
             status_stars = (/j, 0/)
             terminate_integration = .true.
             return
          endif
       enddo
    endif

    ! compute forces

    rn2_e(:)  = rn_e(:)**2
    rnm1_e(:) = 1.d0 / rn_e(:)
    rnm2_e(:) = rnm1_e(:)**2
    rnm3_e(:) = rnm1_e(:) * rnm2_e(:)
    rnm4_e(:) = rnm2_e(:)**2
    rnm5_e(:) = rnm4_e(:) * rnm1_e(:)
    rnm7_e(:) = rnm4_e(:) * rnm3_e(:)
    rnm8_e(:) = rnm4_e(:)**2

    ! used for qd or tf
    if (use_sd .or. use_tf) then
       do l=1,nedge
          do j=1,2
             i = ies(j,l)
             wdr_c(j,l) = dot_product(w_v(:,i), r_e(:,l))
             wxr_c(:,j,l) = cross_product(w_v(:,i), r_e(:,l))
          end do
       end do
    end if

    ! direct gravity
    do l=1,nedge
       v(:) = ((- GRAV) * rnm3_e(l) * mp_e(l)) * r_e(:,l)
       fng_e(:,l) = v(:)
       f_e(:,l) = v(:)
    enddo

    if (use_pn) then
       do l=1, nedge
          eta = mp_e(l) * mm2_e(l)
          phi = GRAV * rnm1_e(l) * m_e(l)
          vn2 = sum(v_e(:,l)**2)
          vr = dot_product(v_e(:,l), r_e(:,l))
          g = Gcm2 * mp_e(l) * rnm3_e(l)
          v(:) =  &
               ((( 4.d0 + 2.d0 * eta) * phi + &
                 (-1.d0 - 3.d0 * eta) * vn2 + &
                 (       1.5d0 * eta) * vr**2 * rnm2_e(l)) * g) * r_e(:,l) + &
                (( 4.d0 - 2.d0 * eta) * vr * g) * v_e(:,l)
          f_e(:,l) = f_e(:,l) + v(:)
          fpn_e(:,l) = v(:)

          ggr_e(l) = CLIGHTm2 * ( &
              (1.d0 - 3.d0 * eta) * 0.5d0 * vn2 &
            + (3.d0 +        eta) * phi)
          fgr_e(l) = 1.d0 + ggr_e(l)
       enddo
    else
       ggr_e(:) = 0.d0
       fgr_e(:) = 1.d0
       fpn_e(:,:) = 0.d0
    endif

    torque_v(:,:) = 0.d0

    ! permanent quadrupole-quadrupole interaction

    if (use_ori .and. use_pqq .and. any(lpqm_v)) then
       error stop '[direct4f] use_pqq not implemented'
    endif

    ! permanent quadrupole-monopole interaction

    if (use_ori .and. use_pqm .and. any(lpqm_v)) then
       do l=1,nedge
          g = G3 * rnm2_e(l) * rnm3_e(l)
          do j=1,ncorn
             i = ies(j,l)
             if (.not.lpqm_v(i)) then
                fpq_c(:,j,l) = 0.d0
                tpq_c(:,j,l) = 0.d0
                cycle
             endif
             ir = ies(3-j,l)
             f = g * m_v(ir)
             fpq_c(:,j,l) = f * ( &
                qr_c(:,j,l) - 2.5d0 * rnm2_e(l) * qrr_c(j,l) * r_e(:,l))
             tpq_c(:,j,l) = f * cross_product(qr_c(:,j,l), r_e(:,l)) * fgr_e(l)

             f_e(:,l) = f_e(:,l) + fpq_c(:,j,l)
             torque_v(:,i) = torque_v(:,i) + tpq_c(:,j,l)
          end do
       end do
    else
       fpq_c(:,:,:) = 0.d0
       tpq_c(:,:,:) = 0.d0
    endif

    if (use_sd .or. use_td) then
       do l=1,nedge
          if (use_td) then
             h = (-12.d0 * GRAV) * rnm1_e(l)
          endif
          do j=1,2
             i = ies(j,l)
             ir = ies(3-j,l)
             f = s5k_v(i) * m_v(ir) * rnm7_e(l)
             if (use_td) then
                ftd_c(:,j,l) = (f * h * m_v(ir)) * r_e(:,l)
                f_e(:,l) = f_e(:,l) + ftd_c(:,j,l)
             else
                ftd_c(:,j,l) = 0.0d0
             endif

             if (use_sd) then
                g = -2.d0 * rn2_e(l) * wdr_c(j,l) * f
                fsd_c(:,j,l) = &
                     (f * (5.d0 * wdr_c(j,l)**2 - rn2_e(l) * wn2_v(i))) * r_e(:,l) &
                     +  g * w_v(:,i)
                f_e(:,l) = f_e(:,l) + fsd_c(:,j,l)

                torque_v(:,i) = torque_v(:,i) + (g * fgr_e(l)) * wxr_c(:,j,l)
             else
                fsd_c(:,j,l) = 0.0d0
             endif
          end do
       end do
    else
       ftd_c(:,:,:) = 0.0d0
       fsd_c(:,:,:) = 0.0d0
    endif

    if (use_td3) then
       do i=1, nstar
          do j1=1, nstar1
             j = iij(i,j1)
             ij = ipe(i,j)
             rij(:) = r_e(:,ij)
             do k1=1, nstar2
                k = iijk(i,j1,k1)

                ! DEBUG
                if ((i==j).or.(j==k).or.(k==i)) STOP 'index bug'

                ik = ipe(i,k)
                rik(:) = r_e(:,ik) ! all signs of r cancel out for f_e
                if (j < k) then ! no need to redo symmetric (dot) products
                   g = dot_product(rij(:), rik(:))
                   gjk_x(j,k) = g
                   f = m_v(j) * m_v(k)
                   mjk_x(j,k) = f
                else
                   g = gjk_x(k,j)
                   f = mjk_x(k,j)
                endif
                f = f * s5k_v(i) * GRAV * rnm3_e(ij) * rnm5_e(ik)
                h = 6.d0 * f * g * rnm2_e(ij)
                v(:) = &
                     ((3.d0 - 15.d0 * g**2 * rnm2_e(ij) * rnm2_e(ik)) * f &
                     ) * rik(:) + &
                     h * rij(:)
                f_e(:, ik) = f_e(:, ik) + v(:)
                ftd3_d(:,i,j1,k1) = v(:)

                ! DEVEL - only record net torque, cancels out due to symmetry
                if (.not. use_pn) cycle
                if (j < k) then ! no need to redo antisymmetric cross product
                   v(:) = (h * (ggr_e(ik) - ggr_e(ij))) * &
                        cross_product(rij(:), rik(:))
                   torque_v(:,i) = torque_v(:,i) + v(:)
                endif
             enddo
          enddo
       enddo
    else
       ftd3_d(:,:,:,:) = 0.0d0
    endif

    if (use_tf .and. any(t_v /= 0.d0)) then
       do l=1, nedge
          w(:) = cross_product(r_e(:,l), v_e(:,l)) * rnm2_e(l)
          g = 12.d0 * GRAV * rnm8_e(l)
          h = -2.d0 * dot_product(r_e(:,l), v_e(:,l)) * rnm2_e(l)
          if (any(t_v(ies(1:ncorn,l)) < 0.d0)) then
             if (use_rxv) then
                nm1 = 1.d0 / max(norm2(w), &
                     omegalim * sqrt(GRAV * m_e(l) * rnm3_e(l)))
             else
                Gm = GRAV * m_e(l)
                phi = -Gm * rnm1_e(l)
                eps = min(0.5d0 * sum(v_e(:,l)**2) + phi, epslim * phi)
                an = -0.5d0 * Gm / eps
                nm1 = sqrt(an**3 / Gm)
             endif
          end if
          do j=1, ncorn
             i = ies(j,l)
             f = t_v(i)
             if (f == 0.0d0) then
                ftf_c(:,j,l) = 0.d0
                continue
             else if (f < 0.0d0) then
                f = -0.5d0 * f * nm1
             endif
             ir = ies(3-j,l)
             f = f * g * s5k_v(i) * m2_v(ir)
             v(:) = f * (wxr_c(:,j,l) - v_e(:,l) + h * r_e(:,l))

             ftf_c(:,j,l) = v(:)
             f_e(:,l) = f_e(:,l) + v(:)

             f = f * fgr_e(l)
             v(:) = (f * wdr_c(j,l)) * r_e(:,l) &
                  + (f * rn2_e(l)) * (w(:) - w_v(:,i))
             torque_v(:,i) = torque_v(:,i) + v(:)
          end do
       end do
    else
       ftf_c(:,:,:) = 0.d0
    endif

    ! Note that f_e points in opposite direction as r_e as per ML02

    acc_j(:,1) = &
         + ( mi_v(1)  +  mi_v(2)) * f_e(:,1) &
         +   mi_v(1)  * (f_e(:,2) + f_e(:,3)) &
         + (-mi_v(2)) * (f_e(:,4) + f_e(:,5))
    if (use_qp) then
       acc_j(:,2) = &
            + ( mi_v(3)  +  mi_v(4)) * f_e(:,6) &
            + (-mi_v(3)) * (f_e(:,2) + f_e(:,4)) &
            +   mi_v(4)  * (f_e(:,3) + f_e(:,5))
       acc_j(:,3) = (mi_d(1) + mi_d(2)) * sum(f_e(:,2:5), dim=2)
    else
       acc_j(:,2) = (m_123 * mi_12 * mi_v(3)) * (f_e(:,2) + f_e(:,4)) &
            + mi_12 * (f_e(:,3) + f_e(:,5)) + (-mi_v(3)) * f_e(:,6)
       acc_j(:,3) = (m_1234 * mi_123 * mi_v(4)) * (f_e(:,3) + f_e(:,5) + f_e(:,6))
    endif

    yp( 1: 9) = reshape(v_j,      (/ndim*norb/))
    yp(10:18) = reshape(acc_j,    (/ndim*norb/))
    yp(19:30) = reshape(torque_v, (/ndim*nstar/))
    if (use_ori) then
       if (use_orq) then
          yp(31:46) = reshape(roq_v, (/nqstar/))
       else
          yp(31:42) = reshape(rot_v, (/nvstar/))
       endif
    endif

    if (local_store) then

       ! actual derivative data
       local_data(  1: size(yp,1)) = yp(:)

       ! basic forces and torques
       local_data( 47: 82) = reshape(fsd_c,  (/ndim*ncorn*nedge/)) !36
       local_data( 83:118) = reshape(ftd_c,  (/ndim*ncorn*nedge/)) !36
       local_data(119:154) = reshape(ftf_c,  (/ndim*ncorn*nedge/)) !36
       local_data(155:226) = reshape(ftd3_d, (/ndim*nstar*nstar1*nstar2/)) !72
       local_data(227:244) = reshape(fpn_e,  (/ndim*nedge/)) !18
       local_data(245:262) = reshape(fng_e,  (/ndim*nedge/)) !18

       ! permanent quadrupole momenent
       local_data(263:298) = reshape(fpq_c,  (/ndim*ncorn*nedge/)) !36
       local_data(299:334) = reshape(tpq_c,  (/ndim*ncorn*nedge/)) !36
    endif

  end function directpf

end module polyadic