module xstraln_r8track
 
  SAVE
 
  ! these globals can be increased on the fly...
 
  integer :: maxpts=1000  ! standard grid array size for structure
  integer :: maxvac=10    ! standard vac segment array size for structure
 
  type :: r8_track
 
     integer :: maxpts    ! max dimension of intercept and vector arrays
     integer :: maxvac    ! max dimensino of vacuum segment arrays
 
     integer :: guard     ! =12345 if track arrays are allocated.
 
     integer :: n_int     ! number of rho and theta surfaces intercepted
 
     ! in cartesian coordinates [x,y,z] [m,m,m]
     REAL*8, dimension(3) :: pta,ptb   ! track start & stop points (as used)
     REAL*8, dimension(3) :: gxyz      ! (ptb-pta)/|ptb-pta|  [no units]
 
     ! data on zones crossed
 
     integer :: nzcrtk               ! number of zones (0 means missed plasma)
     integer, dimension(:), allocatable :: nznttk   ! radial zone indices
     integer, dimension(:), allocatable :: nznbtk   ! beam 2d grid zone indices
     REAL*8, dimension(:), allocatable :: xlngtk      ! segment lengths [m]
 
     ! data on vacuum segments (parts of track that are in vacuum; track
     ! has left plasma but will re-enter).  xltkvac(0) gives length
     ! from pta to first entry into plasma
 
     integer :: nvacseg              ! number of vacuum segments
     integer, dimension(:), allocatable :: nztkvac  ! preceding plasma segment
     REAL*8, dimension(:), allocatable :: xltkvac     ! length of segment
 
     ! data on intercepted surfaces (0's indicate intercept of other type)
     ! distances measured from the *internal* start point (pta).  Caution,
     ! this can be different than what the user gave, because the software
     ! may shorten the segment so that no part of it is outside the xstraln
     ! mapped region [Rmin,Rmax]x[Zmin,Zmax].  If either endpoints are moved
     ! the corresponding quantities (s_to_pta) (s_from_ptb) will be set.
 
     integer, dimension(:), allocatable :: irho_int        ! rho surface
     integer, dimension(:), allocatable :: ibr_int,ith_int ! theta surface
     REAL*8, dimension(:), allocatable :: s_int    ! distances to intercepts
     REAL*8, dimension(:), allocatable :: s_zon    ! distances to segment ctrs
     REAL*8 :: s_to_pta               ! distance from user start point to pta
     REAL*8 :: s_from_ptb             ! distance to user end point from ptb
 
     ! notes: nzcrtk = n_int-1
     ! s_int has n_int elements defined; s_zon has nzcrtk elements defined.
 
     !-----------------------------------------------
     ! internal tracker (adaptive) grid-- more tightly spaced
     ! where target surfaces are close together, and near tangency points
 
     integer :: npts    ! gridsize
     REAL*8, dimension(:), allocatable :: svec   ! grid [m] from 0.0 to |ptb-pta|
     REAL*8, dimension(:,:), allocatable :: spkg ! (cf pspline genxpkg, from svec)
 
     REAL*8, dimension(:), allocatable :: rvec   ! R [m]
     REAL*8, dimension(:), allocatable :: zvec   ! Z [m]
     REAL*8, dimension(:), allocatable :: phivec ! phi [rad] no branch cut jumps
 
     REAL*8, dimension(:), allocatable :: rhovec ! rho (norm. sqrt tor. flux)
     REAL*8, dimension(:), allocatable :: thvec  ! theta [rad] no branch cut jumps
 
     REAL*8, dimension(:), allocatable :: Bxvec  ! Bx [T]
     REAL*8, dimension(:), allocatable :: Byvec  ! By [T]
     REAL*8, dimension(:), allocatable :: Bzvec  ! Bz [T]
 
  end type r8_track
 
  contains
 
    subroutine xstraln_r8track_expand(i)
 
      type (r8_track) :: i
 
      ! increase size of arrays (implies re-initialization of structure)
 
      if(i%maxpts.ge.maxpts) then
         maxpts = 2*maxpts
         maxvac = 2*maxvac
      endif
 
      call xstraln_r8track_init(i,(2*i%maxpts),.FALSE.)
 
    end subroutine xstraln_r8track_expand
 
    subroutine xstraln_r8track_init(i,isize,iforce)
 
      type (r8_track) :: i
      integer, intent(in) :: isize   ! desired size (or if 0 use the global)
      logical, intent(in) :: iforce
 
      ! initialize track data structure
      ! if (iforce) force the initialization; otherwise only do it if
      ! the array size variable does not match the saved global size.
      ! if i%guard indicates a prior allocation, free that memory
      ! with a dealloc call as well.
 
      ! (iforce) is a good idea on first initialization of a new object.
 
      integer jsize
 
      !-----------------------------
 
      if(iforce) then
         i%maxpts=0
         i%maxvac=0
         i%guard=0
 
         if(allocated(i%nznttk)) deallocate(i%nznttk)
         if(allocated(i%nznbtk)) deallocate(i%nznbtk)
         if(allocated(i%xlngtk)) deallocate(i%xlngtk)
 
         if(allocated(i%nztkvac)) deallocate(i%nztkvac)
         if(allocated(i%xltkvac)) deallocate(i%xltkvac)
 
         if(allocated(i%irho_int)) deallocate(i%irho_int)
         if(allocated(i%ibr_int)) deallocate(i%ibr_int)
         if(allocated(i%ith_int)) deallocate(i%ith_int)
         if(allocated(i%s_int)) deallocate(i%s_int)
         if(allocated(i%s_zon)) deallocate(i%s_zon)
 
         if(allocated(i%svec)) deallocate(i%svec)
         if(allocated(i%spkg)) deallocate(i%spkg)
 
         if(allocated(i%rvec)) deallocate(i%rvec)
         if(allocated(i%zvec)) deallocate(i%zvec)
         if(allocated(i%phivec)) deallocate(i%phivec)
 
         if(allocated(i%rhovec)) deallocate(i%rhovec)
         if(allocated(i%thvec)) deallocate(i%thvec)
 
         if(allocated(i%Bxvec)) deallocate(i%Bxvec)
         if(allocated(i%Byvec)) deallocate(i%Byvec)
         if(allocated(i%Bzvec)) deallocate(i%Bzvec)
 
      endif
 
      i%n_int = 0
      i%nzcrtk = 0
      i%npts = 0
      i%nvacseg = 0
 
      i%pta = 0
      i%ptb = 0
      i%gxyz = 0
 
      i%s_to_pta = 0
      i%s_from_ptb = 0
 
      if(isize.le.0) then
         jsize=maxpts
      else
         jsize=max(10,min(maxpts,isize))
      endif
 
      if((i%maxpts .ne. jsize).or.(i%maxvac .ne. maxvac)) then
 
         if(i%guard.eq.12345) call xstraln_r8track_dealloc(i)
 
         i%maxpts = jsize
         i%maxvac = maxvac
         i%guard = 12345
 
         allocate(i%nznttk(jsize)); i%nznttk = 0
         allocate(i%nznbtk(jsize)); i%nznbtk = 0
         allocate(i%xlngtk(jsize)); i%xlngtk = 0
 
         allocate(i%nztkvac(maxvac)); i%nztkvac = 0
         allocate(i%xltkvac(0:maxvac)); i%xltkvac = 0
 
         allocate(i%irho_int(jsize)); i%irho_int = 0
         allocate(i%ibr_int(jsize)); i%ibr_int = 0
         allocate(i%ith_int(jsize)); i%ith_int = 0
         allocate(i%s_int(jsize)); i%s_int = 0
         allocate(i%s_zon(jsize)); i%s_zon = 0
 
         allocate(i%svec(jsize)); i%svec = 0
         allocate(i%spkg(jsize,4)); i%spkg = 0
 
         allocate(i%rvec(jsize)); i%rvec = 0
         allocate(i%zvec(jsize)); i%zvec = 0
         allocate(i%phivec(jsize)); i%phivec = 0
 
         allocate(i%rhovec(jsize)); i%rhovec = 0
         allocate(i%thvec(jsize)); i%thvec = 0
 
         allocate(i%Bxvec(jsize)); i%Bxvec = 0
         allocate(i%Byvec(jsize)); i%Byvec = 0
         allocate(i%Bzvec(jsize)); i%Bzvec = 0
 
      endif
    end subroutine xstraln_r8track_init
 
    subroutine xstraln_r8track_dealloc(i)
 
      type (r8_track) :: i
 
      ! cleanup: deallocate arrays inside track structure
 
      i%maxpts = 0
      i%maxvac = 0
      i%guard = 0
 
      i%n_int = 0
      i%nzcrtk = 0
      i%npts = 0
      i%nvacseg = 0
 
      i%pta = 0
      i%ptb = 0
      i%gxyz = 0
 
      deallocate(i%nznttk)
      deallocate(i%nznbtk)
      deallocate(i%xlngtk)
 
      deallocate(i%nztkvac)
      deallocate(i%xltkvac)
 
      deallocate(i%irho_int)
      deallocate(i%ibr_int)
      deallocate(i%ith_int)
      deallocate(i%s_int)
      deallocate(i%s_zon)
 
      deallocate(i%svec)
      deallocate(i%spkg)
 
      deallocate(i%rvec)
      deallocate(i%zvec)
      deallocate(i%phivec)
 
      deallocate(i%rhovec)
      deallocate(i%thvec)
 
      deallocate(i%Bxvec)
      deallocate(i%Byvec)
      deallocate(i%Bzvec)
 
    end subroutine xstraln_r8track_dealloc
 
end module xstraln_r8track
