!
! ----------------- tuben0_return -------------
! return the neutral density in some tubes from xplasma -- the tube data begins at 'zs0'
! from the tube start for 'nsegtube(itube)' segments with segment length 
! 'deltaseg(itube)'
!
subroutine tuben0_return(ivec, ivtube, ivwant, mn0, ztuben0, mdot, ztubedot, zs0, ierr)
  use iso_c_binding, only: fp => c_double
  use tube_mod, only: deptube_xplasma_name
  use xplasma_obj_instance, only: s, xplasma, xplasma_blkbxid, xplasma_blackbox_retrieve, xplasma_get_item_info
  implicit none

  integer,       intent(in)  :: ivec              ! number of tubes wanted
  integer,       intent(in)  :: ivtube(*)         ! tube indices being fectched [ivec]
  integer,       intent(in)  :: ivwant(*)         ! number of segments desired per tube -- checked against nsegtube [ivec]
  integer,       intent(in)  :: mn0               ! first dimension of ztuben0
  real(fp),        intent(out) :: ztuben0(mn0,*)    ! neutral particle density (#/cm^3) [mn0,ivec]
  integer,       intent(in)  :: mdot              ! first dimension of ztubedot
  real(fp),        intent(out) :: ztubedot(mdot,*)  ! cos(angle_beam_to_tube) (-)  [mdot,ivec]
  real(fp),        intent(out) :: zs0(*)            ! distance from tube start where data begins  [ivec]
  integer,       intent(out) :: ierr              ! nonzero on error
  
  integer,save :: id=0   ! xplasma integer id, saved between calls since it might not change
  character*64 :: sname  ! name stored under id in xplasma

  integer :: ii     ! error flag
  integer :: iv     ! loop index
  integer :: iwant  ! number of segments in tube
  integer :: itube  ! tube index to fetch
  integer :: kt     ! number of tubes stored in xplasma
  integer :: ks     ! number of segments in desired tube
  integer :: k      ! index in real array of tube data

  integer, pointer :: iptr(:)   ! array of indices pointing to start of each tube
  real(fp),  pointer :: rptr(:)   ! packed tube data, for each tube  s0,tuben0(1:nseg),tubedot(1:nseg)

  ierr=0
  if (ivec<=0) return

  ! -- get id --
  if (id>0) then
     call xplasma_get_item_info(s,id,ii,name=sname) ! check to see if we already have the id
     if (ii/=0) then
        id=0
     else if (sname /= deptube_xplasma_name) then
        id=0                                        ! do not test if ii/=0 because sname may not have been set
     end if
  end if

  if (id<=0) then
     call xplasma_blkbxid(s,deptube_xplasma_name,id)
  end if

  if (id<=0) then
     print *, "?tuben0_return: error fetching xplasma blackbox id for '"//deptube_xplasma_name//"'"
     ierr=1 ; return
  end if
  
  ! -- fetch pointers --
  call xplasma_blackbox_retrieve(s, id, ierr, ia_ptr=iptr, r8a_ptr=rptr)
  
  if (ierr/=0) then
     print *, '?tuben0_return: error fetching deptube pointers from black box '//deptube_xplasma_name
     ierr=1 ; return
  end if

  ! -- extract data --
  kt = size(iptr)-1
  
  do iv = 1, ivec
     itube = ivtube(iv)
     iwant = ivwant(iv)

     if (itube<1 .or. itube>kt) then
        print *, '?tuben0_return: tube index is out of bounds'
        print *, '    itube  = ', itube
        print *, '    ntubes = ', kt, ' in xplasma'
        ierr=1 ; return
     end if
     
     k  = iptr(itube)
     ks = (iptr(itube+1)-k)/2
     
     if (iwant/=ks .or. iwant>mn0 .or. iwant>mdot) then
        print *, '?tuben0_return: wrong number of segments requested'
        print *, '    iwant  = ', iwant
        print *, '    nsegs  = ', ks, ' in xplasma for this tube'
        print *, ' dim(n0)   = ', mn0
        print *, ' dim(dot)  = ', mdot
        ierr=1 ; return
     end if
     
     zs0(iv)              = rptr(k)             ! tubes0(itube)
     ztuben0(1:iwant,iv)  = rptr(k+1:k+ks)      ! tuben0(1:iwant,itube)
     ztubedot(1:iwant,iv) = rptr(k+ks+1:k+2*ks) ! tubedot(1:iwant,itube)
  end do
end subroutine tuben0_return

!
! ----------------- tuben0_return_state -------------
! return the neutral density in one tube by fetching from the state variables
!
subroutine tuben0_return_state(itube, iwant, ztuben0, ztubedot, zs0, ierr)
  use iso_c_binding, only: fp => c_double
  use tube_mod
  implicit none

  integer, intent(in)  :: itube        ! tube index
  integer, intent(in)  :: iwant        ! number of segments desired -- checked against nsegtube
  real(fp),  intent(out) :: ztuben0(*)   ! neutral particle density #/cm^3
  real(fp),  intent(out) :: ztubedot(*)  ! cos(angle_beam_to_tube)
  real(fp),  intent(out) :: zs0          ! distance from tube start where data begins
  integer, intent(out) :: ierr         ! nonzero on error
  
  ierr = 0
  
  if (itube<1 .or. itube>nbtube) then
     write(nonlin, *) '?tuben0_return: tube index is out of bounds'
     write(nonlin, *) '    itube  = ', itube
     write(nonlin, *) '    ntubes = ', nbtube
     ierr=1 ; return
  end if
  
  if (iwant /= nsegtube(itube)) then
     write(nonlin, *) '?tuben0_return: wrong number of segments requested'
     write(nonlin, *) '    iwant = ', iwant
     write(nonlin, *) '    nsegs = ', nsegtube(itube)
     ierr=1 ; return
  end if

  zs0               = tubes0(itube)
  ztuben0(1:iwant)  = tuben0(1:iwant,itube)
  ztubedot(1:iwant) = tubedot(1:iwant,itube)
end subroutine tuben0_return_state

