module stardata

  use typedef, only: &
       real64, int32, int64

  implicit none

  save

  type stardatatype
     real(kind=real64), dimension(:), allocatable :: &
          t
     real(kind=real64), dimension(:,:), allocatable :: &
          data
     integer(kind=int32) :: &
          n, m, i
  end type stardatatype

  type(stardatatype), dimension(:), allocatable :: &
       stars

  real(kind=real64), dimension(:), allocatable  :: &
       t0

  integer(kind=int64) :: &
       stardata_id

  real(kind=real64) :: &
       toffset = 0.d0

  integer(kind=int32), parameter :: &
       NDATA = 8

  integer(kind=int32), parameter :: &
       ITIME = 1, &
       IMASS = 2, &
       IRADI = 3, &
       IINEX = 4, &
       IINEY = 5, &
       IINEZ = 6, &
       IAPSI = 7, &
       ITAUQ = 8

  integer(kind=int32), dimension(3), parameter :: &
       IINER = (/IINEX, IINEY, IINEZ /)

contains

  ! ======================================================================
  ! IDs are being used to track data for different runs using same module

  subroutine set_stardata_id(id_)

    use typedef, only: &
         int64

    implicit none

    save

    integer(kind=int64), intent(in) :: &
         id_

    stardata_id = id_

  end subroutine set_stardata_id

  ! ======================================================================

  subroutine get_stardata_id(id_)

    use typedef, only: &
         int64

    implicit none

    save

    integer(kind=int64), intent(out) :: &
         id_

    id_ = stardata_id

  end subroutine get_stardata_id

  ! ======================================================================

  subroutine random_stardata_id(id_)

    use typedef, only: &
         int64, real64

    implicit none

    save

    integer(kind=int64), intent(out) :: &
         id_
    real(kind=real64) :: &
         r

    call random_number(r)

    stardata_id = nint(r * 9223372036854775800.d0, kind=int64)

    id_ = stardata_id

  end subroutine random_stardata_id

  ! ======================================================================

  subroutine set_toffset(toffset_)

    use typedef, only: &
         real64

    implicit none

    save

    real(kind=real64), intent(in) :: &
         toffset_

    toffset = toffset_

  end subroutine set_toffset

  ! ======================================================================

  subroutine add_stardata(data_, t_)

    use typedef, only: &
         real64, int32

    implicit none

    save

    real(kind=real64), intent(in), dimension(:,:)  :: &
         data_
    real(kind=real64), intent(in)  :: &
         t_

    type(stardatatype), dimension(:), allocatable :: &
         stars_
    real(kind=real64), dimension(:), allocatable :: &
         t0_
    integer(kind=int32) :: &
         i, n

    if (.not.allocated(stars)) then
       allocate(stars(1), t0(1))
       n = 1
    else
       n = size(stars, 1) + 1
       allocate(stars_(n), t0_(n))
       do i=1, n-1
          call move_alloc(stars(i)%t, stars_(i)%t)
          call move_alloc(stars(i)%data, stars_(i)%data)
          stars_(i)%n = stars(i)%n
          stars_(i)%m = stars(i)%m
          stars_(i)%i = stars(i)%i
          t0_(i) = t0(i)
       enddo
       call move_alloc(stars_, stars)
       call move_alloc(t0_, t0)
    endif
    stars(n)%i = -1
    stars(n)%n = size(data_, 1)
    stars(n)%m = size(data_, 2)
    allocate(stars(n)%t(stars(n)%n))
    allocate(stars(n)%data(stars(n)%m, stars(n)%n))
    stars(n)%t(:) = data_(:, ITIME)
    stars(n)%data(:,:) = transpose(data_(:,:))

    t0(n) = t_

    if (n > 1) then
       if (stars(n)%m /= stars(n-1)%m) error STOP '[ADD_STARDATA] dimension error m'
    endif

  end subroutine add_stardata

  ! ======================================================================

  subroutine clear_stardata()

    use typedef, only: &
         int32

    implicit none

    integer(kind=int32) :: &
         i

    if (.not.allocated(stars)) return

    do i = 1, size(stars, 1)
       deallocate(stars(i)%t, stars(i)%data)
    enddo
    deallocate(stars)
    deallocate(t0)

    toffset = 0.d0
    stardata_id = -1

  end subroutine clear_stardata

  ! ======================================================================

  function n_stardata() result(count)

    use typedef, only: &
         int32

    implicit none

    integer(kind=int32) :: &
         count

    if (.not.allocated(stars)) then
       count = 0
    else
       count = size(stars, 1)
    endif

  end function n_stardata

  ! ======================================================================

  function star_huntpol(t, j) result(val)

    use typedef, only: &
         int32, real64

    use hunting, only: &
         huntpoli

    implicit none

    real(kind=real64), intent(in) :: &
         t
    integer(kind=int32), intent(in) :: &
         j
    real(kind=real64), dimension(stars(j)%m) :: &
         val

    val = huntpoli(t + t0(j), stars(j)%t, stars(j)%n, stars(j)%data, stars(j)%m, stars(j)%i)

  end function star_huntpol

  ! ======================================================================

  subroutine get_star_data(tt, dd)

  use typedef, only: &
       real64, int32

  implicit none

  real(kind=real64), intent(out), dimension(:,:,:) :: &
       dd
  real(kind=real64), intent(in), dimension(:) :: &
       tt

  integer(kind=int32) :: &
       nt, nd, ns, i, j, n
  real(kind=real64), dimension(:), allocatable :: &
       t

  if (.not.allocated(stars))     error stop '[get_star_data] no stardata available'
  n = size(stars, 1)

  nd = size(dd, 1)
  nt = size(dd, 2)
  ns = size(dd, 3)

  ! could use ierr
  if (size(tt, 1) /= nt)         error stop '[get_star_data] tt, dd diminsion mismatch'
  if (ns /= n)                   error stop '[get_star_data] dd(:,:,x) diminsion equal stardata required'

  allocate(t(nt))
  do i=1,nt
     t(i) = toffset + tt(i)
  enddo

  do j=1, n
     do i=1, nt
        dd(:,i,j) = star_huntpol(t(i), j)
     enddo
  end do

  deallocate(t)

  end subroutine get_star_data

end module stardata