module unary

  implicit none

contains

  ! =======================================================================
  ! single star - just spin/tumble
  ! =======================================================================

  function direct1(trel, y) result(yp)

    use typedef, only: &
         int32, real64

    use stardata, only: &
         IINER, &
         toffset, stars, &
         star_huntpol

    use deriv_data, only: &
         local_store, local_data

    use flags_data, only: &
         use_ori, use_orq, use_ors

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

    implicit none

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

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

    integer(kind=int32), parameter :: &
         nstar = 1, &
         ndim  = 3, &
         nquat = 4

    ! local variables

    real(kind=real64), dimension(nquat) :: &
         orq, roq
    real(kind=real64), dimension(ndim) :: &
         ori, j, w, i, &
         torque, rot
    real(kind=real64), dimension(stars(1)%m) :: &
         star
    real(kind=real64) :: &
         t
    real(kind=real64), dimension(ndim, ndim) :: &
         A
    logical :: &
         lpqm

    t = trel + toffset

    star(:) = star_huntpol(t, 1)
    i(:) = star(IINER)

    j(:) = y(1:3)
    if (use_ori) then
       if (use_orq) then
          orq(:) = y(4:7)
       else
          ori(:) = y(4:6)
       endif
    end if

    if (use_ori) then

       ! alternate approach if use i_x == i_y == i_z
       ! if (abs(i1(1,i)-i1(2,i)) + abs(i1(2,i)-i1(3,i)) > 1.d-12 * I1(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(3) == 0.d0) then
          w(:) = 0.d0
          lpqm = .false.
          goto 1000
       endif

       if (i(1) > 0.d0) then
          if (use_orq) then
             A(:,:) = rotation_f2A(orq(:))
          else
             A(:,:) = rotation_w2A(ori(:))
          endif
          w(:) = matmul(A(:,:), matmul(j(:), A(:,:)) / i(:))
          lpqm = .true.
       else
          w(:) = j(:) / i(3)
          lpqm = .false.
       endif

       if (.not.(lpqm .or. use_ors)) &
            goto 1000

       if (use_orq) then
          roq(:) = rotation_dqdt(orq(:), w(:))
       else
          rot(:) = rotation_dodt(ori(:), w(:))
       endif
       goto 2000
    endif

1000 continue
    if (use_orq) then
       roq(:) = 0.d0
    else
       rot(:) = 0.d0
    endif

2000 continue
    torque = 0.d0

    yp(1:3) = torque(:)
    if (use_ori) then
       if (use_orq) then
          yp(4:7) = roq(:)
       else
          yp(4:6) = rot(:)
       endif
    endif

    if (local_store) then
       local_data( 1:size(yp,1)) = yp(:)
    endif

  end function direct1

end module unary