subroutine trackr8_array(kptype,ntracks,endpts,ibflag,imcflag,itracks,ierr)
 
  use xstraln_r8track
  implicit NONE
 
  ! process an array of track problems
  ! INPUT arguments:
 
  ! ...kptype specifies form of specification of endpts:
  integer,intent(in) :: kptype  ! =0: cylindric, =1: flux, =2: cartesian
 
  integer,intent(in) :: ntracks ! number of tracks to follow
 
  ! ... start (1) and stop (2) endpoints.
  ! ... [m,rad,m] or [-,rad,rad] or [m,m,m] depending on kptype:
 
  REAL*8, intent(in) :: endpts(3,2,ntracks) ! 2 endpoints for each track
 
  logical, intent(in) :: ibflag  ! .TRUE. to compute B vector along track
  logical, intent(in) :: imcflag ! .TRUE. to compute intersections with 2d
                                 ! MC grid (irregular poloidal zones)
 
  ! OUTPUT: the track descriptions & error flag
 
  type (r8_track), intent(inout) :: itracks(ntracks)
 
  integer, intent(out) :: ierr  ! completion code, 0=OK
 
  !----------------------------------
 
  integer j,ierr0
 
  !----------------------------------
 
  ierr=0
 
  do j=1,ntracks
 
     call trackr8(kptype,endpts(1:3,1:2,j),ibflag,imcflag,itracks(j),ierr0)
     ierr=max(ierr,ierr0)
 
  enddo
 
end subroutine trackr8_array
 
!==============================================
 
subroutine trackr8(kptype,endpts,ibflag,imcflag,i,ierr)
 
  use xstraln_mod
  use xstraln_2d
  use xstraln_r8track
 
  use xplasma_calls
 
  implicit NONE
 
  ! process a single track problem
 
  ! ...kptype specifies form of specification of endpts:
  integer,intent(in) :: kptype  ! =0: cylindric, =1: flux, =2: cartesian
 
  ! ... start (1) and stop (2) endpoints.
  ! ... [m,rad,m] or [-,rad,rad] or... [m,m,m] depending on kptype:
  ! ... (R,phi,Z)    (rho,theta,zeta)  (x,y,z)
 
  REAL*8, intent(in) :: endpts(3,2)  ! 2 endpoints for each track
 
  logical, intent(in) :: ibflag  ! .TRUE. to compute B vector along track
  logical, intent(in) :: imcflag ! .TRUE. to compute intersections with 2d
                                 ! MC grid (irregular poloidal zones)
 
  ! OUTPUT: the track description & error flag
 
  type (r8_track), intent(inout) :: i
 
  integer, intent(out) :: ierr  ! completion code, 0=OK
 
  !----------------------------------
  ! local...
 
  ! quantities connected with setting up the end points:
 
  integer jregion(1)
  REAL*8 slength       ! segment length
  REAL*8 dzds,dhds,drdh,dtdh,zdet
  REAL*8 zlento_zmin,zlento_zmax,zlento_near,zlento_far
  REAL*8 radius_a,radius_b
  REAL*8 ztemp,zzmin,zzmax
 
  REAL*8 user_a(3),user_b(3) ! original user endpts (cartesion coords)
 
  !------------
  ! quantities connected with analyzing the track
 
  logical iexpand     ! set .TRUE. to expand work vectors
 
  logical flag        ! set to .TRUE. if more pts needed
 
  integer icur,icurp,icurn,inew,ij,izon,i1,i2,iflag
 
  integer isave_npts,inew_npts  ! working gridsizes
 
  REAL*8, dimension(:), allocatable :: rsave,zsave,phisave
  REAL*8, dimension(:), allocatable :: rhosave,thsave
  REAL*8, dimension(:), allocatable :: zbrvec,zbphivec
  REAL*8, dimension(:,:), allocatable :: wkvec
 
  integer, dimension(:), allocatable :: insert,iregion,kzones,iwarn
  logical, dimension(:), allocatable :: ltangent
 
  REAL*8 :: zcos,zsin,zarg1,zsmall,z0,z1,z2,z8
 
  !-------------
  ! quantities associated with finding segments within analyzed track
  !
 
  logical mcflag
 
  integer :: izr,ibr,ith, izrp,ibrp,ithp, izonb, ivacseg, int_last_enter
  REAL*8 zthn,zthnp,zhi,zdsmin,zvacseg
 
  !------------------------------------------------------------------------
  data z0/0.0/
  data z1/1.00000000000000000000/
  data z2/2.00000000000000000000/
  data z8/8.00000000000000000000/
  !------------------------------------------------------------------------
  ! initialize track data structure; set #zones crossed = 0
 
  call xstraln_r8track_init(i,i%maxpts,.FALSE.)
  iexpand=.FALSE.
 
  mcflag = imcflag
 
  zzmin = zmin
  zzmax = zmax
 
  !------------------------------------------------------------------------
  ! re-enter here if track structure is re-initialized to increase the
  ! number of grid points
 
10 continue
  if(iexpand) then
     call xstraln_r8track_expand(i)
  endif
 
  !------------------------------------------------------------------------
  ! check endpoints -- convert to cartesian for internal use
 
  ierr=0
  if(kptype.eq.0) then
 
     ! (R,phi,Z) input endpoints
 
     i%pta(3)=endpts(3,1)
     i%ptb(3)=endpts(3,2)
     call eq_xcyl(1,endpts(1:1,1),endpts(2:2,1),i%pta(1:1),i%pta(2:2))
     call eq_xcyl(1,endpts(1:1,2),endpts(2:2,2),i%ptb(1:1),i%ptb(2:2))
 
  else if(kptype.eq.1) then
 
     ! (rho,theta,zeta) input endpoints
 
     call eq_xget(1,endpts(1,1:1),endpts(2,1:1),endpts(3,1:1), &
          i%pta(1:1),i%pta(2:2),i%pta(3:3),jregion,ierr)
     if(jregion(1).gt.2) then
        write(lunmsgs,*) &
             ' ?xstraln(track): track endpoint 1, given in flux coordinates'
        write(lunmsgs,*) &
             '  is outside the plasma mapped region: ',endpts(1:3,1)
        ierr=1
     endif
 
     call eq_xget(1,endpts(1,2:2),endpts(2,2:2),endpts(3,2:2), &
          i%ptb(1:1),i%ptb(2:2),i%ptb(3:3),jregion,ierr)
     if(jregion(1).gt.2) then
        write(lunmsgs,*) &
             ' ?xstraln(track): track endpoint 2, given in flux coordinates'
        write(lunmsgs,*) &
             '  is outside the plasma mapped region: ',endpts(1:3,2)
        ierr=1
     endif
 
  else if(kptype.eq.2) then
 
     ! Cartesian input endpoints
 
     i%pta = endpts(1:3,1)
     i%ptb = endpts(1:3,2)
 
  else
     write(lunmsgs,*) ' ?xstraln(track): kptype=',kptype,' not in range (1:3)'
     ierr=1
  endif
 
  if(ierr.ne.0) RETURN
 
  !------------------------------------------------------------------------
  ! adjust endpoints if not in the [Rmin,Rmax]x[Zmin,Zmax] annular region
 
  !------------------------------------------------------------------------
  ! #1.  confine to region between two planes Z=Zmin and Z=Zmax
 
  slength = sqrt( &
       (i%ptb(1)-i%pta(1))**2 + &
       (i%ptb(2)-i%pta(2))**2 + &
       (i%ptb(3)-i%pta(3))**2 )
 
  if(slength.eq.z0) then
     ! zero length track means no intersection
     RETURN
  endif
 
  ! save original end points: segment might get shortened...
 
  user_a = i%pta
  user_b = i%ptb
 
  i%gxyz = (i%ptb-i%pta)/slength    ! direction of segment
 
  dzds = i%gxyz(3)                  ! dead zero means a horizontal track
 
  ! in this block we adjust A and B necessary, to force that the entire
  ! track segment lie between Zmin and Zmax
  if(dzds.eq.z0) then
     if(i%pta(3).gt.zmax) GO TO 999    ! horizontal track, Z > Zmax
     if(i%pta(3).lt.zmin) GO TO 999    ! horizontal track, Z < Zmin
  else
     zlento_zmin=(zmin-i%pta(3))/dzds    ! distance, i%pta to Z=Zmin
     zlento_zmax=(zmax-i%pta(3))/dzds    ! distance, i%pta to Z=Zmax
     zlento_near=min(zlento_zmin,zlento_zmax)
     zlento_far =max(zlento_zmin,zlento_zmax)
     if((zlento_far.le.z0).or.(zlento_near.ge.slength)) then
        GO TO 999                      ! segment never between Zmin and Zmax
     endif
 
     if(zlento_near.gt.z0) then
        i%pta = i%pta + zlento_near*i%gxyz   ! move A up to Z=Zmin or Z=Zmax
        ztemp=i%pta(3)
        i%pta(3)=max(zzmin,min(zzmax,ztemp))   ! be sure
        slength = slength-zlento_near
        zlento_far=zlento_far-zlento_near
        zlento_near=z0
     endif
 
     if(slength.gt.zlento_far) then
        i%ptb = i%pta + zlento_far*i%gxyz    ! move B back to Z=Zmin or Z=Zmax
        ztemp=i%ptb(3)
        i%ptb(3)=max(zzmin,min(zzmax,ztemp))   ! be sure
        slength=zlento_far
     endif
 
     ! now:  zlento_near <= 0 < slength <= zlento_far
  endif
 
  !------------------------------------------------------------------------
  ! #2.  deal with concentric cylinders R=Rmin and R=Rmax
  ! first compute starting radius and dhds = dL(horizontal)/dL
 
  radius_a=sqrt(i%pta(1)**2 + i%pta(2)**2)
  dhds = sqrt(i%gxyz(1)**2 + i%gxyz(2)**2)  ! dhds = 0 means vertical track
 
  if(dhds.eq.z0) then
 
     if((radius_a.lt.Rmin).or.(radius_a.gt.Rmax)) then
        GO TO 999         ! vertical track not btw Rmin and Rmax
     else
        continue       ! vertical track btw Rmin and Rmax and we are done
     endif
 
  else
 
     ! track is not dead vertical
 
     if(radius_a.eq.z0) then
        ! we are starting on axis; distance to cylinders is simply Rmin, Rmax
        if(slength.lt.Rmin) then
           GO TO 999    ! track never makes it to Rmin
        else
           if(slength.gt.Rmax) then
              i%ptb = i%pta + Rmax*i%gxyz/dhds  ! trim track to Rmax
              slength = Rmax/dhds
           endif
           i%pta = i%pta + Rmin*i%gxyz/dhds     ! move start to Rmin
           slength = sqrt( &                    ! and check length
                (i%ptb(1)-i%pta(1))**2 + &
                (i%ptb(2)-i%pta(2))**2 + &
                (i%ptb(3)-i%pta(3))**2 )
           if(slength.le.z0) GO TO 999
           radius_b=sqrt(i%ptb(1)**2 + i%ptb(2)**2)
           if(radius_b.le.Rmin) GO TO 999
        endif
 
     else
 
        ! track starts away from cylinder axis
        ! compute (at pt A) d(tangential distance)/d(horizontal distance)
        !                   d(radial distance)/d(horizontal distance)
 
        dtdh=(-i%pta(2)*i%gxyz(1)+i%pta(1)*i%gxyz(2))/(radius_a*dhds)
        drdh=( i%pta(1)*i%gxyz(1)+i%pta(2)*i%gxyz(2))/(radius_a*dhds)
 
        ! cylinder intercepts are solutions to quadratic equation
        ! deal with outer cylinder first
 
        zdet = Rmax**2 - radius_a**2*dtdh**2
 
        if(zdet.le.z0) then
           GO TO 999        ! track is outside or at best tangent to R=Rmax
        else
           ! line containing segment does come inside R=Rmax
           zlento_near=(-radius_a*drdh-sqrt(zdet))/dhds  ! near intercept dist.
           zlento_far =(-radius_a*drdh+sqrt(zdet))/dhds  ! far intercept dist.
           if((zlento_far.le.z0).or.(zlento_near.ge.slength)) then
              GO TO 999                      ! segment never inside Rmax
           endif
 
           if(zlento_near.gt.z0) then
              i%pta = i%pta + zlento_near*i%gxyz   ! move A up to R=Rmax
              slength = slength-zlento_near
              zlento_far=zlento_far-zlento_near
              zlento_near=z0
              ! need to update pt A radius and direction info
              radius_a=sqrt(i%pta(1)**2 + i%pta(2)**2)
              dtdh=(-i%pta(2)*i%gxyz(1)+i%pta(1)*i%gxyz(2))/(radius_a*dhds)
              drdh=( i%pta(1)*i%gxyz(1)+i%pta(2)*i%gxyz(2))/(radius_a*dhds)
           endif
 
           if(slength.gt.zlento_far) then
              i%ptb = i%pta + zlento_far*i%gxyz    ! move B back to R=Rmax
              slength=zlento_far
           endif
        endif
 
        !------------------------------------------------------------------
        ! deal with inner cylinder.  The only time this has an effect is
        ! when a track starting outside R=Rmin hits the cylinder.  If the
        ! track is starting inside the Rmin cylinder or is moving away
        ! from it, the "near" intercept is sure to be negative.
 
        zdet = Rmin**2 - radius_a**2*dtdh**2
 
        if(zdet.le.z0) then
           continue        ! track never hits inner cylinder:  OK.
        else
           zlento_near=(-radius_a*drdh-sqrt(zdet))/dhds  ! near intercept dist.
           if(zlento_near.gt.z0) then
              ! Rmin intercept is ahead of start point
              if(slength.gt.zlento_near) then
                 i%ptb = i%pta + zlento_near*i%gxyz    ! move B back to R=Rmin
                 ztemp=i%ptb(3)
                 i%ptb(3)=max(zzmin,min(zzmax,ztemp))    ! be sure
                 slength=zlento_near
              endif
           else
              continue     ! near intercept not in way of segment.
           endif
        endif
 
     endif  ! (radius_a =? 0)
  endif  ! (track is/is not vertical)
 
  !-----------------------------------------------------------------------
  ! OK-- track segment, entirely contained in [Rmin,Rmax]x[Zmin,Zmax],
  ! has been defined; proceed to find track intersection lengths.
  !
  !  method: build a 1d grid in length s along the track segment.
  !  ...an "adaptive grid" is built to assure accuracy where the
  !  target surfaces are located, and extra accuracy at edge and
  !  in tangency regions...
 
  iexpand = .FALSE.    ! this will be set if more grid points are needed.
 
  call loc_allocate
 
  !
  ! initial evaluation:  max(11,nrho_target) points
  !
 
  i%npts=max(11,nrho_target)
 
  icur=1
 
  ! on first call, move in slightly from endpoints
  ! this prevents a segment that is at the boundary of
  ! the xplasma mapped region from being declared outside
 
  zsmall=1.0d-6
  wkvec(1,icur)=zsmall*slength
  wkvec(i%npts,icur)=slength-wkvec(1,icur)
 
  call sgen(2,i%npts-1,z0,slength)
  call seval(ierr)
  if(ierr.ne.0) then
     write(lunmsgs,*) ' ?xstraln(track) failure on 1st seval call.'
     call loc_free
     go to 999
  endif
 
  ! restore endpts
 
  wkvec(1,icur)=z0
  wkvec(i%npts,icur)=slength
 
  if(ierr.ne.0) return
  call stats(flag)
 
  !----------------------------
  ! add to grid until all conditions satisfied:
  !   a) sufficient resolution near tangencies
  !   b) sufficient resolution near edges
  !   c) at least one grid point between each pair of target surfaces
  !   d) at least 4 grid points
 
  do
     if(.not.flag) exit            ! all done
 
     icurp=icur
     icurn=4-icur                  ! 3,1,3,1,...
     icur=2
 
     allocate(rsave(i%npts),zsave(i%npts),phisave(i%npts))
     allocate(rhosave(i%npts),thsave(i%npts))
     allocate(insert(i%npts))
 
     rsave=i%rvec(1:i%npts)
     zsave=i%zvec(1:i%npts)
     phisave=i%phivec(1:i%npts)
     rhosave=i%rhovec(1:i%npts)
     thsave=i%thvec(1:i%npts)
     insert=0
 
     isave_npts=i%npts
     i%npts=0
 
     inew_npts=1
     wkvec(1,icurn)=wkvec(1,icurp)
 
     !--------------------------------
     ! next iteration on grid...
     !
     ! in what follows:
     !   isave_npts = number of old grid points
     !   i%npts = number of new grid points
     !   inew_npts = i%npts + isave_npts = number after merger of grids
 
     ! wkvec(:,icurp) = {set of old points}
     ! wkvec(:,icur) = {set of new points}
     ! wkvec(:,icurn) = {set of old points} U {set of new points}
 
     do ij=2,isave_npts
        izon=ij-1
 
        i1=0
        i2=i%npts
        if(kzones(izon).gt.1) then
           i1=i%npts+1
           i2=i%npts+1+kzones(izon)
        else if(ltangent(izon)) then
           i1=i%npts+1
           i2=i%npts+3
        endif
 
        if(i2.gt.i%maxpts) then
           iexpand=.TRUE.
           exit
        endif
 
        if(i1.gt.0) then
           call sgen(i1,i2,wkvec(ij-1,icurp),wkvec(ij,icurp))
           inew=i2-i1+1
        else
           inew=0
        endif
        i%npts=i2
 
        if(inew_npts+inew+1.gt.i%maxpts) then
           iexpand=.TRUE.
           exit
        endif
 
        insert(ij)=inew
 
        if(inew.gt.0) then
           wkvec(inew_npts+1:inew_npts+inew,icurn)=wkvec(i1:i2,icur)
        endif
        inew_npts=inew_npts+inew+1
        wkvec(inew_npts,icurn)=wkvec(ij,icurp)
 
     enddo
 
     if(.not.iexpand) then
 
        ! new grid seems OK; evaluate track data
 
        call seval(ierr)
 
        iflag=0
        if(ierr.ne.0) then
           write(lunmsgs,*)  &
                ' ?xstraln(track): failure on nth seval call.'
           iflag=1  ! this should not happen.
        else
           ! merge data associated with new grid
 
           call trk_merge(rsave,isave_npts,i%rvec,inew_npts,i%npts,insert)
           call trk_merge(zsave,isave_npts,i%zvec,inew_npts,i%npts,insert)
           call trk_merge(phisave,isave_npts,i%phivec,inew_npts,i%npts,insert)
           call trk_merge(rhosave,isave_npts,i%rhovec,inew_npts,i%npts,insert)
           call trk_merge(thsave,isave_npts,i%thvec,inew_npts,i%npts,insert)
 
           ! remove possible branch cut jumps
 
           call dbranch(i%phivec,inew_npts)
           call dbranch(i%thvec,inew_npts)
 
        endif
     endif
 
     ! cleanup
     deallocate(rsave,zsave,phisave,rhosave,thsave,insert)
 
     if(iexpand) then
        call loc_free
        go to 10    ! restart with more buffer space
     endif
 
     if(iflag.eq.1) then
        call loc_free
        go to 999   ! give up after error.
     endif
 
     ! check merged grid & re-enter loop
 
     icur = icurn
     i%npts = inew_npts
     if(no_fRZ) then
        call rho_patch(i%npts,wkvec(1:i%npts,icur),i%rhovec)
     endif
 
     call stats(flag)
  enddo
 
  !----------------------------------------------------------
  ! have grid... now set up piecewise linear interpolating functions
 
  i%svec(1:i%npts)=wkvec(1:i%npts,icur)
 
  !----------------------------------------------------------
  ! setup pspline fast lookup on adaptive, unevenly spaced grid
  !
  ! enabled Nov2002 RGA for eflux
  !
  zsmall=1.0d-5
  call genxpkg(i%npts,i%svec,i%spkg,0,1,0,zsmall,-3,iflag)
  if(iflag.ne.0) then
     write(lunmsgs,*) ' ?xstraln(track): genxpkg failure.'
     write(lunmsgs,*) '  i%npts = ',i%npts
     write(lunmsgs,*) '  i%svec = ',i%svec(1:i%npts)
     call loc_free
     go to 999
  endif
 
  if(ibflag) then
     ! get B field; convert (BR,Bphi,BZ) -> (Bx,By,Bz)
 
     allocate(zbrvec(i%npts),zbphivec(i%npts),iwarn(i%npts))
     zarg1 = 1
     call xs_brz_r8(i%npts,zarg1,i%phivec,i%rvec,i%zvec,zarg1, &
          zbrvec,i%bzvec,zbphivec,iwarn)
 
     do ij=1,i%npts
        zcos=cos(i%phivec(ij))
        zsin=sin(i%phivec(ij))
        i%bxvec(ij)=zbrvec(ij)*zcos-zbphivec(ij)*zsin
        i%byvec(ij)=zbrvec(ij)*zsin+zbphivec(ij)*zcos
     enddo
 
     deallocate(zbrvec,zbphivec,iwarn)
  else
     ! just clear the B vectors
 
     do ij=1,i%npts
        i%bxvec(ij)=z0
        i%byvec(ij)=z0
        i%bzvec(ij)=z0
     enddo
  endif
 
  call loc_free
 
  !-----------------------------------------------------------------------
  ! OK:  have rho(l), theta(l) ... treat as piecewise linear functions and
  ! compute radial zone and 2d beam zone intersections...
 
  zvacseg=z0
  ivacseg=0
  i%xltkvac(0)=z0
 
  i%n_int=1
  i%irho_int(1)=0
  i%ibr_int(1)=0
  i%ith_int(1)=0
  i%s_int(1)=z0
 
  i%nzcrtk=0
  int_last_enter = -1
 
  izr=nrho_target
  ibr=inbr
  ith=0
  call izones(1,i%rhovec(1),i%thvec(1),izr,ibr,ith,zthn)
  zdsmin=dsmin/z8
 
  do ij=2,i%npts
     izrp=izr
     ibrp=ibr
     ithp=ith
     zthnp=zthn
 
     call izones(ij,i%rhovec(ij),i%thvec(ij),izr,ibr,ith,zthn)
 
     if((izr.eq.izrp).and.(ibr.eq.ibrp).and.(ith.eq.ithp)) then
     else
        zhi=z1/(i%svec(ij)-i%svec(ij-1))  ! 1/h -- inverse step size
        call asegment(ij-1,ij,izrp,izr,ibrp,ibr,ithp,ith,zthnp,zthn)
     endif
 
  enddo
 
  !-----------------------------------------------------------------------
  ! exit.  record distances by which original user endpoints were moved
  !        go get internal endpoints used in calculation.
 
999 continue
 
  i%s_to_pta = sqrt( &
       (user_a(1)-i%pta(1))**2 + &
       (user_a(2)-i%pta(2))**2 + &
       (user_a(3)-i%pta(3))**2 )
 
  i%s_from_ptb = sqrt( &
       (i%ptb(1)-user_b(1))**2 + &
       (i%ptb(2)-user_b(2))**2 + &
       (i%ptb(3)-user_b(3))**2 )
 
 
  contains
 
    !=====================================
 
    subroutine loc_allocate
 
      allocate(wkvec(i%maxpts,3))
      allocate(iregion(i%maxpts),kzones(i%maxpts))
      allocate(ltangent(i%maxpts))
 
    end subroutine loc_allocate
 
    !=====================================
 
    subroutine loc_free
 
      deallocate(wkvec)
      deallocate(iregion,kzones,ltangent)
 
    end subroutine loc_free
 
    !=====================================
 
    subroutine sgen(i1,i2,s1m,s2p)
 
      !  fill in wkvec(i1:i2,icur) with points evenly spaced
      !  btw s1m and s2p but not including s1m,s2p.
 
      integer, intent(in) :: i1,i2
      REAL*8, intent(in) :: s1m,s2p
 
      !----------------------
 
      integer ij
 
      do ij=i1,i2
         wkvec(ij,icur)=((i2-ij+1)*s1m + (ij-i1+1)*s2p)/(i2-i1+2)
      enddo
 
    end subroutine sgen
 
    !=====================================
 
    subroutine seval(ierr)
      integer, intent(out) :: ierr   ! 0=OK
 
      ! compute R,Z,phi,rho,theta  vs.  s (in wkvec)
 
      integer ij
      REAL*8 zx,zy
 
      REAL*8, parameter :: ztol = 1.0d-2  ! polar map OK
 
      integer iwarn(i%npts)
 
      ierr=0
 
      ! evaluate mappings
 
      do ij=1,i%npts
         zx=i%pta(1)+i%gxyz(1)*wkvec(ij,icur)
         zy=i%pta(2)+i%gxyz(2)*wkvec(ij,icur)
         i%zvec(ij)=i%pta(3)+i%gxyz(3)*wkvec(ij,icur)
         i%rvec(ij)=sqrt(zx**2+zy**2)
         if(i%rvec(ij).eq.z0) then
            i%phivec(ij)=z0
         else
            i%phivec(ij)=atan2(zy,zx)
         endif
      enddo
 
      if(no_fRZ) then
         call eq_inv(i%npts,i%rvec,i%zvec,i%phivec,ztol, &
              i%rhovec,i%thvec,iwarn,ierr)
      else
         zarg1=1
         call xs_fxget_r8(i%npts,zarg1,i%phivec,i%rvec,i%zvec, &
              i%rhovec,i%thvec,iwarn)
      endif
 
      ! check/remove branch cuts (theta)
      call dbranch(i%thvec,i%npts)
 
      ! check/remove branch cuts (phi)
      call dbranch(i%phivec,i%npts)
 
    end subroutine seval
 
    subroutine rho_patch(inum,xl,xrho)
 
      !  extrapolate rho(L) in regions beyond plasma bdy
      !  only needed when xstraln is in no_fRZ mode...
 
      integer, intent(in) :: inum
      REAL*8 :: xl(inum)           ! length vector L
      REAL*8 :: xrho(inum)         ! rho(L)
 
      !-----------------------------
      integer in1,in2,iout1,iout2,iexit,i
      REAL*8, parameter :: zrho_outmin = 1.0001d0
      REAL*8 :: zs, zsa, ztrans
      integer :: in1a,in2a
      integer :: inleft, inright
      integer :: iregion(inum)   ! =1 -- inside; >1 -- outside
      real*8, parameter :: zxtest = 0.999999999999d0
      !-----------------------------
 
      do i=1,inum
         if(xrho(i).gt.zxtest) then
            iregion(i)=2
         else
            iregion(i)=1
         endif
      enddo
 
      iexit=0
 
      if(iregion(1).gt.1) then
         !  patch leading edge
         in1=1
         do
            in1=in1+1
            if(in1.gt.inum) then
               iexit=1  ! nothing inside is found
               xrho(1:inum)=max(zrho_outmin,xrho(1:inum))
               exit
            else if(iregion(in1).eq.1) then
               iout1=1
               iout2=in1-1
               in2=in1+1
               exit
            endif
         enddo
 
         if(iexit.eq.1) return
 
         if(iregion(in2).gt.1) then
            zs = 0
         else
            zs = (xrho(in1)-xrho(in2))/(xl(in1)-xl(in2))
         endif
         do i=iout1,iout2
            xrho(i)=max(zrho_outmin,xrho(in1)+zs*(xl(i)-xl(in1)))
         enddo
         inleft = in1
      else
         inleft = 1
      endif
 
      if(iregion(inum).gt.1) then
         in1=inum
         do
            in1=in1-1
            if(iregion(in1).eq.1) then
               iout1=in1+1
               iout2=inum
               in2=max(1,in1-1)
               exit
            endif
         enddo
 
         if(in1.eq.in2) then
            zs = 0
         else if(iregion(in2).gt.1) then
            zs = 0
         else
            zs = (xrho(in1)-xrho(in2))/(xl(in1)-xl(in2))
         endif
         do i=iout1,iout2
            xrho(i)=max(zrho_outmin,xrho(in1)+zs*(xl(i)-xl(in1)))
         enddo
         inright=in1
      else
         inright=inum
      endif
 
      ! OK look for interior regions of track that are outside the plasma
 
      iout1=inleft
      do
         iout1=iout1+1
         if(iout1.ge.inright) exit  ! done
         if(iregion(iout1).eq.1) cycle  ! keep looking...
 
         ! found an interior region of track that is outside plasma
         ! get extrapolation info from left side
 
         in1=iout1-1
         in2=max(inleft,in1-1)
         if(in1.eq.in2) then
            zs = 0
         else if(iregion(in2).gt.1) then
            zs = 0
         else
            zs = (xrho(in1)-xrho(in2))/(xl(in1)-xl(in2))
         endif
 
         ! look for right hand side of region
 
         iout2=iout1
         do
            iout2=iout2+1
            if(iregion(iout2).eq.1) exit
         enddo
         in1a=iout2
         iout2=iout2-1
         in2a=min(inright,in1a+1)
 
         ! extrapolation info from right side
 
         if(in1a.eq.in2a) then
            zsa=0
         else if(iregion(in2a).gt.1) then
            zsa=0
         else
            zsa = (xrho(in1a)-xrho(in2a))/(xl(in1a)-xl(in2a))
         endif
 
         !  calculate the patch; blend extrapolations from left and right
         do i=iout1,iout2
            ztrans=(xl(i)-xl(in1))/(xl(in1a)-xl(in1))
            xrho(i) = max(zrho_outmin, &
                 (1-ztrans)*(xrho(in1)+zs*(xl(i)-xl(in1))) + &
                 ztrans*(xrho(in1a)+zsa*(xl(i)-xl(in1a))) )
         enddo
 
         !  resume search
 
         iout1=iout2+1
      enddo
 
    end subroutine rho_patch
 
    !=====================================
 
    subroutine dbranch(zthvec,inpts)
 
      integer, intent(in) :: inpts
      REAL*8, intent(inout) :: zthvec(inpts)
 
      ! remove branch cut jumps from an angle sequence
 
      REAL*8 :: zdang,zadjust
 
      zadjust=z0
      ! check/remove branch cuts (theta)
      do ij=1+1,inpts
         zthvec(ij)=zthvec(ij)+zadjust
         zdang=zthvec(ij)-zthvec(ij-1)
         if(zdang.gt.xpi) then
            zadjust=zadjust-twopi
            zthvec(ij)=zthvec(ij)-twopi
         else if(zdang.lt.-xpi) then
            zadjust=zadjust+twopi
            zthvec(ij)=zthvec(ij)+twopi
         endif
      enddo
 
    end subroutine dbranch
 
    !=====================================
 
    subroutine stats(flag)
      logical, intent(out) :: flag
 
      ! compute statistics; flag=.TRUE. if more work to do
 
      integer iminm(i%npts),iedgz,ij,isave,idir,ibdy,ikzth,ifirst,ilast,ict
      logical iminsee,imaxsee
      REAL*8 zrhotes1,zrhotes2,zrhotes3,zrhotes4,zdthmax
      REAL*8 zrhomin,zrhomax,zrhostep,zthstep,zdrho(4),zds
 
      REAL*8, dimension(:), allocatable :: zdist,zdum
 
      !  *** limit max =1 re-entries into the plasma
      !  *** this will need to change if we ever go to non-axisymmetric geometry!
 
      integer, parameter :: max_renter = 1
      integer ienter
 
      !-------------------------
 
      flag = .FALSE.
 
      iminm=0
 
      !-------------------------
      ! first see if track can be shortened
      ! by cutting track short if a vacuum segment strikes a limiter.
 
      ifirst=0
      do ij=1,i%npts
         if(i%rhovec(ij).le.rho_target(nrho_target)) then
            ifirst=ij
            exit
         endif
      enddo
 
      if(ifirst.gt.0) then
 
         ! count entry events; truncate track if too many...
 
         ienter=-1
         do ij=1,i%npts-1
            if (ienter<0) then
               if (i%rhovec(ij).gt.rho_target(nrho_target)) cycle  ! have not entered plasma
               ienter=0
            end if
            if((i%rhovec(ij).gt.rho_target(nrho_target)).and. &
                 (i%rhovec(ij+1).le.rho_target(nrho_target))) then
               ienter = ienter+1
               if(ienter.gt.max_renter) then
                  write(lunmsgs,*) ' %trackr(xstraln): too many re-entries, track truncated.'
                  !  if this happens a lot:  axisymmetry?
                  !  if it happens rarely: rho(R,Z) is imperfect, can seem to happen with axisymmetry
                  i%npts = ij
                  exit
               endif
            endif
         enddo
 
         do
            ilast=0
            do ij=ifirst+1,i%npts
               if(i%rhovec(ij).gt.rho_target(nrho_target)) then
                  ilast=ij-1
                  exit
               endif
            enddo
            if(ilast.eq.0) exit  ! no closing vacuum segment
 
            ifirst=i%npts+1
            do ij=ilast+1,i%npts
               if(i%rhovec(ij).le.rho_target(nrho_target)) then
                  ifirst=ij
                  exit
               endif
            enddo
 
            ict=ifirst-ilast-1
            allocate(zdist(ict),zdum(ict))
 
            zdist=z0
            call eq_rzdist(ict, &
                 i%rvec(ilast+1:ifirst-1),i%zvec(ilast+1:ifirst-1), &
                 i%phivec(ilast+1:ifirst-1),zdist,0,zdum,zdum,zdum,ierr)
            do ij = ilast+1,ifirst-1
               if(zdist(ij-ilast).gt.z0) then
                  ! limiter struck in vacuum segment: truncate track
                  i%npts=ij
                  exit
               endif
            enddo
 
            deallocate(zdist,zdum)
 
            if(ifirst.ge.i%npts) exit
 
            ! cycle to check more vacuum segments.
         enddo
      endif
 
      !-------------------------
      ! find zones adjacent to local minima or maxima & inside plasma
      ! mark zones next to min/max marked points whose steplengths are
      ! greater than a reference step length (dsref/2)
 
      ! endpoints count if they are between the target surfaces
 
      if((rho_target(1).le.i%rhovec(1)).and. &
           (i%rhovec(1).le.rho_target(nrho_target))) then
         iminm(1)=1
      endif
      if((rho_target(1).le.i%rhovec(i%npts)).and. &
           (i%rhovec(i%npts).le.rho_target(nrho_target))) then
         iminm(i%npts)=1
      endif
 
      do ij=2,i%npts-1
         if((i%rhovec(ij).ge.i%rhovec(ij-1)).and. &
              (i%rhovec(ij).ge.i%rhovec(ij+1))) then
            ! local maximum
            iminm(ij)=1
            if(i%rhovec(ij-1).gt.i%rhovec(ij+1)) iminm(ij-1)=1
            if(i%rhovec(ij+1).gt.i%rhovec(ij-1)) iminm(ij+1)=1
         else if((i%rhovec(ij).le.i%rhovec(ij-1)).and. &
              (i%rhovec(ij).le.i%rhovec(ij+1))) then
            ! local minimum
            iminm(ij)=1
            if(i%rhovec(ij-1).lt.i%rhovec(ij+1)) iminm(ij-1)=1
            if(i%rhovec(ij+1).lt.i%rhovec(ij-1)) iminm(ij+1)=1
         endif
      enddo
 
      do ij=1,i%npts-1
         ltangent(ij)=.FALSE.
         if(iminm(ij)+iminm(ij+1).gt.0) then
            if(wkvec(ij+1,icur)-wkvec(ij,icur).gt.dsref/2) then
               ltangent(ij)=.TRUE.
               flag=.TRUE.
            endif
         endif
      enddo
 
      !-------------------------
      ! check rho-spacing of zones.  aim for the following resolution:
      !   del = spacing of target surfaces
      !
      !   del/8 at boundaries
      !   del/2 inside boundaries
      !   del*2 outside boundaries
      !
      ! also check theta-spacing of zones.
      !   also aim for delta(theta) <~ 0.1 radians
      !
      ! but steps should not get less than dsmin
 
      zdrho(1)=(rho_target(2)-rho_target(1))
      zdrho(2)=rho_target(nrho_target)-rho_target(nrho_target-1)
      zdrho(3)=max(zdrho(1),zdrho(2))
 
      zdthmax=xpi/10
 
      !-------------------------
      ! find number of target surfaces crossed in each ds zone
 
      isave=1
      idir=1
 
      do ij=1,i%npts-1
         zrhomin=min(i%rhovec(ij),i%rhovec(ij+1))
         zrhomax=max(i%rhovec(ij),i%rhovec(ij+1))
         zrhostep=zrhomax-zrhomin
         zthstep=abs(i%thvec(ij+1)-i%thvec(ij))
         zds=wkvec(ij+1,icur)-wkvec(ij,icur)
         ikzth = 1 + zthstep/zdthmax
 
         if(zrhostep.eq.z0) then
            kzones(ij)=1
            cycle
         endif
 
         if((zrhomax.lt.rho_target(1)-zdrho(1)/z2).or. &
              (zrhomin.gt.rho_target(nrho_target)+zdrho(2)/z2)) then
 
            ! for zones well beyond the boundaries
 
            kzones(ij)=max(ikzth, &
                 int( z1 + zrhostep/(z2*zdrho(3))))
 
         else if(zrhomax.gt.rho_target(nrho_target)) then
 
            ! straddling outer boundary
 
            kzones(ij)=max(ikzth, int( z1 + zds/dsref), &
                 int( z1 + z8*zrhostep/(zdrho(2))))
 
         else if(zrhomin.lt.rho_target(1)) then
 
            ! straddling inner boundary
 
            kzones(ij)=max(ikzth, int( z1 + zds/dsref), &
                 int( z1 + z8*zrhostep/(zdrho(1))))
 
         else
 
            iminsee=.false.
            imaxsee=.false.
            zdrho(4)=rho_target(nrho_target)-rho_target(1) ! safe upper bd.
            isave=isave-idir
            do
               isave=isave+idir
               if(isave.ge.nrho_target) then
                  imaxsee=.true.
                  idir=-1
               else if(isave.le.1) then
                  iminsee=.true.
                  idir=1
               else if(rho_target(isave).lt.zrhomin) then
                  iminsee=.true.
                  idir=1
               else if(rho_target(isave).gt.zrhomax) then
                  imaxsee=.true.
                  idir=-1
               else
                  ztemp = (rho_target(isave+1)-rho_target(isave-1))/2
                  zdrho(4)=min(zdrho(4),ztemp)
               endif
               if(iminsee.and.imaxsee) exit
            enddo
 
            ! straddling 0 or more interior zones
 
            kzones(ij)=max(ikzth, int( z1 + zds/dsref), &
                 int( z1 + z2*zrhostep/(zdrho(4))))
         endif
         kzones(ij)=min(4,int(1+zds/dsmin),kzones(ij))
         if(kzones(ij).gt.1) flag=.TRUE.
      enddo
      if(.not.flag) then
         if(i%npts.lt.4) then
            kzones(1:i%npts-1)=2
            flag=.TRUE.
         endif
      endif
 
    end subroutine stats
 
    !=====================================
    ! merge vectors...
 
    subroutine trk_merge(oldvec,n_old,newvec,n_new,n_add,insert)
 
      integer, intent(in) :: n_old   ! old vector size
      integer, intent(in) :: n_new   ! new vector size after merger
      integer, intent(in) :: n_add   ! number of new points being added
 
      REAL*8, intent(in) :: oldvec(n_old)  ! old vector
      REAL*8, intent(inout) :: newvec(n_new)  ! new vector
 
      integer, intent(in) :: insert(n_old)  ! where to insert new elements
 
      ! on input:  newvec(:) contains just the new points
      ! on output:  newvec(:) contains new & old points merged
 
      !---------------------------------
      REAL*8 :: buf(n_add)  ! temp storage for new pts
      integer ijk,iprevn,iprevb,inew
      !---------------------------------
 
      buf = newvec(1:n_add)
 
      newvec(1)=oldvec(1)
      iprevn=1
      iprevb=0
 
      do ijk=2,n_old
         inew=insert(ijk)
         if(inew.gt.0) then
            newvec(iprevn+1:iprevn+inew) = buf(iprevb+1:iprevb+inew)
            iprevb=iprevb+inew
         endif
         iprevn=iprevn+inew+1
         newvec(iprevn)=oldvec(ijk)
      enddo
 
    end subroutine trk_merge
 
    !=====================================
 
    subroutine izones(ia, zrho,zth, izr, ibr, ith, zthu)
 
      integer, intent(in) :: ia  ! vector address of pt.
      REAL*8, intent(in) :: zrho   ! radial coordinate
      REAL*8, intent(in) :: zth    ! poloidal angle coordinate
 
      integer, intent(inout) :: izr   ! radial zone
      integer, intent(inout) :: ibr   ! beam zone row
      integer, intent(inout) :: ith  ! theta zone w/in beam zone row
 
      !----------------
      REAL*8 zthu
      !----------------
 
      zthu = zth
      do while(zthu < -xpi)
         zthu=zthu+twopi
      enddo
      do while(zthu > xpi)
         zthu=zthu-twopi
      enddo
 
      if(zrho.lt.rho_target(1)) then
         izr=0  ! out of range low (?)
         ibr=0
         ith=0
      else if(zrho.ge.rho_target(nrho_target)) then
         izr=nrho_target   ! out of range high
         ibr=inbr+1
         ith=0
      else
         ! find rho zone
         izr=max(1,min((nrho_target-1),izr))
         do
            if(rho_target(izr).gt.zrho) then
               izr=izr-1
            else if(rho_target(izr+1).lt.zrho) then
               izr=izr+1
            else
               exit
            endif
         enddo
         ! beam zone row
         ibr=izonrb(izr)
         ! find theta zone
         call izonth(i%rvec(ia),i%zvec(ia),ith,ibr,zthu)
      endif
    end subroutine izones
 
    !=====================================
 
    subroutine izonth(zzr,zzz,ith,ibr,zthu)
 
      REAL*8, intent(in) :: zzr,zzz    ! (R,Z) position (needed for ibr=1)
      integer, intent(inout) :: ith  ! beam th zone index (may be updated)
      integer, intent(in) :: ibr     ! beam zone row index
      REAL*8, intent(in) :: zthu       ! angle, normalized to range [-pi,pi]
 
      !---------------
      REAL*8 zrpt,zzpt,zcr1,zcr2,zdot
      integer ii
      !---------------
 
      if(.not.mcflag) then
         ith=0                   ! ignore theta zones
      else if(ibr.gt.inbr) then
         ith=0
      else if(ibr.lt.1) then
         ith=0
      else
         ith=1+(inthbr(ibr)-1)*(zthu+xpi)/twopi
         ith=max(1,min((inthbr(ibr)-1),ith))
         do
            if(thbr(ith,ibr).gt.zthu) then
               if(ith.eq.1) exit
               ith=ith-1
            else if(thbr(ith+1,ibr).lt.zthu) then
               if(ith.eq.inthbr(ibr)-1) exit
               ith=ith+1
            else
               exit
            endif
         enddo
         if(ibr.eq.1) then
            !  innermost zone: use (R,Z) vector information to determine
            !  zone membership; do not rely on nominal theta from fast fit
            zrpt=zzr-Raxis
            zzpt=zzz-Zaxis
            ith=ith-1
            do ii=1,inthbr(1)-1
               ith=ith+1
               if(ith.eq.inthbr(1)) ith=1
               zdot=zrpt*(rdel1(ith)+rdel1(ith+1)) + &
                    zzpt*(zdel1(ith)+zdel1(ith+1))
               zcr1=-zzpt*rdel1(ith)+zrpt*zdel1(ith)
               zcr2=-zzpt*rdel1(ith+1)+zrpt*zdel1(ith+1)
               if(((zcr1*zcr2).le.z0).and.(zdot.ge.z0)) then
                  exit
               endif
            enddo
         endif
      endif
 
    end subroutine izonth
 
    !=====================================
 
    subroutine asegment(i1,i2,iz1,iz2,ib1,ib2,ith1,ith2,zth1,zth2)
 
      integer, intent(in) :: i1,i2   ! segment endpoints: note i2=i1+1
      integer, intent(in) :: iz1,iz2 ! radial zone indices at i1,i2
      integer, intent(in) :: ib1,ib2 ! beam zone row indices at i1,i2
      integer, intent(in) :: ith1,ith2  ! beam theta zone indices at i1,i2
      REAL*8, intent(in) :: zth1,zth2  ! normalized theta (in range [pi,2pi])
                                     ! at i1,i2
 
      !--------------------
      integer iz,izp,isurf,inumz,idir,ij,ijp,ith,ik
      integer, dimension(:), allocatable :: iza,iba,ith1a,ith2a
      REAL*8, dimension(:), allocatable :: zth1a,zth2a,zl1a,zl2a
      REAL*8 zrho1,zrho2,zdrho,zdrhoi,zrat,zthu,zdthtot,zzr,zzz
      !--------------------
 
      zzr=z0
      zzz=z0
      if(iz1.eq.iz2) then
 
         ! (implies ib1.eq.ib2 also)
         ! find theta break & mark in output arrays:
 
         call asegth(i1,i2,i%svec(i1),i%svec(i2),iz1,ib1, &
              ith1,ith2,zth1,zth2)
 
      else
 
         ! there are radial zone boundaries inside this step
         ! find them
 
         zrho1=i%rhovec(i1)
         zrho2=i%rhovec(i2)
         zdrho=zrho2-zrho1  ! can be negative
         zdrhoi=z1/zdrho    ! can be negative
         zdthtot=i%thvec(i2)-i%thvec(i1)   ! can be negative
         if(iz1.lt.iz2) then
            idir=1
            inumz=iz2-iz1+1   ! no. of radial zones this step
         else
            idir=-1
            inumz=iz1-iz2+1   ! no. of radial zones this step
         endif
 
         !  iza: zone index    iba: beam zone row index
         !  ith1a,zth1a: beam th index & value, entering this zone
         !  ith2a,zth2a: beam th index & value, exiting this zone
         !  zl1a,zl2a: length along track on entry and exit of this zone
 
         allocate(iza(inumz),iba(inumz),ith1a(inumz),ith2a(inumz))
         allocate(zth1a(inumz),zth2a(inumz),zl1a(inumz),zl2a(inumz))
         ij=0
         do iz=iz1,iz2,idir
            ijp=ij
            ij=ij+1
            iza(ij)=iz
            iba(ij)=izonrb(iz)
 
            if(ij.eq.1) then
               ! start of step
               ith1a(ij)=ith1
               ith=ith1
               zth1a(ij)=zth1
               zl1a(ij)=i%svec(i1)
            else
               ! intermediate zone boundary
               izp=iz-idir
               isurf=max(iz,izp)
 
               zrat=(rho_target(isurf)-zrho1)*zdrhoi
               zl1a(ij)=i%svec(i1)+zrat*(i%svec(i2)-i%svec(i1))
               zl2a(ijp)=zl1a(ij)
 
               zthu=i%thvec(i1)+(zl1a(ij)-i%svec(i1))*zdthtot*zhi
               do while(zthu < -xpi)
                  zthu=zthu+twopi
               enddo
               do while(zthu > xpi)
                  zthu=zthu-twopi
               enddo
               zth1a(ij)=zthu
               zth2a(ijp)=zthu
 
               if(min(iba(ijp),iba(ij)).eq.1) then
                  zzr=i%rvec(i1)+ &
                       (zl1a(ij)-i%svec(i1))*zhi*(i%rvec(i2)-i%rvec(i1))
                  zzz=i%zvec(i1)+ &
                       (zl1a(ij)-i%svec(i1))*zhi*(i%zvec(i2)-i%zvec(i1))
               endif
               call izonth(zzr,zzz,ith,iba(ijp),zthu)
               ith2a(ijp)=ith
               call izonth(zzr,zzz,ith,iba(ij),zthu)
               ith1a(ij)=ith
 
               if(ij.eq.inumz) then
                  ! end of step
                  ith2a(ij)=ith2
                  zth2a(ij)=zth2
                  zl2a(ij)=i%svec(i2)
               endif
            endif
         enddo
 
         ! process all boundaries...
 
         do ij=1,inumz
            if(ith1a(ij).ne.ith2a(ij)) then
               call asegth(i1,i2,zl1a(ij),zl2a(ij),iza(ij),iba(ij), &
                    ith1a(ij),ith2a(ij),zth1a(ij),zth2a(ij))
            endif
            if(ij.lt.inumz) then
               if(zl2a(ij).gt.i%s_int(i%n_int)+zdsmin) then
 
                  ! radial boundary intersection
 
                  i%n_int=i%n_int+1
                  ik=i%n_int
                  i%s_int(ik)=zl2a(ij)
                  i%irho_int(ik)=max(iza(ij),iza(ij+1))
                  if(iba(ij).ne.iba(ij+1)) then
                     i%ibr_int(ik)=max(iba(ij),iba(ij+1))
                  else
                     i%ibr_int(ik)=0
                  endif
                  i%ith_int(ik)=0
 
                  if(iza(ij).lt.nrho_target) then
                     i%nzcrtk = i%nzcrtk+1
                     ik=i%nzcrtk
                     i%xlngtk(ik)=i%s_int(i%n_int)-i%s_int(i%n_int-1)
                     i%s_zon(ik) =(i%s_int(i%n_int)+i%s_int(i%n_int-1))/2.
                     i%nznttk(ik)=iza(ij)
                     izonb=ibzns(ith2a(ij),iba(ij))
                     i%nznbtk(ik)=izonb
                  endif
 
               endif
               ! check for vacuum segment
 
               if((iza(ij).eq.nzones).and.(iza(ij+1).eq.nzones+1)) then
                  ! leaving the target
                  if (i%n_int /= int_last_enter) then
                     zvacseg=zl2a(ij)  ! normal, a boundary or segment between enter and exit
                     ivacseg=i%nzcrtk
                  else
                     if (i%n_int==0 .or. i%nvacseg==0) then
                        zvacseg = 0.  ! previous entry was from the initial vacuum start
                     else
                        zvacseg = i%s_int(i%n_int) - i%xltkvac(i%nvacseg)  ! previous zvacseg
                        i%nvacseg = i%nvacseg-1
                     end if
                     i%n_int=max(0, i%n_int-1)   ! remove boundary of last entry
                     int_last_enter = -1         ! just in case
                  end if
               else if((iza(ij).eq.nzones+1).and.(iza(ij+1).eq.nzones)) then
                  ! entering the target
                  int_last_enter = i%n_int
                  if(zvacseg.gt.z0) then
                     i%nvacseg=i%nvacseg+1
                     i%nztkvac(i%nvacseg)=ivacseg
                     i%xltkvac(i%nvacseg)=zl2a(ij)-zvacseg
                  else
                     i%xltkvac(0)=zl2a(ij)
                  endif
               endif
 
            endif
         enddo
      endif
 
    end subroutine asegment
 
    !=====================================
 
    subroutine asegth(i1,i2,zl1,zl2,izr,ibr,ith1,ith2,zth1,zth2)
 
      !  find theta boundaries within a radial zone
 
      integer i1,i2    ! index to i%svec segment boundaries
      REAL*8 zl1,zl2     ! substep within i%svec(i1) to i%svec(i2)
      integer izr      ! radial index
      integer ibr      ! beam radial zone row index
      integer ith1,ith2  ! theta zone indices
      REAL*8 zth1,zth2   ! normalized theta (in range [-pi,pi])
 
      !----------------
      integer ii,ik,inthzons,inth2,idir,ishift,inumth,ithloc
      REAL*8, dimension(:), allocatable :: zthbds
      integer, dimension(:), allocatable :: ithbds
      REAL*8 zoffset,zdthnorm,zlen
      REAL*8 zzr1,zzz1,zzr2,zzz2
      !----------------
 
      inthzons=inthbr(ibr)-1
      inth2=inthzons/2
 
      if(ibr.gt.1) then
         ! which direction in index space...
         if(ith2.gt.ith1) then
            idir=1
            ishift=1
         else
            idir=-1
            ishift=0
         endif
         inumth=idir*(ith2-ith1)
         if((idir*(zth2-zth1)).gt.xpi) then
            idir=-idir
            ishift=1-ishift
            inumth=inthzons-inumth
         endif
 
      ! th boundaries hit, & bdy at end of segment (zth2)
      ! detect branchcut if present
 
         allocate(zthbds(inumth+1),ithbds(inumth))
         zoffset=z0
         ithloc=ith1-idir
         do ii=1,inumth+1
            ithloc=ithloc+idir
            if(ithloc.gt.inthzons) then
               ithloc=ithloc-inthzons
               zoffset=zoffset+twopi
            else if(ithloc.lt.1) then
               ithloc=ithloc+inthzons
               zoffset=zoffset-twopi
            endif
            if(ii.le.inumth) then
               ithbds(ii)=ithloc+ishift
               zthbds(ii)=thbr(ithloc+ishift,ibr)+zoffset
            else
               zthbds(inumth+1)=zth2+zoffset
            endif
         enddo
 
         ithloc=ith1-idir
         do ii=1,inumth
            ithloc=ithloc+idir
            if(ithloc.gt.inthzons) then
               ithloc=ithloc-inthzons
            else if(ithloc.lt.1) then
               ithloc=ithloc+inthzons
            endif
            zdthnorm = (zthbds(ii)-zth1)/(zthbds(inumth+1)-zth1)
            zlen=zl1+zdthnorm*(zl2-zl1)
            if(zlen.gt.(i%s_int(i%n_int)+zdsmin)) then
 
               i%n_int=i%n_int+1
               ik=i%n_int
               i%s_int(ik)=zlen
               i%irho_int(ik)=0
               i%ibr_int(ik)=0
               i%ith_int(ik)=ithbds(ii)
 
               i%nzcrtk = i%nzcrtk+1
               ik=i%nzcrtk
               i%xlngtk(ik)=i%s_int(i%n_int)-i%s_int(i%n_int-1)
               i%s_zon(ik) =(i%s_int(i%n_int)+i%s_int(i%n_int-1))/2.
               i%nznttk(ik)=izr
               izonb=ibzns(ithloc,ibr)
               i%nznbtk(ik)=izonb
 
            endif
         enddo
 
         deallocate(zthbds,ithbds)
 
      else
 
         zzr1=i%rvec(i1)+(zl1-i%svec(i1))*zhi*(i%rvec(i2)-i%rvec(i1))-Raxis
         zzz1=i%zvec(i1)+(zl1-i%svec(i1))*zhi*(i%zvec(i2)-i%zvec(i1))-Zaxis
         zzr2=i%rvec(i1)+(zl2-i%svec(i1))*zhi*(i%rvec(i2)-i%rvec(i1))-Raxis
         zzz2=i%zvec(i1)+(zl2-i%svec(i1))*zhi*(i%zvec(i2)-i%zvec(i1))-Zaxis
 
         call asegth0(izr,ith1,zl1,zzr1,zzz1,ith2,zl2,zzr2,zzz2)
      endif
 
    end subroutine asegth
 
    !=========================================
 
    subroutine asegth0(izr,ith1,zl1,zzr1,zzz1,ith2,zl2,zzr2,zzz2)
 
      ! special theta segmenting method for innermost zone row
      ! (due to polar singularity, standard method is inaccurate)
 
      integer, intent(in) :: izr         ! radial zone row
      integer, intent(in) :: ith1,ith2   ! theta zones (beam zone row 1)
      REAL*8, intent(in) :: zl1,zl2        ! lengths to segment end pts
      REAL*8, intent(in) :: zzr1,zzr2      ! R dist. from mag. axis
      REAL*8, intent(in) :: zzz1,zzz2      ! Z dist. from mag. axis
 
      !----------------
      REAL*8 zlloc,zdist,zdrn,zdzn,zdet,zs,zl,zlen
      integer idir,ithb,ii,inthzons,ik,ithcur
      !----------------
 
      inthzons=inthbr(1)-1
      zdist=sqrt((zzr2-zzr1)**2+(zzz2-zzz1)**2)
      zdrn=(zzr2-zzr1)/zdist
      zdzn=(zzz2-zzz1)/zdist
 
      idir=0
      do
         if(idir.eq.0) then
            idir=1
            ithb=ith1+1-idir
            ithcur=ith1-idir
         else if(idir.eq.1) then
            idir=-1
            ithb=ith1-idir
            ithcur=ith1-idir
         else
            exit  ! nothing found: assume a near miss
         endif
 
         do ii=1,inthzons
 
            !  find intersection between lines; if the intersection
            !  lies in both segments, accept it; there could be
            !  multiple intercepts...
 
            ithb=ithb+idir
            if(ithb.lt.1) ithb=ithb+inthzons
            if(ithb.gt.inthzons) ithb=ithb-inthzons
 
            ithcur=ithcur+idir
            if(ithcur.lt.1) ithcur=ithcur+inthzons
            if(ithcur.gt.inthzons) ithcur=ithcur-inthzons
 
            zdet=zdzn*rdel1(ithb)-zdrn*zdel1(ithb)
            if(zdet.eq.z0) exit
            zs=(zdzn*zzr1-zdrn*zzz1)/zdet
            zl=(zdel1(ithb)*zzr1-rdel1(ithb)*zzz1)/zdet
            if(zs.lt.z0) exit
            if(zs.gt.z1) exit
            if(zl.lt.z0) exit
            if(zl.gt.zdist) exit
 
            !  solution acceptable
 
            zlen = zl1+(zl/zdist)*(zl2-zl1)
            if (zlen.gt.(i%s_int(i%n_int)+zdsmin)) then
 
               i%n_int=i%n_int+1
               ik=i%n_int
               i%s_int(ik)=zlen
               i%irho_int(ik)=0
               i%ibr_int(ik)=0
               i%ith_int(ik)=ithb
 
               i%nzcrtk = i%nzcrtk+1
               ik=i%nzcrtk
               i%xlngtk(ik)=i%s_int(i%n_int)-i%s_int(i%n_int-1)
               i%s_zon(ik) =(i%s_int(i%n_int)+i%s_int(i%n_int-1))/2.
               i%nznttk(ik)=izr
               izonb=ibzns(ithcur,1)
               i%nznbtk(ik)=izonb
 
            endif
          enddo
 
         if(ii.gt.1) exit
 
      enddo
 
    end subroutine asegth0
 
end subroutine trackr8
