 
module map_larmor_mod
!mg comment 03/13/12:
!some confusion with coordinate system here:
!(x,y,z) is a (r,phi,z) 
  !
  ! this is a series of subroutines for finding the point which satisfies
  !   x_point = x_gc + (bvec_point x v_point)/(b_point * omega_point)
  !
  ! the subroutine map_larmor_gc needs to be called first to fill guiding
  ! center and gyrophase angle dependent data before the other subroutines
  ! can be called.  this will be done automatically for map_larmor and
  ! map_larmor_plot if imap_larmor_badgc=.true. so this flag should be
  ! set when guiding center data is modified (nbigc, ...).
  ! the gyrophase dependent data (zm_gc) is always recomputed so there is
  ! no need to set imap_larmor_badgc when cospho,sinpho is changed.
  !
  ! the convention for vpll employed here is that the sign is relative
  ! to the poloidal current direction.  also, the poloidal flux psi
  ! returned by xplasma is always positive so it is adjusted based on
  ! the poloidal current direction (nsnccw*nsjdotb).
  !
  ! a possibly incomplete list of trcom variables which are used in these
  ! subroutines is:
  !   xion  th  rmjion  yion  cospho  sinpho  xbjaco  vpll  vperp  vion
  !   gflr_op(10) 1-> 0.0 for newton iteration, otherwise secant
  !               2-> convergence criteria for x,v equation (suggest 1.e-4)
  !               3-> convergence criteria for a,b (suggest 1.e-5)
  !               4-> 0.0 to use map_larmor_loss, otherwise use limlos (orbit)
  !               5-> maximum step size (3.0)
  !               6-> maximum number of iterations in map_larmor_newton (try 9)
  !               7-> /=0. -> accept points which did not converge but
  !                           are bracketed in the v equation  (try 1.0)
  !               8-> acceptance criteria on x equation for closeness (try 0.2)
  !               9-> acceptance criteria on v equation for closeness (try 0.1)
  !              10-> nonzero -> try the other newton or secant method
  !                   if the original one fails to get close
  !              11-> threshold above which to try rescaling the search region.
  !                   this is a threshold relative to (rlar/r)*(xksid**2/1.-xksid**2)
  !                   which appears in a first order analysis of the v equation (try 5.)
  !
  !        gflr_ll -> maximum value allowed for a,b parameters (default 30.0)
  !        gflr_rl -> particles further then this (or gflr_ll) in larmor radii
  !                   from the limiter are considered in the plasma
  !     gflr_xv(2) -> criteria for judging validity of x equation and v equation
  !                   at closest limiter position.  particle is good if
  !                   abs(x equation)<gflr_xv(1) and (v equation)>gflr_xv(2)
  !
 
  use iso_c_binding, only: fp => c_double
  use xplasma_calls
  use physconst_mod, only: zel,vc,zmtocm
 
  implicit none
  private :: fp
  logical, private :: idiagnostic=.false.            ! .true. for diagnostic messages not normally
                                                     ! desired
  logical, private :: idiag_test=.false.             ! more diagnostic
 
  logical, parameter, private :: ibffget=.false.     ! .true. to call eqxyz_bffget directly in
                                                     ! subroutine map_larmor_bffget
 
  integer, parameter, private :: ivecd = 15  ! maximum number of points per xplasma call -- check
                                             ! map_larmor_jacobian, map_larmor, map_larmor_search
                                             ! before decreasing
  integer, parameter, private :: ifcnd = 2   ! number of functions to evaluate per xplasma call
 
  real(fp), parameter, private :: zconv_psi = 1.e8_fp  ! convert webers to maxwells
  real(fp), parameter, private :: zz_e_c    =  zel/(vc*zmtocm)  ! e/c
  real(fp), parameter, private :: ztol      = 1.e-6_fp  ! guess at a tolerance for xplasma
  real(fp), parameter, private :: zmicro_step = 0.4_fp  ! step reduction for each micro iteration in search
 
  !
  ! the following variables are available for use in any subroutine
  ! so they may change across subroutine calls
  !
  integer, private, save :: ifcns(ifcnd)           ! list of xplasma functions
  integer, private, save :: ierr                   ! error code from xplasma
 
  real(fp), private, save :: zbvec(3, ivecd)           ! b vector returned by xplasma call
  real(fp), private, save :: zfvals(ivecd, ifcnd)      ! xplasma function values
  real(fp), private, save :: zz_xx(ivecd),zz_yy(ivecd),zz_zz(ivecd) ! holds coordinates in meters for xplasma
 
  !
  ! data which is setup by map_larmor_gc for use in the other subroutines,
  ! private to minimize name pollution -- can be made public
  !
  logical,save :: imap_larmor_badgc=.true.   ! .true. when guiding center data is bad
 
  integer,save,private :: jregion_gc         ! the guiding center region code from xplasma
                                             ! this is either 1 or 4
 
  real(fp),save,private :: zx_gc(3)              ! x,y,z of guiding center -- on y-plane
  real(fp),save,private :: zbvec_gc(3)           ! guiding center bx,by,bz
  real(fp),save,private :: zb_gc                 ! total b field at guiding center
  real(fp),save,private :: zpsi_gc               ! poloidal flux in webers at guiding center
  real(fp),save,private :: zphi_gc               ! electric potential in volts at guiding center
  real(fp),save         :: zxi_gc(3)             ! unit vector orthogonal to magnetic surface
  real(fp),save         :: zperp_gc(3)           ! unit vector on magnetic surface orthogonal to b
  real(fp),save         :: zpar_gc(3)            ! unit vector in direction of b
  real(fp),save,private :: zm_gc(3)              ! direction of ideal larmor radius displacement
  real(fp),save,private :: z_p_phi_gc            ! toroidal angular momentum at guiding center
  real(fp),save,private :: z_e_gc                ! kinetic energy at guiding center
  real(fp),save,private :: zrlar_gc              ! ideal larmor radius at guiding center
  real(fp),save,private :: zscale_gc             ! scale factor for dimensionless search in distance
  real(fp),save,private :: zmin_scale_gc         ! minimum distance in scale radii for a stable
                                             ! finite difference
 
  integer,save :: ic_loss(6) = (/0,0,0,0,0,0/) ! count number of calls to sections
                                               ! in map_larmor_loss
  integer,save :: ic_map(4)  = (/0,0,0,0/)     ! count number of map_larmor calls
                                               ! for each value of ioops+1
  real(fp),save,private :: ip_dir = 1              ! direction of toroidal current
  real(fp),save         :: zvr_ml, zvy_ml, zvtor_ml ! velocity in rhat, zhat, phihat directions to compare
                                                ! with depsub
 
contains
 
  !
  ! ============================== map_larmor_gc ======================================
  !
  subroutine map_larmor_gc
    !
    ! this subroutine grabs some guiding center data needed by the other
    ! subroutines -- imap_larmor_badgc is .true. if there was an error
    !
    use nbi_com
    !
    logical :: im                    ! .true. for nonorthogonality messages
    logical :: ixmomry               ! .true. when jacobian is based on xmomry called within
                                     ! this subroutine
 
    integer :: jregion(ivecd)        ! xplasma code for region of point
 
    real(fp) :: zz, zz2                  ! temporary
    real(fp) :: zzjaco(2,2)              ! jacobian
 
    imap_larmor_badgc=.true.         ! data is being updated
    ip_dir = nsnccw*nsjdotb
    !
    ! put gc at phi = 0.
    !
    zx_gc = (/rmjion, 0._fp, yion/)
 
    ifcns(1) = id_psi
    ifcns(2) = id_phi   ! electrostatic potential
 
    !
    ! it would be faster calling eq_bffget but we want to stay consistent
    ! with the call used for the real particle point since the displacement
    ! from the guiding center is small
    !
    zz_xx(1) = zx_gc(1)/100._fp  ! the xplasma cm to m got me again -- argh!!
    zz_yy(1) = zx_gc(2)/100._fp
    zz_zz(1) = zx_gc(3)/100._fp
    if(ifcns(2).gt.0) then
       call map_larmor_bffget(1, zz_xx, zz_yy, zz_zz, ztol, zbvec, ifcnd, ifcns, &
            ivecd, zfvals, jregion, ierr)
    else
       call map_larmor_bffget(1, zz_xx, zz_yy, zz_zz, ztol, zbvec, 1, ifcns(1:1), &
            ivecd, zfvals(:,1), jregion, ierr)
       zfvals(:,2)=0.0_fp
    endif
 
    if (ierr /= 0 .or. jregion(1)>3) then
       if (jregion(1)>3) then
          write(nout,'(a)') '? guiding center is outside r,z grid'
       end if
       write(nout,'(a)') '? map_larmor_gc: error in call to eq_bffget'
       call map_larmor_info
       return
    end if
 
    jregion_gc = jregion(1)
 
    zpsi_gc = ip_dir*zconv_psi*zfvals(1,1)  ! in maxwells
    zphi_gc = zfvals(1,2)                   ! keep in units of volts
 
    ! --- form parallel unit vector ---
    ! note: btoroidal/b = zpar_gc(2)
    !
    zbvec_gc = (/ zbvec(1,1), zbvec(2,1), zbvec(3,1) /)  ! b field at gc in teslas (bvec=x,y,z)
    zb_gc    = sqrt(sum(zbvec_gc*zbvec_gc))              ! total b field at gc in teslas
    zpar_gc  = zbvec_gc/zb_gc                            ! unit vector in parallel direction
    zrlar_gc = rmassb*vperp/(xzbeami*zz_e_c*1.e4_fp*zb_gc)  ! include tesla->gauss factor
    zscale_gc     = zrlar_gc
    zmin_scale_gc = raxis*ztol*10._fp/max(zscale_gc,1.e-10_fp*raxis) ! guess
 
    !
    ! --- form xi unit vector ---
    ! vector orthogonal to drvec/dtheta
    !
    ixmomry = .false.   ! first try the original jacobian for xi vector
    zzjaco  = xbjaco    ! copy jacobian
10  continue            ! return to this point if xi and b not orthogonal and try again
 
    zz = sqrt(zzjaco(1,2)**2 + zzjaco(2,2)**2)
    if (zz /= 0._fp) then
       zxi_gc = (/zzjaco(2,2)/zz, 0._fp, -zzjaco(1,2)/zz/) ! (dz/dtheta, 0., -dr/dtheta) / ||....||
    else
       zxi_gc = (/1._fp, 0._fp, 0._fp/)  ! on axis -- any direction will do
    end if
 
    zz = sum(zxi_gc*zpar_gc)    ! sin(angle error from 90 degrees) between xi and b vectors
    if (abs(zz)>.02_fp) then       ! test for orthogonality of xi and b (.02 -> 1.15 degrees)
       im = (xion<1.0_fp) .or. (abs(zz)>.2_fp)  ! reduce number of messages
       if (im) then
          write(nout,'(a)') '! map_larmor_gc: b vector and xi vector not orthogonal'
          write(nout,'(a,f8.3,a)') '   angle = ', 57.296_fp*asin(zz), ' degrees'
       end if
       if (.not. ixmomry) then
          ixmomry = .true.
          if (im) write(nout,'(a)') ' trying to recover with xmomry...'
          call xmomry(xion, th, 1, zz, zz2, zzjaco)   ! don't change original stuff
          goto 10
       end if
       if (im) then
          write(nout,'(a)') '  enforcing orthogonality and continuing'
          call map_larmor_info
       end if
       zxi_gc = zxi_gc - zz*zpar_gc               ! orthogonal to b
       zxi_gc = zxi_gc / sqrt(sum(zxi_gc*zxi_gc)) ! make a unit vector
    end if
 
    !
    ! --- form perp unit vector ---
    ! perp = xi x par
    !
    zperp_gc(1) = zxi_gc(2)*zpar_gc(3) - zxi_gc(3)*zpar_gc(2)
    zperp_gc(2) = zxi_gc(3)*zpar_gc(1) - zxi_gc(1)*zpar_gc(3)
    zperp_gc(3) = zxi_gc(1)*zpar_gc(2) - zxi_gc(2)*zpar_gc(1)
 
    ! --- form m vector, unit vector in ideal guiding center direction ---
    zm_gc = cospho*zxi_gc + sinpho*zperp_gc
 
    !
    ! --- canonical angular momentum  m*r*vtor-q*psi/c  ---
    ! by convention in depsub  vtor = ip_dir * vpll * |btor|/|b|  for the guiding center
    !
    ! do not save canonical angular momentum explicitly because psi can have a really
    ! different magnitude then the velocity term
    !
    z_p_phi_gc = ip_dir*rmassb*rmjion*vpll*abs(zpar_gc(2))
 
    ! --- kinetic energy in ev ---
    z_e_gc = vion*vion*v2toev
 
    imap_larmor_badgc=.false.   ! data has been updated
 
  end subroutine map_larmor_gc
 
  !
  ! ============================== map_larmor ======================================
  !
  subroutine map_larmor(ioops, iregion, zr, zy, zvpar, zvperp, alpha_v, ierr_write)
    !
    ! find the point on the gyro which satisfies  v = omega * (delta_x x bpar/b)
    ! using a newton or secant iteration method
    !
    ! map_larmor_gc must be called before this subroutine
    !
    use nbi_com
 
    logical            :: i_newton     ! .true. for newton iterations, .false. for secant iterations
    logical            :: i_method     ! newton/secant method actually used in iteration
    logical            :: i_close      ! .true. when iterations got close but did not converge
 
    real(fp)               :: zxv_conv     ! convergence criteria for x,v equations
    real(fp)               :: zab_conv     ! end iterations when change in position relative to a
                                       ! larmor radius is smaller than this
 
    logical, intent(in) :: ierr_write  ! .true. to write out error messages (not called from limlos)
                                       ! .false. to suppress error messages ( called from limlos )
 
    integer, intent(out) :: ioops      ! 0 -> no errors, larmor radius from gyro was found (or did not need it)
                                       ! 1 -> failure, used classical displacement from guiding center
                                       ! 2 -> major failure, used guiding center with no displacement
                                       ! 3 -> iterations failed but used a v equation bracketed point
    integer, intent(out) :: iregion    ! region code for point returned
                                       ! since eqxyz_bffget is no longer being called, the region code
                                       ! is either 1 or 4 with 2,3 equivalent to 1
 
    integer :: jr(ivecd)               ! region code
    integer :: jmax                    ! actual number of iterations+1
    integer :: jlimit                  ! maximum number of iterations allowed
    integer :: i_try                   ! number of retries of iteration
    integer :: ifail(ivecd)            ! nonzero for equation failure
 
    real(fp), intent(out) :: zr, zy        ! position of the particle on the gyro
    real(fp), intent(out) :: zvpar, zvperp ! parallel and perpendicular components of the velocity
                                       ! vector on the gyro relative to the b field direction
                                       ! at the guiding center (with an nsnccw factor thrown in)
 
    real(fp) :: za(ivecd), zb(ivecd)       ! position parallel and orthogonal to b in larmor radii
    real(fp) :: zxeq(ivecd), zveq(ivecd)   ! values of dimensionless x and v equation for each iteration
    real(fp) :: zxp(3,ivecd)               ! x,y,z position of point for each iteration
    real(fp) :: zvec(3,ivecd)              ! vx,vy,vz velocity of the particle for each iteration
    real(fp) :: zfact                      ! the factor before the linear term in the dominant quadratic
                                       ! equation which governs the v equation
    real(fp) :: zomega                     ! guiding center gyrofrequency
 
    logical :: idiag,ibad
    integer :: j1max,i,j
    real(fp) :: zzr=.01_fp, zzk=1.1_fp, zz2=2._fp  ! for testing
    real(fp) :: z1, z2, z3
    integer :: jregion_wrk(ivecd)
    real(fp) :: zcos, zsin     ! b field direction --
    real(fp) :: v_g,vg_perp !unit vectors for coordinate system (e_g,e_ll,e_perp) == (zxi_gc,zpar_gc,zperp_gc)
    !  where,
    !  e_g is the unit vector normal to the flux surface
    !  e_ll is a unit vector defining the direction of the magnetic field, b/|b|
    !  e_perp = e_g x e_ll
    !  v = (v_g , v_ll, vg_perp) 

    real(fp),intent(out) :: alpha_v  !  angle  between normal to the flux surface  and vperp = (v_g,vg_perp)
                   !   0 degrees-> straight up;
                   ! +/-  45 degrees-> on a cone in an upward direction
                   ! +/-  90 degrees-> in a horizontal direction...
                   ! +/- 135 degrees-> on a cone in a downward direction
                   ! +/- 180 degrees-> straight down

    jlimit    = min(ivecd,int(gflr_op(6)+.01_fp))  ! maximum number of iterations+1 -- keep jmax<=ivecd
    i_newton  = gflr_op(1) == 0._fp  ! first choice
    zxv_conv  = gflr_op(2)
    zab_conv  = gflr_op(3)
 
    ! --- check ---
    if (imap_larmor_badgc) then
       call map_larmor_gc                       ! update guiding center data
       if (imap_larmor_badgc) call nbi_exit     ! exit if fails
    else
       zm_gc = cospho*zxi_gc + sinpho*zperp_gc  ! repeat in case cospho,sinpho changed
    end if
 
    ! --- begin ---
    i_method = i_newton   ! the desired method to use
    i_try    = 0          ! count retries
5   continue              ! restart point for other method retry
    ioops = 0
    jmax = jlimit         ! maximum number of iterations+1, changed in call
 
    ! -- look for actual point --
    ! start at classical larmor position
    za(1) = 1.0_fp  ! position based on b field at guiding center
    zb(1) = 0.0_fp
 
    call map_larmor_newton(jmax, i_method, zab_conv, zxv_conv, za, zb, jr, &
         zxeq, zveq, zxp, zvec, i_close)
 
    zscale_gc = zrlar_gc                                       ! may have been changed below
    zmin_scale_gc = raxis*ztol*10._fp/max(zscale_gc,1.e-10_fp*raxis) ! guess
 
    ! -- errors in iteration --
    j1max = jmax                    ! save for debugging
    if (jmax == 0) then             ! first point was bad
       if (ierr_write) write(nout,*) '! map_larmor: first iteration failed'
    else if (jr(jmax)>3) then       ! oops, iterations failed
       if (ierr_write) write(nout,*) '! map_larmor: iterations took gyro position outside r,z grid'
    else if (abs(zxeq(jmax))>zxv_conv .or. abs(zveq(jmax))>zxv_conv) then
       if (ierr_write) then
          write(nout,*) '! map_larmor: iterations did not converge'
          write(nout,'(a,g14.6,a,g14.6,a,g14.6)') ' xeq = ', zxeq(jmax), '  veq = ', zveq(jmax),&
               '  xksid = ', xksid
       end if
    else
       if (i_try/=0 .and. ierr_write) write(nout,*) '! map_larmor: other method worked!'
       goto 10            ! success
    end if
 
    ! -- almost last chance --
    if (i_close .and. gflr_op(7)/=0._fp.and. jr(jmax)<4) then
       if (ierr_write) then
          write(nout,*) '! map_larmor: using the velocity bracketed point'
          if (i_try /= 0) write(nout,*) '! map_larmor: other method worked well enough -- using velocity bracket'
       end if
       ioops = 3
       goto 10            ! somewhat success
    end if
 
    ! -- really last chance, except ...--
    if (i_try==0 .and. jmax==jlimit .and. gflr_op(10)/=0._fp) then
       !
       ! first method had no early exit but failed to converge so now
       ! attempt to try the other method newton/secant, sometimes works
       !
       i_try = i_try+1
       i_method = .not. i_newton
       if (ierr_write) write(nout,*) '! map_larmor: first method failed (newton or secant) attempting other one'
       goto 5
    end if
 
    ioops = 1             ! darn
    if (ierr_write) then
       if (i_try /= 0) write(nout,*) '! map_larmor: alternative method failed'
       call map_larmor_info  ! informative output for errors
    end if
 
    ! --- get zxp,zvec the old fashioned way, from the guiding center ---
    jmax  = 1      ! the first point holds the data
    za(1) = 1.0_fp  ! one larmor radius
    zb(1) = 0.0_fp
    call map_larmor_points(1, za, zb, jr, zxeq, zveq, zxp, zvec, ifail=ifail)
    if (jr(1)>3 .or. ifail(1)/=0) then
       if (ierr_write) then
          write(nout,*) '! map_larmor: classical larmor displacement is outside r,z grid'
          write(nout,*) '              using guiding center'
       end if
       ioops = 2           ! this is not good
       jr(1) = jregion_gc  ! force the region code to the guiding center
       zxp(:,1)  = zx_gc
       
       !mg 03/13/12      zvec(:,1) = nsjdotb*vpll*zpar_gc + vperp*zperp_gc
       zvec(:,1) = nsjdotb*vpll*zpar_gc + vperp*zm_gc !mg 03/13/12: zxi_gc component has been forgotten
    end if
 
    !
    ! large xksid recovery, use larmor radius value of v equation
    ! to judge the direction of the new scale.  the scale assumes
    ! that the 1/r length scale dominates over the bphi/b length
    ! scale which may not be correct but the derivative of bphi/b
    ! is not as readily available.  see map_larmor_plot if you
    ! want to add this.
    !
    zfact = (zrlar_gc/rmjion)*(xksid**2)/max(1.e-5_fp,(1-xksid**2))
    if (i_try==0 .and. zfact>gflr_op(11) .and. ioops==1) then
       !
       ! look for new dominant scale factor
       !
       i_try = i_try+1                               ! get only one shot
       zomega = (xzbeami*zz_e_c*1.e4_fp*zb_gc)/rmassb   ! gyro frequency
       if (zveq(1)>0._fp) then
          zscale_gc = (rmjion/2._fp)*(1-xksid**2)/(xksid**2)   ! rlar/(2*zfact), probably very small
       else
          zscale_gc = 2._fp*(vpll/zomega)**2/rmjion            ! 2*rlar*zfact, relatively big
       end if
       zmin_scale_gc = raxis*ztol*10._fp/max(zscale_gc,1.e-10_fp*raxis) ! ugly coding but works
       if (ierr_write) write(nout,'(a,g14.5)') &
            ' ! map_larmor: attempt new search with rescaling to scale=',zscale_gc
       goto 5     ! restart with new scale
    end if
 
    !
    ! -- done grab data --
    !
10  zr = sqrt(zxp(1,jmax)*zxp(1,jmax)+zxp(2,jmax)*zxp(2,jmax))
    !zxp(1:3,jmax) = (x,y,z)
    zy = zxp(3,jmax)
    ! --- get poloidal aiming angle
    zz_xx(1)=zxp(1,jmax)/100._fp
    zz_yy(1)=zxp(2,jmax)/100._fp
    zz_zz(1)=zxp(3,jmax)/100._fp
    !
    call map_larmor_bffget(1,zz_xx , zz_yy, zz_zz, &
         & ztol, zbvec, ifcnd, ifcns, ivecd, zfvals, jregion_wrk, ierr)
    ! now we have bx,by,bz = zbvec(1:3)
    ! find b field direction --
    if (zr == 0._fp) then
       zcos = 1._fp
       zsin = 0._fp
    else
       zcos = zxp(1,jmax)/zr
       zsin = zxp(2,jmax)/zr
    end if
    
    !  unit vectors for coordinate system (e_g,e_ll,e_perp) == (zxi_gc,zpar_gc,zperp_gc)
    !  where,
    !  e_g is the unit vector normal to the flux surface
    !  e_ll is a unit vector defining the direction of the magnetic field, b/|b|
    !  e_perp = e_g x e_ll
    !  v = (v_g, v_ll, vg_perp) 

    v_g=zzdot(zvec,zxi_gc) !velocity component in orthogonal to magnetic surface
    vg_perp=zzdot(zvec,zperp_gc) !velocity component on magnetic surface orthogonal to b
    alpha_v=atan2(vg_perp,v_g)

!this would give an angle s.t. alpha=0 always means perpendicular to the flux surface and outward; 
!alpha = +/- pi/2 are tangential to the flux surface in opposite directions one with an inward flr 
!displacement to gc and the other with outward flr displacement, and alpha= +/- pi perpendicular and inward.  
!the angle then would have physical meaning w.r.t. real space gradients in gc distribution and/or proximity 
!to loss cone regions with gc f=0. 
 
    zvpar  = nsjdotb*sum(zpar_gc*zvec(:,jmax))               ! parallel component with sign
                                                             ! relative to ip
    zvperp = sqrt(max(0._fp,sum(zvec(:,jmax)*zvec(:,jmax))-zvpar*zvpar))
 
    iregion = jr(jmax)
 
    if (ierr_write) ic_map(ioops+1) = ic_map(ioops+1)+1     ! collect statistics outside of limlos
 
    ! -- for comparison with depsub --
    zvr_ml   = (zxp(1,jmax)*zvec(1,jmax)  + zxp(2,jmax)*zvec(2,jmax))/zr
    zvy_ml   = zvec(3,jmax)
    zvtor_ml = (-zxp(2,jmax)*zvec(1,jmax) + zxp(1,jmax)*zvec(2,jmax))/zr
 
    ! -- debugging --
    if (.false. .and. ioops/=0 .and. ioops/=3 .and. rlar>zzr .and. abs(xksid)<zzk .and. ierr_write) then
       idiag = idiagnostic
       idiagnostic = .true.
       j1max = max(j1max,jmax)
       write(*,*) 'ierr_write=', ierr_write
       write(*,*) 'ioops=', ioops
       write(*,*) 'a=', za(j1max), '  b=', zb(j1max)
       write(*,*) 'enter 0->small, limit->large'
       read(*,*) zr
       if (zr .le. 0._fp) then
          call map_larmor_plot(.7_fp,1.3_fp,61,-.3_fp,.3_fp,61)
       else
          call map_larmor_plot(0._fp,zr,61,-zr/2,zr/2._fp,61)
       end if
       idiagnostic = idiag
       write(*,*) 'enter anything'
       read(*,*) zr
    end if
  end subroutine map_larmor
 
  !
  ! ============================== map_larmor_loss ======================================
  ! returns ilim=0 if the particle's guiding center + larmor radius  is inside the limiter
  ! and ilim=1 if outside the limiter
  !
  ! the trcom variables gflr_rl and gflr_xv(2) are used to take short cuts:
  !   - the particle is ok if the nearest limiter is > gflr_rl*(the guiding center larmor radius)
  !   - the particle is ok if within the nearest limiter, the x equation is satisfied to within
  !     gflr_xv(1) and the v equation is greater than gflr_xv(1)>0.
  !
  subroutine map_larmor_loss(ilim)
    use nbi_com
 
    integer, intent(out) :: ilim    ! 0-> particle is ok, 1-> particle is lost
 
    integer :: iopt                 ! xplasma option to get nearest limiter point
    integer :: ierr                 ! xplasma return code
    integer :: jr(1)                ! region code
    integer :: ioops                ! error  code from map_larmor
    integer :: iregion              ! region code from map_larmor
    integer :: ifail(1)             ! nonzero on equation failure
 
    real(fp) :: zz_r(1), zz_z(1), zz_phi(1)  ! vector inputs to xplasma for guiding center r,z
    real(fp) :: zdist(1)                     ! distance to nearest limiter in meters (-for plasma side)
    real(fp) :: zz_rlim(1), zz_zlim(1), zz_philim(1) ! nearest limiter point
    real(fp) :: zzdiff(3)                    ! x,y,z vector from guiding center to nearest limiter
    real(fp) :: zxi, zperp                   ! projections of delta vector on gyro plane
    real(fp) :: zz                           ! distance to limiter in gyro plane from guiding center
    real(fp) :: zdist_gc                     ! signed distance to limiter from guiding center
    real(fp) :: zzcospho, zzsinpho           ! save gyro angle
    real(fp) :: za(1), zb(1)                 ! displacement in xi,par plane in larmor radii
    real(fp) :: zxeq(1), zveq(1)             ! x,v equations at point
    real(fp) :: zxp(3,1), zvec(3,1)          ! map_larmor_points data not used
    real(fp) :: zr, zy, zvpar, zvperp        ! map_larmor data not used
    real(fp) :: alpha_v    
 
    ilim = 0    ! assume particle is ok
    ic_loss(1) = ic_loss(1) + 1
 
    ! -- from limlos --
    if(xion.lt.xiliml) return
    if((rlar.eq.0.0_fp).and.(xion.lt.1.0_fp)) return
 
    ! --- check gc data ---
    if (imap_larmor_badgc) call map_larmor_gc    ! need the guiding center data
    if (imap_larmor_badgc) then
       ilim=1          ! gc is outside the r,z grid??
       return
    end if
 
    !
    ! -- get distance to nearest limiter from guiding center
    !
    zz_r(1)   = rmjion/100._fp  ! guiding center r,z,phi coordinates
    zz_z(1)   = yion/100._fp
    zz_phi(1) = 0._fp
 
    ic_loss(2) = ic_loss(2) + 1
    iopt=1     ! we want the point returned
    call eq_rzdist(1, zz_r, zz_z, zz_phi, zdist, iopt, zz_rlim, zz_zlim, zz_philim, ierr)
    if (ierr /= 0) goto 20
    zdist_gc =100._fp* zdist(1) ! <0 if inside limiter
 
    !
    ! -- test for particle too far away from limiter --
    !
    zz = max(gflr_ll, gflr_rl)
    if (zdist_gc < -zz*rlar) return                    ! too many larmor radii to
                                                       ! limiter
    !
    ! find gyro angle based on projection of closest limiter point
    ! onto xi,perp gyro plane -- should be near zero degrees
    !
    zzdiff = (/ 100._fp*zz_rlim(1)-rmjion, 0._fp, 100._fp*zz_zlim(1)-yion/)  ! vector from gc to nearest limiter
    zxi   = sum(zzdiff*zxi_gc)                            ! projection on xi_hat
    zperp = sum(zzdiff*zperp_gc)                          ! projection on perp_hat
    zz = sqrt(zxi*zxi+zperp*zperp)                        ! projected distance in gyro plane
    if (zz==0._fp.or. zdist_gc>0._fp) then
       ilim=1                          ! guiding center is on or outside the limiter
       return                          ! an impossibility so of course it happened
    end if
 
    zzcospho = cospho                  ! restore later, all exits through label 10 from here on
    zzsinpho = sinpho
    cospho = zxi/zz
    sinpho = zperp/zz
    zm_gc = cospho*zxi_gc + sinpho*zperp_gc  ! needed by map_larmor_points
 
    !
    !  -- relatively quick test at a limiter displacement --
    ! the particle is good if
    !   - the displacement kept the particle within the limiter (should always be true)
    !   - the x equation (measurement of b field orthogonality to gyro) < gflr_xv(1)
    !   - the v equation (energy conservation = 0. at gyro location)    > gflr_xv(2)>0.
    !       v equation is quadratic with a positive second derivative and negative minimum
    !
    ! the result of this test may be used later on
    !
    ic_loss(3) = ic_loss(3) + 1
    za(1) = 0.98_fp*zz/max(zscale_gc,1.e-10_fp*raxis)  ! displace towards limiter with almost limiter distance
    zb(1) = 0.0_fp
    call map_larmor_points(1, za, zb, jr, zxeq, zveq, zxp, zvec, ifail=ifail)
    !
    ! can't trust region code because of map_larmor_bffget so find out
    ! if this point is outside the limiter and set region code jr(1)
    !
    if (jr(1)<=3 .and. ifail(1)==0) then
       zz_r(1) = sqrt(zxp(1,1)**2+zxp(2,1)**2)/100._fp
       zz_z(1) = zxp(3,1)/100._fp
       iopt=0                   ! we don't want the point returned
       call eq_rzdist(1, zz_r, zz_z, zz_phi, zdist, iopt, zz_rlim, zz_zlim, zz_philim, ierr)
       if (ierr /= 0) goto 20
       if (zdist(1)>0) jr(1)=3     ! outside the limiter
       !debug write(77,'(f12.4,2x,f12.4,2x,f12.4)') rlar,zz,100.*zdist
    end if
 
    if (ifail(1)/=0) goto 10   ! energy equation is not even valid so particle is considered ok

    if (jr(1)<3 .and. abs(zxeq(1))<gflr_xv(1) .and. zveq(1)>max(0._fp,gflr_xv(2))) goto 10  ! particle is ok
 
    !
    ! -- still don't know, find actual location of particle --
    !
    ic_loss(4) = ic_loss(4) + 1
    call map_larmor(ioops, iregion, zr, zy, zvpar, zvperp, alpha_v ,.false.)
 
    if (ioops==1 .or. ioops==2) then   ! got classical larmor displacment or guiding center
       if (jr(1)<4) then
          ! nothing is working so now we assume that the position found
          ! in the previous step is at the limiter (though it was only
          ! at 0.98 of the distance) and directly towards it and use
          ! the energy conservation equation (v equation) and ignore b
          ! field orthogonality (x equation) to judge whether the particle
          ! is inside the limiter or not -- in other words, guess ...
          ic_loss(5) = ic_loss(5) + 1
          if (zveq(1)<0._fp) ilim=1   ! lost the particle
          write(nout,*) '! map_larmor_loss: flr estimated for particle loss'
          goto 10
       else
          ic_loss(6) = ic_loss(6) + 1
          write(nout,*) '! map_larmor_loss: warning, using poor flr information for particle loss'
          if (ioops==1) then
             write(nout,*) '! map_larmor_loss: loss based on classical flr displacement'
          else
             write(nout,*) '! map_larmor_loss: loss based on guiding center'
          end if
       end if
    end if
 
    if (ioops/=2) then
       ! -- have a displacement, is this inside or outside the limiter --
       zz_r(1) = zr/100._fp
       zz_z(1) = zy/100._fp
       iopt=0                   ! we don't want the point returned
       call eq_rzdist(1, zz_r, zz_z, zz_phi, zdist, iopt, zz_rlim, zz_zlim, zz_philim, ierr)
       if (ierr /= 0) goto 20
       if (zdist(1)>0) ilim=1   ! outside the limiter
    else
       ! -- guiding center was returned --
       if (-zdist_gc<rlar) ilim=1   ! within one larmor radius of guiding center
    end if
 
10  cospho = zzcospho      ! normal exit -- restore gyro angle
    sinpho = zzsinpho
    return
 
20  continue               ! error exit for eq_rzdist
    write(nout,*) '? map_larmor_loss: error in call to eq_rzdist'
    call map_larmor_info
    call nbi_exit
  end subroutine map_larmor_loss
 
  !
  ! ============================== map_larmor_plot ======================================
  !
  subroutine map_larmor_plot(zamin,zamax,ina,zbmin,zbmax,inb)
    !
    ! do some plotting -- really a junk routine so check before using
    !
    use nbi_com
 
    integer, intent(in) :: ina, inb   ! number of test points in a,b directions
    real(fp), intent(in) :: zamin, zamax  ! extent of coverage of test points orthogonal to b field
                                      ! in units of ideal larmor radius
    real(fp), intent(in) :: zbmin, zbmax  ! extent of coverage of test points along b field
                                      ! in units of ideal larmor radius
 
    logical i_close                   ! closeness flag
 
    integer :: ia, ib                 ! loop variables
    integer :: i,j, io                ! temporary
    integer :: jregion(ivecd)         ! xplasma code for region of point
 
    real(fp) :: zda, zdb                  ! step size in a,b direction
 
    real(fp) :: za(ivecd),zb(ivecd)       ! displacement in b plane
    real(fp) :: zxeq(ivecd)               ! dimensionless  delta_x . b ---> should be zero
    real(fp) :: zveq(ivecd)               ! dimensionless  v_from_deltax**2 - v_from_e_conservation**2
                                      ! which should be zero
    real(fp) :: zvec(3,ivecd)             ! velocity vector
    real(fp) :: zxp(3,ivecd)              ! position vector
 
    ! -- grab data variables --
    integer i1,i2
    real(fp)    :: z1, z2, zrho(ivecd), zchi(ivecd), zphi(ivecd), ztens(3,3,ivecd)
    real(fp)    :: zp2,zt2,zb2,zpg(3),ztg(3),zbg(3),zppg(3),zppgmag,zbbgmag,zttgmag
 
    ! -- plot variables --
    integer, parameter :: imaxab = 201
    integer, parameter :: ikolors(6) = (/2,3,4,5,6,7/)
    integer, parameter :: isymbs(6)=(/2,3,4,5,6,7/)
    integer, parameter :: imaxlevl=200
    integer :: igrid, ilevl, incl, iconst
    character*55 :: title
    real(fp) :: zout_xeq(imaxab,imaxab), zout_veq(imaxab,imaxab)
    real(fp) :: za_data(imaxab), zb_data(imaxab)
    real(fp) :: zxmin, zxmax, zxinc, zvmin, zvmax, zvinc
    real(fp) :: zcl(imaxlevl)
    integer :: ixl, ivl, io70, io71, io72, io74, ilun_list(4)
 
    !
    ! ------- check -------
    !
    if (imap_larmor_badgc) then
       call map_larmor_gc    ! update guiding center data
       if (imap_larmor_badgc) call nbi_exit
    else
       zm_gc = cospho*zxi_gc + sinpho*zperp_gc  ! repeat in case cospho,sinpho changed
    end if

    call find_io_unit_list(4,ilun_list)  ! never tested
    io70 = ilun_list(1)
    io71 = ilun_list(2)
    io72 = ilun_list(3)
    io74 = ilun_list(4)

    ! ---- grab some data ---
    if (idiag_test) then
       i1=141
       i2=281
       write(io74,'(a)') '#    r              z              xi            th               1/r          1/dl' &
            // '           1/omegal      bt/b             b          xgrad(bt/b)    ygrad(bt/b)    zgrad(bt/b)' &
            // '          xgrad(b)      ygrad(b)        zgrad(b)'
       do i=1,i1
          z1 = 20._fp+ 140._fp*real(i-1,fp)/real(i1-1,fp)
          do j=1,i2
             z2 = -140 + 280._fp*real(j-1,fp)/real(i2-1,fp)
 
             zz_xx(1) = z1/100._fp
             zz_yy(1) = 0._fp
             zz_zz(1) = z2/100._fp
             call eqxyzmap_bgrad(1, zz_xx, zz_yy, zz_zz, zrho, zchi, zphi, ztol, &
                  zbvec, ztens, jregion, ierr)
             if (jregion(1)<4) then
                zp2 = zbvec(1,1)**2+zbvec(3,1)**2
                zt2 = zbvec(2,1)**2
                zb2 = zp2 + zt2
                zpg = (zbvec(1,1)*ztens(:,1,1) + zbvec(3,1)*ztens(:,3,1))/zp2
                ztg = ztens(:,2,1)/zbvec(2,1)
                zbg = (zbvec(1,1)*ztens(:,1,1) + zbvec(2,1)*ztens(:,2,1) &
                     + zbvec(3,1)*ztens(:,3,1))/zb2
                zppg = (zp2/zb2)*(ztg-zpg)
                zppgmag = sqrt(sum(zppg*zppg))  ! in inverse meters from xplasma call
                zttgmag = sqrt(sum(ztg*ztg))
                zbbgmag = sqrt(sum(zbg*zbg))
                write(io74,'(15(g14.7,1x))') z1, z2, zrho(1), zchi(1)*180._fp/3.1415926_fp, &
                     1._fp/z1, zppgmag/100._fp, zbbgmag/100._fp, &
                     sqrt(zt2/zb2), sqrt(zb2), zppg/zppgmag, zbg/zbbgmag
             end if
          end do
          write(io74,*) '   '
       end do
       call nbi_exit
    end if
 
    if (ina>imaxab .or. inb>imaxab) then
       write(nout,*) '? map_larmor: number of points exceeded the dimensions'
       write(nout,*) ' na  = ', ina
       write(nout,*) ' nb  = ', inb
       write(nout,*) ' dim = ', imaxab
       call nbi_exit
    end if
 
    !
    ! ----------- loop over points in guiding center plane ------------
    !
    write(io70,*) ' '
    write(io70,*) '# xion  = ', xion
    write(io70,*) '# th    = ', th
    write(io70,*) '# ang   = ', (180._fp/3.14159265359_fp)*atan2(sinpho,cospho)
    write(io70,*) '# scale = ', zscale_gc
    write(io70,*) '#   b           a             xeq         veq'
    zda = (zamax-zamin)/max(real(ina-1),1._fp)
    zdb = (zbmax-zbmin)/max(real(inb-1),1._fp)
    zxmin = 0._fp
    zxmax = 0._fp
    zvmin = 0._fp
    zvmax = 0._fp
    do ib = 1, inb
       zb(1) = zbmin + (ib-1)*zdb
       zb_data(ib) = zb(1)
 
       do ia = 1, ina
          za(1) = zamin + (ia-1)*zda
          za_data(ia) = za(1)
 
          call map_larmor_points(1, za, zb, jregion, zxeq, zveq, zxp, zvec)
 
          zout_xeq(ia,ib) = zxeq(1)
          zout_veq(ia,ib) = zveq(1)
          write(io70,'(4(g13.6,1x))') zb(1), za(1), zxeq(1), zveq(1)
 
          if (zxeq(1)<zxmin) zxmin = zxeq(1)
          if (zxeq(1)>zxmax) zxmax = zxeq(1)
          if (zveq(1)<zvmin) zvmin = zveq(1)
          if (zveq(1)>zvmax) zvmax = zveq(1)
 
       end do  ! end loop over a points
       write(io70,*) ' '
    end do     ! end loop over b points
    write(io70,*) ' '
 
    ! -- contour limits --
    zxmin = (int(zxmin*10._fp)-1)/10._fp
    zxmax = (int(zxmax*10._fp)+1)/10._fp
    zvmin = (int(zvmin*10._fp)-1)/10._fp
    zvmax = (int(zvmax*10._fp)+1)/10._fp
    ixl = min(imaxlevl,int(10._fp*(zxmax-zxmin)+1.001_fp))
    ivl = min(imaxlevl,int(10._fp*(zvmax-zvmin)+1.001_fp))
 
    zxmin = -1._fp
    zxmax =  1._fp
    zvmin = -1._fp
    zvmax =  1._fp
    ixl = 21
    ivl = 21
 
    !
    ! --------- plot ---------
    !
!%!    call initt(0)
!%!    call erase
!%!    call binitt
!%!    call slimx(100,700)
!%!    call slimy(100,700)
!%!    call agclini(0,0,0)
!%!    call agcolrs(6,ikolors,0)
!%!    call agcsymb(6,isymbs,0)
    igrid=2
    ilevl=8
!%!    call agcfram(za_data, ina, zb_data, inb, igrid)
 
    ! --- x equation contours ---
    i = 1
    do ia=1,ixl
       zcl(ia) = zxmin+((ia-1._fp)/10._fp)
       if (abs(zcl(ia))<.01_fp) i = ia
    end do
    zcl(i) = zcl(ixl)
    zcl(ixl) = 0._fp
    incl=ixl
!%!    call rcontr(incl, zcl, ixl, zout_xeq, imaxab, za_data, 1, ina, 1, zb_data, 1, inb, 1)
 
    ! --- v equation contours ---
    i = 1
    do ia=1,ivl
       zcl(ia) = zvmin+((ia-1._fp)/10._fp)
       if (abs(zcl(ia))<.01_fp) i = ia
    end do
    zcl(i) = zcl(ivl)
    zcl(ivl) = 0._fp
    incl=ivl
!%!    call rcontr(incl, zcl, ivl, zout_veq, imaxab, za_data, 1, ina, 1, zb_data, 1, inb, 1)
 
    !
    ! --- iterations ---
    !
    ! newton
    !
    za(1) = 1.0_fp
    zb(1) = 0.0_fp
    j = min(ivecd,int(gflr_op(6)+.01_fp))
    call map_larmor_info
    print *, ' newton '
    call map_larmor_newton(j, .true., gflr_op(3), gflr_op(2), za, zb, jregion, zxeq, zveq, &
         zxp, zvec, i_close)
    write(io71,*) ' '
    write(io71,'(a,f10.6,a,f8.3,a,f8.3)') '#  xi = ', xion, '    th = ', th*(180/xpi), '  angle = ', &
         (180._fp/3.14159265359_fp)*atan2(sinpho,cospho)
    write(io71,'(a,f12.6,a,f11.7,a,f12.6)') '#  e(kev)= ', vion*vion*v2toev/1000._fp, '  xksid = ', xksid, &
         '  rlar = ', zrlar_gc
    write(io71,*) '#   k    a           b             xeq         veq'
    write(*,*) 'close = ', i_close
    do i = 1, j
       write(*,'(i2,1x,4(g13.6,1x))') i, za(i), zb(i), zxeq(i), zveq(i)
       write(io71,'(i2,1x,4(g13.6,1x))') i, za(i), zb(i), zxeq(i), zveq(i)
!%!       call movea(za(i), zb(i))
!%!       call bsyms(za(i), zb(i),isymbs(1))
    end do
 
    !
    ! secant
    !
    za(1) = 1.0_fp
    zb(1) = 0.0_fp
    j = min(ivecd,int(gflr_op(6)+.01_fp))
    print *, ' secant '
    call map_larmor_newton(j, .false., gflr_op(3), gflr_op(2), za, zb, jregion, zxeq, zveq, &
         zxp, zvec, i_close)
    write(io72,*) ' '
    write(io72,'(a,f10.6,a,f8.3,a,f8.3)') '#  xi = ', xion, '    th = ', th*(180/xpi), '  angle = ', &
         (180._fp/3.14159265359_fp)*atan2(sinpho,cospho)
    write(io72,'(a,f12.6,a,f11.7,a,f12.6)') '#  e(kev)= ', vion*vion*v2toev/1000._fp, '  xksid = ', xksid, &
         '  rlar = ', zrlar_gc
    write(io72,*) '#   k    a           b             xeq         veq'
    write(*,*) 'close = ', i_close
    do i = 1, j
       write(*,'(i2,1x,4(g13.6,1x))') i, za(i), zb(i), zxeq(i), zveq(i)
       write(io72,'(i2,1x,4(g13.6,1x))') i, za(i), zb(i), zxeq(i), zveq(i)
!%!       call movea(za(i), zb(i))
!%!       call bsyms(za(i), zb(i),isymbs(2))
    end do
 
    ib = 744  ! a nice position for the first line of the title
    write(title,'(a,f10.6,a,f8.3,a,f8.3)') 'flr at  xi=', xion, ' th=', th*(180/xpi), ' angle=', &
         (180._fp/3.14159265359_fp)*atan2(sinpho,cospho)
!%!    call agtcgt(ib,title)
 
    ib = 725  ! a nice position for the second line of the title
    write(title,'(a,f12.6,a,f11.7,a,f12.6)') 'e(kev)=', vion*vion*v2toev/1000._fp, ' xksid=', xksid, &
         ' rlar=', zrlar_gc
    !write(title, '(a,f6.4,a,f8.3,a,f8.3)') 'flr at xi=',xion,' th=',th*180./3.1415926, &
    !     ' ang=',atan2(sinpho,cospho)*180./3.1415926
!%!    call agtcgt(ib,title)
 
!%!    call agtcvs('parallel distance (larmor)')
!%!    call agtchs('perpendicular distance(larmor)')
 
!%!    call upause(0)
!%!    call finitt(0,0)
  end subroutine map_larmor_plot
 
  !
  ! ============================== map_larmor_newton ======================================
  !
  subroutine map_larmor_newton(istep, i_newton, zconv_ab, zconv_xv, za_arg, zb_arg, &
       jcode, zx_arg, zv_arg, zxp_arg, zvec_arg, i_close)
    !
    ! given a starting point za_arg(1), zb_arg(1), take (istep-1) more steps
    ! using a newton iteration scheme for the first and secant or newton steps for
    ! the remainder.  return in istep the last valid data point.  istep could be
    ! zero if the first point is bad.
    !
    ! also stop when the change in a,b is less then zconv_ab or the magnitude of
    ! the x and v equations is less than zconv_xv
    !
    ! the i_close return variable indicates that the algorithm did not converge but
    ! a v equation bracket was performed on the last iterated point and the 'a' variable
    ! in the bracket was not too close to the boundary.  so while the point may not
    ! satisfy convergence, it is probably better then a classical larmor offset.
    !
    use nbi_com
 
    integer, parameter :: kmax=6      ! number of iterations to find a metric minimum in
                                      ! step advancement
    real(fp) :: zdelta_factor_start=0.2_fp  ! initial control for the bracketing width
 
    logical, intent(in) :: i_newton   ! .true. to take only newton steps
    logical, intent(out) :: i_close   ! .true. when the algorithm failed but the returned point is
                                      ! is a bracket on the v equation
    integer, intent(inout) :: istep   ! number of valid points requested and returned
    integer, intent(out) :: jcode(:)  ! region code for the iterated points
 
    real(fp), intent(in)    :: zconv_ab             ! terminate the loop if change in a,b is too small
    real(fp), intent(in)    :: zconv_xv             ! terminate the loop if the x,v equations are small
    real(fp), intent(inout) :: za_arg(:), zb_arg(:) ! iterated points, first one is input
    real(fp), intent(out)   :: zx_arg(:), zv_arg(:) ! value of x,v equation at each step
    real(fp), intent(out)   :: zxp_arg(3,*)         ! value of position vector at each step
    real(fp), intent(out)   :: zvec_arg(3,*)        ! value of velocity vector at each step
 
    logical :: i_newton_step  ! .true. when the jacobian is found from
                              ! a finite difference approximation
    logical :: i_center       ! .true. if the center point is needed in evaluating x,v equations
    logical :: i_vbracket     ! .true. if a v equation bracket has been performed
    logical :: i_search_fail  ! .true. when a search fails
    logical :: ilast          ! last chance bracketing flag
    logical :: i_bad          ! .true. when v bracket returned is not within r,z grid
 
    integer :: i_exceed       ! number of times a or b went out of bounds
    integer :: i_start        ! number of restarts for this iteration
    integer :: k              ! loop variable
    integer :: k1,km          ! k+1,k-1
    integer :: kstep_ok       ! number of successive successful step advancements
    integer :: i_status       ! return code from step advancement
    integer :: istep_max      ! the initially requested size of istep -- maximum number of points
    integer :: kcheck(3)      ! last 3 step advancements in a row
    integer :: klast          ! step at which to start testing for convergence rate
 
    real(fp) :: zmax_larmor       ! don't let a,b get beyond this
    real(fp) :: zmax_larmor_test  ! algorithm will start to fail with a,b beyond this
    real(fp) :: zmicro            ! last successful step size in micro search iteration
    real(fp) :: zfinite_diff_step ! finite difference step size in a,b
    real(fp) :: zza, zzb          ! temporary holding
 
    real(fp) :: zzfda, zzfdb      ! finite difference step size
    real(fp) :: zzx, zzv          ! x and v equation at the point
    real(fp) :: zzj(2,2)          ! jacobian at the point d(x,v)/d(a,b)
    real(fp) :: zrowx, zrowv      ! temporary
    real(fp) :: zda, zdb          ! change in a,b in this iteration
    real(fp) :: zz,zf             ! temporary
    real(fp) :: zdelta            ! v equation bracketing width
    real(fp) :: zdelta_factor     ! dynamically controls the bracketing width
    real(fp) :: zlast             ! criteria for slow convergence of v function
 
    zmax_larmor_test = gflr_ll-0.5_fp  ! warning limit -- keep offset bigger then zdelta_factor_start
    zmax_larmor      = gflr_ll            ! hard limit
    zmicro = 0.8_fp*zmicro_step              ! start with moderate size step
 
    zfinite_diff_step = .05_fp/max(.1_fp,zscale_gc)  ! on order of 5cm scale size * 1e-2
    zzfda = zfinite_diff_step                  ! finite difference step size
    zzfdb = zfinite_diff_step
 
    istep_max = istep      ! istep will be changed, don't change the loop limit
    i_newton_step = .true. ! need the jacobian on first step
    i_center = .true.      ! need the center point data on the first step
    i_exceed = 0           ! number of excursions beyond (a,b) boundaries
    i_vbracket = .false.   ! allow v bracketing
    ilast      = .false.   ! last chance bracketing flag
    zlast      = max(.01_fp,10._fp*zconv_xv)   ! convergence rate criteria
    klast      = max(5,(2*istep_max/3))  ! start checking convergence
    kstep_ok = 0
    kcheck   = 0
    zdelta_factor = zdelta_factor_start
    i_close = .false.      ! in case of failure
    do k = 1, istep_max    ! loop over the points, but don't fill past istep<=istep_max
       !
       ! evaluate x,v equations at this point, also jacobian
       ! if not the last point
       !
       i_start = 0
10     continue  ! restart point
 
       !
       ! -- must have a>0 and a,b bounded --
       !
       if (za_arg(k)<zconv_ab .or. za_arg(k)>zmax_larmor_test) then
          za_arg(k) = min(zmax_larmor,max(zconv_ab,za_arg(k)))
          i_exceed  = i_exceed+1
       end if
 
       if (zb_arg(k)<-zmax_larmor_test .or. zb_arg(k)>zmax_larmor_test) then
          zb_arg(k) = min(zmax_larmor,max(-zmax_larmor,zb_arg(k)))
          i_exceed  = i_exceed+1
       end if
 
       !
       ! -- check on bounds and excessive starting --
       !
       i_start=i_start+1
       if (i_exceed>3 .or. i_start>7) then
          if (idiagnostic)  print *, 'map_larmor_newton: error exit i_start=', i_start, ' i_exceed=', i_exceed
          istep = k-1       ! last point is the best
          return
       end if
 
       ! -- get data --
       if (i_newton_step) then
          call map_larmor_jacobian(i_center, za_arg(k), zb_arg(k), zzfda, zzfdb, jcode(k), &
               zx_arg(k), zv_arg(k), zxp_arg(1,k), zvec_arg(1,k), zzj)
       else if (i_center) then
          ! -- don't need the jacobian but do need the point --
          call map_larmor_points(1, za_arg(k:k), zb_arg(k:k), jcode(k:k), zx_arg(k:k), &
               zv_arg(k:k), zxp_arg(1,k), zvec_arg(1,k))
       end if
 
       zzx = zx_arg(k)       ! simplify the equation below
       zzv = zv_arg(k)
 
       !
       ! -- checks --
       !
       if ((abs(zzx)<zconv_xv) .and. (abs(zzv)<zconv_xv) .and. jcode(k)<4) then
          ! -- reached convergence in x,v equations --
          istep = k
          return       ! done ok
       else if (k==istep_max) then
          if (jcode(k)<4) then
             istep = k
          else
             istep = k-1
          end if
          goto 100       ! not ok but ran out of iterations
       else if (jcode(k)>=4) then
          !
          ! -- point is outside the r,z grid, start closer to guiding center --
          !
          if (idiagnostic) print *, 'map_larmor_newton: point is out of bounds, moving inwards'
          za_arg(k) = 0.7_fp*za_arg(k)
          zb_arg(k) = 0.7_fp*zb_arg(k)
          i_newton_step = .true.      ! need the jacobian
          i_center = .true.           ! need the center point
          goto 10                     ! restart iteration
       else if (k==1) then
          if (jcode(k)<4 .and. zv_arg(k)<0._fp) then
             !
             ! -- check for wrong side of minimum on first step --
             ! the v equation is dominantly quadratic with a positive second derivative (i think)
             ! so make sure we are on the correct side of the minimum otherwise may end up
             ! converging to a=0,b=0 and never satisfying the convergence condition
             !
             zz = za_arg(k)*zzj(2,1)+zb_arg(k)*zzj(2,2)  ! in direction of  (a,b) dot gradv
             if (zz<0) then
                if (.not. i_vbracket) then
                   if (idiagnostic) print *, 'map_larmor_newton: attempting to bracket beyond minimum'
                   call map_larmor_vbracket(za_arg(k), zb_arg(k), zzv, 10._fp*zconv_ab, &
                        zmax_larmor, zdelta_factor_start, 2.19_fp, i_bad)
                   i_vbracket = .true.
                   goto 10                   ! restart first iteration
                else if (za_arg(k)>zmax_larmor_test) then
                   if (idiagnostic) print *, 'map_larmor_newton: wrong side of minimum at edge of boundary'
                   istep = k-1
                   return
                end if
             end if
          end if
       end if
 
       !
       ! -- the k point is within the r,z grid (jcode(k)<4) --
       !
       !
       ! ---- get jacobian for secant method ---
       !
       if (.not. i_newton_step) then
          km = k-1                        ! k=2,3,...
          zda = za_arg(k) - za_arg(km)    ! delta_x
          zdb = zb_arg(k) - zb_arg(km)
          zz  = zda*zda+zdb*zdb
 
          zrowx =  (zx_arg(k) - zx_arg(km)) - (zzj(1,1)*zda+zzj(1,2)*zdb)  ! deltaf - j*delta_x
          zrowv =  (zv_arg(k) - zv_arg(km)) - (zzj(2,1)*zda+zzj(2,2)*zdb)
          if (zz /= 0._fp) then
             zzj(1,1) = zzj(1,1) + zrowx*zda/zz
             zzj(2,1) = zzj(2,1) + zrowv*zda/zz
             zzj(1,2) = zzj(1,2) + zrowx*zdb/zz
             zzj(2,2) = zzj(2,2) + zrowv*zdb/zz
          end if
       end if
 
       !
       ! ------- attempt to advance the point -------
       !
       k1 = k+1     ! advance to next point
       za_arg(k1) = za_arg(k)
       zb_arg(k1) = zb_arg(k)
       zx_arg(k1) = zx_arg(k)
       zv_arg(k1) = zv_arg(k)
       call map_larmor_step(k, za_arg(k1), zb_arg(k1), zx_arg(k1), zv_arg(k1), zzj,&
            zconv_ab, zconv_xv, kmax, i_status, jcode(k1), zxp_arg(1,k1), zvec_arg(1,k1),&
            gflr_op(5))
 
       !
       ! a series of advancements which requires a lot of iterations
       ! can indicate being caught in a local minimum, try v bracketing
       !
       kcheck(1+mod(k,3)) = max(i_status,0)     ! store number of steps to convergence
       if (.not. i_vbracket) then
          if ((kcheck(1)+kcheck(2)+kcheck(3))>9) then
             if (idiagnostic) print *, 'map_larmor_newton: steps are too difficult'
             i_status=-1                ! force a v equation bracket
          else if ((.not. ilast) .and. i_status>0 .and. (k >= klast)&
               .and. abs(zv_arg(k1))>zlast) then
             if (idiagnostic) print *, 'map_larmor_newton: converging too slow'
             zdelta_factor = zdelta_factor_start/20._fp  ! go for a tight bracket
             ilast = .true.
             i_status=-1                ! force a v equation bracket
          end if
       end if
 
       !
       ! --------- deal with result of advancement -------
       !
       if (i_status>=0) then
          ! -- step was ok --
          i_center      = .false.    ! center point data came from map_larmor_step
          i_newton_step = i_newton   ! if secant, don't need the jacobian
          if (i_status<3) then
             kstep_ok = kstep_ok + 1             ! only 1 or 2 iterations in step advancement
             if (kstep_ok>1) i_vbracket=.false.  ! enable v equation bracketing after 2 good steps
          end if
       else
          ! -- failure --
          kcheck   = 0                          ! reset the checking
          kstep_ok = 0                          ! bad step
          za_arg(k1) = za_arg(k)                ! start at point before advancement
          zb_arg(k1) = zb_arg(k)
          i_center      = .true.                ! need the data at the next point
          i_newton_step = .true.                ! restart the jacobian
 
          if (idiagnostic) print *, 'map_larmor_newton: step failure at step k=', k
          if (.not. i_vbracket) then
             !
             ! this section attempts to restart the search near the v=0 contour,
             ! growth rate reflects 2.19**3= 10.5
             !
             zdelta = zdelta_factor*max(.1_fp, min(1._fp,za_arg(k)))
             if (idiagnostic) print *, 'map_larmor_newton: attempting v bracket'
             call map_larmor_vbracket(za_arg(k1), zb_arg(k1), zv_arg(k), 10._fp*zconv_ab, &
                  zmax_larmor, zdelta, 2.19_fp, i_bad)
             if (i_bad) then
                if (idiagnostic) print *, 'map_larmor_newton: v bracket ended at edge of r,z grid, abort'
                istep = k
                return
             end if
             zdelta_factor = zdelta_factor_start  ! reset to original
             i_vbracket = .true.
          else
             !
             ! -- last hope, search for a better point --
             !
             if (za_arg(k)>zmax_larmor_test) then
                if (idiagnostic) print *, 'map_larmor_newton: too close to edge of boundary for searching'
                istep = k
                return
             end if
             if (idiagnostic) print *, 'map_larmor_newton: attempting search'
             zf = (zx_arg(k)**2+zv_arg(k)**2)/2._fp  ! metric at k point
             call map_larmor_search(za_arg(k1), zb_arg(k1), zf, zmicro, i_search_fail)
             if (i_search_fail) then
                if (idiagnostic) print *, 'map_larmor_newton: search failure'
                istep = k              ! could not advance
                goto 100
             end if   ! end of search failure test
          end if      ! end of vbracket/search test
       end if         ! end of step advancement result test
 
    end do
    ! code should not pass through here
 
100 continue
    ! --- failed to get convergence, try to get close by returning a v bracket ---
    k  = istep                                ! last point which would have been returned
    k1 = min(k+1, istep_max)                  ! if room, put in next point spot
    za_arg(k1) = za_arg(k)                    ! start at last good point
    zb_arg(k1) = zb_arg(k)
    zza = za_arg(k1)                          ! save in case it does not work
    zzb = zb_arg(k1)
    zdelta_factor = zdelta_factor_start/20._fp  ! go for a tight bracket
    zdelta = zdelta_factor*max(.1_fp, min(1._fp,za_arg(k)))
 
    ! -- bracket --
    call map_larmor_vbracket(za_arg(k1), zb_arg(k1), zv_arg(k), 10._fp*zconv_ab, &
         zmax_larmor, zdelta, 2.19_fp, i_bad)
 
    ! -- get data --
    call map_larmor_points(1, za_arg(k1:k1), zb_arg(k1:k1), jcode(k1:k1), zx_arg(k1:k1), &
         zv_arg(k1:k1), zxp_arg(1,k1), zvec_arg(1,k1))
 
    if (i_bad .or. jcode(k1)>3 .or. za_arg(k1)>(zmax_larmor_test - 2._fp*zdelta) &
         .or. za_arg(k1)<(10._fp*zconv_ab+2._fp*zdelta)) then
       if (idiagnostic) print *, 'map_larmor_newton: unable to return a v bracketed point'
       i_close = .false.   ! maybe too close to edge, probably not valid
    else if ( abs(zx_arg(k1))>gflr_op(8) ) then
       ! exceeded the x equation bracketing, maybe jumped away from the
       ! ideal point due to multivaluedness in the v equation.  see if
       ! last good point is good enough
       za_arg(k1) = zza                    ! original last point
       zb_arg(k1) = zzb
       call map_larmor_points(1, za_arg(k1:k1), zb_arg(k1:k1), jcode(k1:k1), zx_arg(k1:k1), &
            zv_arg(k1:k1), zxp_arg(1,k1), zvec_arg(1,k1))
       if (jcode(k1)<4 .and. abs(zx_arg(k1))<gflr_op(8) .and. abs(zv_arg(k1))<gflr_op(9)) then
          if (idiagnostic) print *, 'map_larmor_newton: bad v bracket but last point is close enough'
          istep   = k1
          i_close = .true.
       else
          if (idiagnostic) print *, 'map_larmor_newton: bad v bracket and last point is not close'
          i_close = .false.
       end if
    else
       ! the x equation is good enough and the v equation is bracketed on
       ! a fine enough interval so this point is close
       if (idiagnostic) print *, 'map_larmor_newton: algorithm failed but returning a v bracketed point'
       istep   = k1
       i_close = .true.
    end if
 
    return
 
  end subroutine map_larmor_newton
 
  !
  ! ============================== map_larmor_info ======================================
  !
  subroutine map_larmor_info
    !
    ! debugging output
    !
    use nbi_com
 
    write(nout,'(a,f10.6,a,f8.3,a,f8.3)') '  xi = ', xion, '    th = ', th*(180/xpi), '  angle = ', &
         (180._fp/3.14159265359_fp)*atan2(sinpho,cospho)
    write(nout,'(a,f12.6,a,f11.7,a,f12.6)') '  e(kev)= ', vion*vion*v2toev/1000._fp, '  xksid = ', xksid, &
         '  rlar = ', zrlar_gc
  end subroutine map_larmor_info
 
  !
  ! ============================== map_larmor_vbracket =============================
  !
  subroutine map_larmor_vbracket(za, zb, zv, zamin, zamax, zdelta, zgrow, i_bad)
    !
    ! this subroutine will test increasing/decreasing (a,b) points in order to bracket
    ! the zero crossing of the v equation and improve the initial za argument.
    !
    ! when the zero crossing is between (za, zamax) the intervals
    ! (za, za*zgrow), (za*zgrow, za*zgrow**2), ... will be searched
    !
    ! require zdelta>0. and zgrow>1. and (za,zb) in r,z grid -- this is not checked
    !
    integer, parameter  :: np=3          ! number of points to check per step
 
    logical, intent(out) :: i_bad        ! .true. if the upper bracket is not
                                         !  within the r,z grid
    real(fp), intent(inout) :: za, zb        ! initial and final starting point
    real(fp), intent(in)    :: zv            ! v equation at the starting point
    real(fp), intent(in)    :: zamin, zamax  ! limits for searching for bracket
    real(fp), intent(in)    :: zdelta        ! stop when bracketed by this amount
    real(fp), intent(in)    :: zgrow         ! must be >1., this is the growth rate for
                                         ! expanding the search -- suggest 10.**(1/3)
 
    logical :: i_change                  ! .true. when a bound has changed in this step
    integer :: i                         ! loop variable
    integer :: jr(np)                    ! xplasma region code
    real(fp) :: amin, amax                   ! current bounds on v crossing
    real(fp) :: samax                        ! the actual search maximum
    real(fp) :: agood                        ! the largest good point found so far
    real(fp) :: adelta                       ! step size between points
    real(fp) :: zsa(np), zsb(np)             ! xplasma points
    real(fp) :: zsxeq(np), v(np)             ! x, v equations from xplasma
    real(fp) :: zsxvec(3,np), zsvvec(3,np)   ! xplasma vectors
 
    ! --- set direction for search ---
    if (zv<0._fp) then
       amin  = za
       amax  = zamax
    else
       amin  = zamin
       amax  = za
    end if
 
    i_change = .true.
    samax    = za
    agood    = amin      ! assume this is good
 
    do while (((amax-amin)>zdelta) .and. i_change)
       samax = min(samax*zgrow, amax)    ! expand the upper bound of the search
       adelta = (samax-amin)/(np+1)
       do i = 1, np
          zsa(i) = amin + i*adelta       ! does not include amin or samax
          zsb(i) = zb
       end do
 
       ! get v equation array
       call map_larmor_points(np, zsa, zsb, jr, zsxeq, v, zsxvec, zsvvec)
 
       ! -- find largest good point --
       do i = 1, np
          if (jr(i)<4) agood = max(agood, zsa(i))
       end do
 
       ! -- get lower bound --
       i_change = .false.
       do i = 1, np
          if (jr(i)<4) then
             if (v(i)>0._fp) exit
             amin = zsa(i)       ! v(i)<=0. new lower bound
             i_change = .true.
          end if
       end do
 
       ! -- get upper bound --
       do i = np, 1, -1
          if (jr(i)<4 .and. v(i)<0._fp) exit
          amax = zsa(i)          ! v(i)>=0. or outside grid gives new upper bound
          i_change = .true.
       end do
 
    end do
    za = min((amin+amax)/2._fp, agood)  ! this point should be near v=0 and in the r,z grid
    i_bad = zv<0._fp.and. agood<amax
    if (i_bad) then
       ! see if amax is outside of r,z grid
       zsa(1) = amax
       zsb(1) = zb
       call map_larmor_points(1, zsa, zsb, jr, zsxeq, v, zsxvec, zsvvec)
       i_bad = jr(1)>3
    end if
  end subroutine map_larmor_vbracket
 
  !
  ! ============================== map_larmor_search ======================================
  !
  subroutine map_larmor_search(za_arg1, zb_arg1, zmetric, zmicro, i_search_fail)
    !
    ! this subroutine searches around the point za_arg1,zb_arg1 for a smaller
    ! x,v equation.  it's primary purpose is to attempt to make progress
    ! when the iteration scheme fails so the iteration can be restarted.
    ! the best point is returned in za_arg1, zb_arg1.
    !
    ! at za_arg1,zb_arg1, zmetric=(zx**2+zv**2)/2. required as input
    !
    ! hexs are searched starting at a radius of zmicro/zmicro_step then decreasing
    ! in zmicro_step steps
    !
    use nbi_com
 
    real(fp), intent(inout)  :: za_arg1, zb_arg1     ! start point and returned best point
    real(fp), intent(inout)  :: zmetric              ! metric to be improved on
    real(fp), intent(inout)  :: zmicro               ! last step size used in the micro iteration
    logical, intent(out) :: i_search_fail        ! return .true. on failure to find a better point
 
    integer, parameter :: idsearch = 6           ! dimension of points per xplasma call
    integer :: i,k,j                             ! loop variable and best point and temp
    integer :: jr(idsearch)                      ! region code at search points
 
    real(fp) :: zz                                   ! test point metric
    real(fp) :: zsa(idsearch), zsb(idsearch)         ! search points
    real(fp) :: zsxeq(idsearch), zsveq(idsearch)     ! x,v equations at search points
    real(fp) :: zsxvec(3,idsearch), zsvvec(3,idsearch) ! x,v vectors at each search point
 
    integer, parameter :: jmicro_max = 5         ! maximum number of micro iterations
    integer, parameter :: jhex = 6               ! number of points in each micro iteration <= idsearch
 
    real(fp), parameter :: zsin60   = 0.8660254037844386_fp  ! sin(60.)
    real(fp), parameter :: zxhex(6) = (/1.0_fp, 0.5_fp, -0.5_fp, -1.0_fp, -0.5_fp, 0.5_fp/)
    real(fp), parameter :: zyhex(6) = (/0.0_fp, zsin60, zsin60, 0.0_fp, -zsin60, -zsin60/)
 
    !
    ! -------------------- micro search ----------------------
    ! search around incrementally decreasing hexs -- not scaled with point
    ! zmicro is the radius of the search hex -- try expanding first before contracting
    !
    k = 0                                      ! best point
    zmicro = min(gflr_op(5),zmicro/zmicro_step) ! hex radius start bigger then before
    if (idiagnostic) print *, 'map_larmor_search: za=', za_arg1, ' metric=', zmetric
    do j = 1, jmicro_max                  ! loop over hex radii
       if (zmicro < 10._fp*gflr_op(3)) exit  ! it is hopeless
       if (idiagnostic) print *, 'map_larmor_search: zmicro = ', zmicro
       ! -- set up points --
       do i = 1, jhex
          zsa(i) = max(0._fp, za_arg1+zmicro*zxhex(i))  ! won't consider zsa=0. case below
          zsb(i) = zb_arg1+zmicro*zyhex(i)
       end do
 
       ! -- get x,v equations --
       call map_larmor_points(jhex, zsa, zsb, jr, zsxeq, zsveq, zsxvec, zsvvec)
 
       ! -- search for best point --
       do i = 1, jhex
          if (jr(i)<=3) then         ! data must be ok
             zz = (zsxeq(i)*zsxeq(i) + zsveq(i)*zsveq(i))/2._fp
             if (zz<0.95_fp*zmetric .and. zsa(i)>0._fp) then
                zmetric = zz         ! make this the new best point
                k = i
             end if
          end if
       end do
 
       if (k>0) then
          za_arg1 = zsa(k)            ! found it and done
          zb_arg1 = zsb(k)
          i_search_fail = .false.
          if (idiagnostic) print *, 'map_larmor_search: success -- metric=', zmetric
          return
       end if
 
       zmicro = zmicro*zmicro_step    ! reduce the search size
    end do    ! end micro loop
 
    i_search_fail = .true.    ! exceeded the iteration count
    if (idiagnostic) print *, 'map_larmor_search: failure'
 
    return
 
  end subroutine map_larmor_search
 
  !
  ! ============================== map_larmor_step =====================================
  ! given the (a,b) point, x,v function values and jacobian, this subroutine
  ! advances the point and returns the new location, function values, region code
  ! and vector position and velocity or returns a negative status for a failure
  !
  ! advancement follows the numerical recipes scheme for first trying
  ! the full newton step then backtracking to find a minimum of the
  ! metric (zx**2+zv**2)/2.
  !
  ! there is a step size maximum of the smaller of zstep_max and sqrt(a*a+b*b)
  !
  subroutine map_larmor_step(itest, za, zb, zx, zv, zjac, zconv_ab, zconv_xv, kmax, &
       i_status, jr, zxp, zvp, zstep_max)
 
    integer, intent(in) :: itest     ! 1-> do not iterate looking for a good point
                                     !     return error if full larmor step fails
    real(fp), intent(inout) :: za, zb    ! starting point and new point on return
    real(fp), intent(inout) :: zx, zv    ! initial x,v equations and new ones on return
    real(fp), intent(in) :: zjac(2,2)    ! jacobian at the initial point
 
    real(fp), intent(in)  :: zconv_ab    ! convergence criteria for a,b
    real(fp), intent(in)  :: zconv_xv    ! convergence criteria for x,v
    integer, intent(in) :: kmax      ! maximum number of iterations for minimizing metric
 
    integer, intent(out) :: i_status ! -1 -> failure to advance
                                     !  1 -> full newton advance
                                     !  n -> fractional advance with this many steps
    integer, intent(out) :: jr       ! final location region code
    real(fp), intent(out) :: zxp(3)      ! final location vector position
    real(fp), intent(out) :: zvp(3)      ! final location vector velocity
    real(fp), intent(in)  :: zstep_max   ! maximum step size allowed
 
    real(fp), parameter :: alpha = 1.e-4_fp ! test for rapid enough convergence of metric
 
    integer :: jsr(1)
    integer :: k                     ! loop variable
 
    real(fp) :: zz,zm,p2,s               ! temporary
    real(fp) :: pa, pb                   ! ideal newton step
    real(fp) :: fnew, fold               ! metrics for new and original point
    real(fp) :: gfa, gfb                 ! gradient of metric at original point
    real(fp) :: slope                    ! slope of metric along p direction
    real(fp) :: alambda                  ! scaled distance along p direction
    real(fp) :: zda, zdb                 ! change in a,b from old to new point
    real(fp) :: anew                     ! new alambda
    real(fp) :: r1,r2,d,u1,u2,a,b        ! temporary for finding minimum of cubic approximation
    real(fp) :: alambda2, fnew2          ! save guesses in iteration for cubic approximation
    real(fp) :: zsa(1), zsb(1), zsx(1), zsv(1), zsxp(3,1), zsvp(3,1)
    real(fp) :: zzj_inv(2,2)             ! inverse of the jacobian
 
    i_status = -1       ! pessimistic
    if (za<=0._fp) return  ! this would be really bad
 
    !
    ! --- get ideal step ---
    !
    zz = zjac(1,1)*zjac(2,2)-zjac(1,2)*zjac(2,1)   ! the determinant
    if (zz == 0) then
       if (idiagnostic) print *, 'map_larmor_step: singular jacobian'
       return
    end if
    zzj_inv(1,1) =  zjac(2,2)/zz   ! jacobian inverse, step size maximum handles
    zzj_inv(2,2) =  zjac(1,1)/zz   ! near singular jacobian
    zzj_inv(1,2) = -zjac(1,2)/zz
    zzj_inv(2,1) = -zjac(2,1)/zz
 
    pa = - (zx*zzj_inv(1,1)+zv*zzj_inv(1,2))   ! ideal step in a,b
    pb = - (zx*zzj_inv(2,1)+zv*zzj_inv(2,2))
    p2 = pa*pa+pb*pb
 
    zm = min(zstep_max*zstep_max, 0.01_fp+0.9_fp*(za*za+zb*zb)) ! maximum allowed step size
    if (p2>zm) then          ! shorten the step
       s = sqrt(zm/p2)
       pa = s*pa             ! normalize
       pb = s*pb
    end if
 
    ! --- get metric and grad(metric) ---
    fold = (zx*zx+zv*zv)/2._fp  ! f = metric at initial point
    gfa  = (zx*zjac(1,1)+zv*zjac(2,1))   ! grad[f] = f dot j
    gfb  = (zx*zjac(1,2)+zv*zjac(2,2))
 
    slope = gfa*pa+gfb*pb    ! slope of line f((a,b) + lambda*(pa,pb)) = grad[f] dot p
 
    alambda = min(1.0_fp, 0.9_fp*za/max(1.e-10_fp,-pa))  ! try full step first but make sure it does not take a<0.
    do k = 1, kmax
       ! --- advance ---
       zda = alambda*pa      ! step size in (a,b)
       zdb = alambda*pb
       if ((abs(zda)+abs(zdb)) < zconv_ab) then
          if (idiagnostic) print *, 'map_larmor_step: step size is too small'
          return
       end if
       zsa(1) = za + zda     ! finally the new point
       zsb(1) = zb + zdb
       if (zsa(1)<=0._fp) then  ! make sure does not happen even though bounded above
          if (idiagnostic) print *, 'map_larmor_step: step into a<0.'
          return
       end if
       call map_larmor_points(1, zsa, zsb, jsr, zsx, zsv, zsxp, zsvp)
 
       ! --- test for metric convergence ---
       fnew = (zsx(1)*zsx(1)+zsv(1)*zsv(1))/2._fp  ! new metric
       if (jsr(1)<4 .and. (fnew <= (fold + alpha*alambda*slope))) then
          za  = zsa(1)
          zb  = zsb(1)
          zx  = zsx(1)
          zv  = zsv(1)
          zxp = zsxp(:,1)
          zvp = zsvp(:,1)
          jr  = jsr(1)
          i_status = k
          return           ! success
       end if
 
       if (itest==1) return ! do not give a chance to improve
 
       ! -- fixup, attempt to find the minimum in this direction --
       if (k==1) then
          if (slope>=0._fp) then
             if (idiagnostic) print *, 'map_larmor_step: bad slope'
             return        ! should be decreasing
          end if
 
          anew = - (alambda*alambda)*slope / (2._fp*(fnew - fold - alambda*slope))  ! new estimate for lambda
       else
          ! --- estimate the minimum from a cubic approximation ---
          r1 = fnew  - fold - alambda*slope
          r2 = fnew2 - fold - alambda2*slope
          d = 1._fp/(alambda-alambda2)
          u1 = r1/alambda**2
          u2 = r2/alambda2**2
          a = d*(u1-u2)
          b = d*(-alambda2*u1+alambda*u2)
          if (a == 0._fp) then
             if (b == 0._fp) return     ! give up
             anew = -slope/(2._fp*b)
          else
             d = b*b-3._fp*a*slope
             if (d<0._fp) return        ! give up
             anew = (-b+sqrt(d))/(3._fp*a)
          end if
          anew = min(anew, 0.5_fp*alambda)  ! force to be less then half the previous estimate
       end if
       alambda2 = alambda                ! save previous lambda guess
       fnew2    = fnew                   ! save previous metric
       alambda  = max(anew, 0.1_fp*alambda) ! new guess at metric minimum -- don't guess too small
    end do
    return                   ! exceeded the maximum number of iterations
  end subroutine map_larmor_step
 
  !
  ! ============================== map_larmor_jacobian =====================================
  !
  subroutine map_larmor_jacobian(i_center, za_arg1, zb_arg1, za_fd, zb_fd, &
       jregion_arg1, zxeq_arg1, zveq_arg1, zxp_arg1, zvelocity_arg1, zjacobian_arg1)
    !
    ! calculates the x and v equations at a single point along with a finite
    ! difference approximation to the jacobian.  an error corresponds to jregion_arg1>=4
    !
    use nbi_com
 
    logical, intent(in) :: i_center       ! .true. when the center point data needs to be calculated
                                          ! .false. if the center point data has already been
                                          !         calculated and is in the input data -- just need
                                          !         the jacobian
 
    real(fp), intent(in) :: za_arg1           ! displacement orthogonal to b in larmor radii
    real(fp), intent(in) :: zb_arg1           ! displacement parallel to b in larmor radii
    real(fp), intent(in) :: za_fd, zb_fd      ! difference in a,b to use for finite difference
                                          ! approximation as requested by caller
 
    integer, intent(inout) :: jregion_arg1   ! xplasma return code for this point, the data is
                                             ! invalid if jregion_arg1>=4 for this point or if the
                                             ! jacobian could not be computed
    real(fp), intent(inout) :: zxeq_arg1         ! the value of the x equation at this point
    real(fp), intent(inout) :: zveq_arg1         ! the value of the v equation at this point
    real(fp), intent(inout) :: zxp_arg1(3)       ! position vector at this point
    real(fp), intent(inout) :: zvelocity_arg1(3) ! velocity vector at this point
    real(fp), intent(out) :: zjacobian_arg1(2,2) ! the finite difference approximation to the
                                             ! jacobian matrix d(x,v)/d(a,b)
 
    integer :: iia(2), iib(2)                ! indices of points for finite differencing
    integer :: jjr(5)                        ! region code at 5 points around the target
    real(fp)  :: za_diff, zb_diff                ! difference in a,b to use for finite difference
                                             ! approximation
    real(fp) :: zja(5), zjb(5), zjx(5), zjv(5)   ! values at 5 points around the target point
    real(fp) :: zjp(3,5)                         ! position vector
    real(fp) :: zjvec(3,5)                       ! velocity vector
 
    ! -- stabalize numerical finite difference --
    za_diff = max(zmin_scale_gc, za_fd)
    zb_diff = max(zmin_scale_gc, zb_fd)
 
    ! -- form a cross of points around a,b ---           ^
    zja(1) = za_arg1              !                      |
    zjb(1) = zb_arg1              !         5            b
    zja(2) = za_arg1 - za_diff    !         |              a - >
    zjb(2) = zb_arg1              !         |
    zja(3) = za_arg1 + za_diff    !  2 ---  1 --- 3
    zjb(3) = zb_arg1              !         |
    zja(4) = za_arg1              !         |
    zjb(4) = zb_arg1 - zb_diff    !         4
    zja(5) = za_arg1              !
    zjb(5) = zb_arg1 + zb_diff    !
 
    if (i_center) then
       call map_larmor_points(5, zja, zjb, jjr, zjx, zjv, zjp, zjvec)  ! get center point also
    else
       call map_larmor_points(4, zja(2:), zjb(2:), jjr(2:), zjx(2:), zjv(2:), zjp(1,2), zjvec(1,2))
       jjr(1) = jregion_arg1
       zjx(1) = zxeq_arg1
       zjv(1) = zveq_arg1      ! don't need point or velocity data
    end if
    if (jjr(1) >= 4) goto 10   ! -- oops, center point is out of there --
 
    ! --- determine which points to use for finite difference ---
    iia  = (/2,3/)    ! ideally use a central difference
    iib  = (/4,5/)
 
    if (jjr(2) >= 4) iia(1) = 1  ! throw out point 2, use center point instead
    if (jjr(3) >= 4) iia(2) = 1
    if (jjr(4) >= 4) iib(1) = 1
    if (jjr(5) >= 4) iib(2) = 1
 
    if ((iia(1) == iia(2)) .or. (iib(1) == iib(2))) goto 10  ! center point is ok but outside
                                                             ! points are not -- weird
 
    ! -- finite difference --
    zjacobian_arg1(1,1) = (zjx(iia(2)) - zjx(iia(1)))/(zja(iia(2)) - zja(iia(1)))
    zjacobian_arg1(1,2) = (zjx(iib(2)) - zjx(iib(1)))/(zjb(iib(2)) - zjb(iib(1)))
    zjacobian_arg1(2,1) = (zjv(iia(2)) - zjv(iia(1)))/(zja(iia(2)) - zja(iia(1)))
    zjacobian_arg1(2,2) = (zjv(iib(2)) - zjv(iib(1)))/(zjb(iib(2)) - zjb(iib(1)))
 
    if (i_center) then
       ! -- fill up the rest of the output --
       jregion_arg1 = jjr(1)
       zxeq_arg1    = zjx(1)
       zveq_arg1    = zjv(1)
 
       zvelocity_arg1(1) = zjvec(1,1)
       zvelocity_arg1(2) = zjvec(2,1)
       zvelocity_arg1(3) = zjvec(3,1)
 
       zxp_arg1(1) = zjp(1,1)
       zxp_arg1(2) = zjp(2,1)
       zxp_arg1(3) = zjp(3,1)
    end if
    return   ! normal exit
 
    !
    ! --- error exit ---
    !
10  continue
    jregion_arg1 = max(4, jjr(1))  ! at least a 4 for a bad exit
    zxeq_arg1 = 0._fp  ! be a good citizen and return values
    zveq_arg1 = 0._fp
    zjacobian_arg1 = 0._fp
    zxp_arg1(1) = 0._fp
    zxp_arg1(2) = 0._fp
    zxp_arg1(3) = 0._fp
    zvelocity_arg1(1) = 0._fp
    zvelocity_arg1(2) = 0._fp
    zvelocity_arg1(3) = 0._fp
 
  end subroutine map_larmor_jacobian
 
  !
  ! ============================== map_larmor_points ========================================
  !
  subroutine map_larmor_points(inp, za_arg, zb_arg, jregion_arg, zxeq_arg, zveq_arg, &
       zxp_arg, zvelocity_arg, ifail)
    !
    ! finds the x and v equations at a set of points in the
    ! plane of displacement
    !
    use nbi_com
    save
 
    integer, intent(in)  :: inp               ! number of points
    integer, intent(out) :: jregion_arg(*)    ! xplasma return code for each point
    real(fp), intent(in)  :: za_arg(*), zb_arg(*) ! displacement in larmor radii perpendicular
                                              ! and parallel to b
    real(fp), intent(out) :: zxeq_arg(*)          ! dimensionless  delta_x . b ---> should be zero
    real(fp), intent(out) :: zveq_arg(*)          ! dimensionless  v_from_deltax**2 - v_from_e_conservation**2
    real(fp), intent(out) :: zvelocity_arg(3,*)   ! the velocity vector
    real(fp), intent(out) :: zxp_arg(3,*)         ! x,y,z of real particle
 
    integer, optional, intent(out) :: ifail(*)  ! if present set this to 1 if the energy at the point is negative

    integer :: i                     ! loop over all the points
 
    real(fp) :: zdelta_x(3,ivecd)        ! displacement of real particle from guiding center
    real(fp) :: zr_p(ivecd)              ! r of real particle
    real(fp) :: zphi_hat(3,ivecd)        ! phi unit vector at real particle
    real(fp) :: zpsi_p                   ! poloidal flux in webers at real particle
    real(fp) :: zphi_p                   ! electric potential in volts at real particle
    real(fp) :: zbvec_p(3)               ! bx,by,bz at real particle
    real(fp) :: zb_p                     ! total b at real particle
    real(fp) :: zpar_p(3)                ! unit vector in direction of b at real particle
    real(fp) :: zomega_p                 ! gyrofrequency at real particle
    real(fp) :: zvphi_p                  ! velocity in toroidal direction at real particle
    real(fp) :: zvion2_p                 ! velocity**2 at real particle
    real(fp) :: zvion_p                  ! velocity at real particle
    real(fp) :: zdx_x_par(3)             ! delta_x x unit_parallel
    real(fp) :: zdx_d_par                ! delta_x . unit_parallel
    real(fp) :: zphi_d_par               ! unit_phi . unit_parallel
    real(fp) :: zdx_x_par_d_phi          ! (delta_x x unit_parallel) . unit_phi
    real(fp) :: zdx_x_par2               ! ||delta_x x unit_parallel||**2
    real(fp) :: zvpar                    ! parallel component of velocity
    real(fp) :: zvelocity(3)             ! total velocity vector
    real(fp) :: zz,yy                    ! temporary
 
    if (inp > ivecd) then
       write(nout,*) '? error in map_larmor, inconsistent size of inp,ivecd'
       call nbi_exit
    end if
 
    !
    ! --- load up the x,y,z coordinates of the points ---
    !
    do i = 1, inp
       yy = zscale_gc*za_arg(i)
       zz = zscale_gc*zb_arg(i)
       zdelta_x(1,i) = yy*zm_gc(1) + zz*zpar_gc(1)
       zdelta_x(2,i) = yy*zm_gc(2) + zz*zpar_gc(2)
       zdelta_x(3,i) = yy*zm_gc(3) + zz*zpar_gc(3)
 
       zxp_arg(1,i) = zx_gc(1) + zdelta_x(1,i)
       zxp_arg(2,i) = zx_gc(2) + zdelta_x(2,i)
       zxp_arg(3,i) = zx_gc(3) + zdelta_x(3,i)
 
       zr_p(i) = sqrt(zxp_arg(1,i)*zxp_arg(1,i)+zxp_arg(2,i)*zxp_arg(2,i))
 
       zphi_hat(1,i) = - zxp_arg(2,i)/zr_p(i)   ! zhat x (x,y,0) / sqrt(x**2+y**2)
       zphi_hat(2,i) =   zxp_arg(1,i)/zr_p(i)
       zphi_hat(3,i) = 0._fp
 
       zz_xx(i) = zxp_arg(1,i)/100._fp
       zz_yy(i) = zxp_arg(2,i)/100._fp
       zz_zz(i) = zxp_arg(3,i)/100._fp
    end do
 
    if (present(ifail)) ifail(1:inp)=0
 
   !
    ! --- get a vector of the b field and psi,phi functions ---
    !
    call map_larmor_bffget(inp, zz_xx, zz_yy, zz_zz, ztol, zbvec, ifcnd, ifcns, &
         ivecd, zfvals, jregion_arg, ierr)
 
    if (ierr /= 0) then
       write(nout,'(a)') '? error in call to eq_bffget in map_larmor'
       write(nout,'(a,e13.6)') '  gc xi  = ', xion
       write(nout,'(a,e13.6)') '  gc th  = ', (180._fp/xpi)*th
       write(nout,'(a,e13.6)') '  gc r   = ', rmjion
       write(nout,'(a,e13.6)') '  gc z   = ', yion
       call nbi_exit
    end if
 
    !
    ! --- calculate the x,v equations ---
    !
    do i = 1, inp
       if (jregion_arg(i)<=3) then
          zpsi_p = ip_dir*zconv_psi*zfvals(i,1)  ! in maxwells
          zphi_p = zfvals(i,2)                   ! keep in units of volts
 
          ! --- form parallel unit vector ---
          zbvec_p(1) = zbvec(1,i)
          zbvec_p(2) = zbvec(2,i)
          zbvec_p(3) = zbvec(3,i)
          zb_p = sqrt(zzdot(zbvec_p,zbvec_p))                ! total b field at p in teslas
          zpar_p(1) = zbvec_p(1)/zb_p                        ! unit vector in parallel direction
          zpar_p(2) = zbvec_p(2)/zb_p
          zpar_p(3) = zbvec_p(3)/zb_p
          zomega_p = (xzbeami*zz_e_c*1.e4_fp*zb_p)/rmassb       ! omega -- include tesla->gauss factor
 
          ! --- conservation with gc ---
          zvphi_p  = (z_p_phi_gc - zz_e_c*xzbeami*(zpsi_gc-zpsi_p))/(rmassb*zr_p(i)) ! conserved momentum
          zvion2_p = (z_e_gc + xzbeami*(zphi_gc-zphi_p))/v2toev                      ! conserved energy
          if (zvion2_p<0._fp) then
             if (present(ifail)) ifail(i)=1
             zvion2_p = 1._fp/v2toev   ! keep larger than 1 ev
          end if
          zvion_p  = sqrt(zvion2_p)
 
          zdx_x_par(1) = zdelta_x(2,i)*zpar_p(3) - zdelta_x(3,i)*zpar_p(2)
          zdx_x_par(2) = zdelta_x(3,i)*zpar_p(1) - zdelta_x(1,i)*zpar_p(3)
          zdx_x_par(3) = zdelta_x(1,i)*zpar_p(2) - zdelta_x(2,i)*zpar_p(1)
 
          zdx_d_par       = zzdot(zdelta_x(1,i),zpar_p)    ! will eventually be zero
          zphi_d_par      = zzdot(zphi_hat(1,i),zpar_p)    ! nonzero if there is a toroidal field
          zdx_x_par_d_phi = zzdot(zphi_hat(1,i),zdx_x_par)
          zdx_x_par2      = zzdot(zdx_x_par,zdx_x_par)
 
          zz    = zomega_p/zvion_p       ! dimensions of 1./distance otherwise boring
          zvpar = (zvphi_p - zomega_p*zdx_x_par_d_phi)/zphi_d_par
          zxeq_arg(i)  = zz*zdx_d_par
          zveq_arg(i)  = zz*zz*zdx_x_par2 + (zvpar/zvion_p)**2 - 1._fp
          zvelocity(1) = zvpar*zpar_p(1) + zomega_p*zdx_x_par(1)
          zvelocity(2) = zvpar*zpar_p(2) + zomega_p*zdx_x_par(2)
          zvelocity(3) = zvpar*zpar_p(3) + zomega_p*zdx_x_par(3)
 
          !write(65,'(15(g14.7,1x))') za_arg(i), zb_arg(i), zomega_p, zz*zz*zdx_x_par2/(za_arg(i)**2), &
          !     zz*zz*zdx_x_par2, (zvpar/zvion_p)**2, zvion_p, 1.-(zvpar/zvion_p)**2, &
          !     (vion**2-vpll**2)/zvion_p, zveq_arg(i), zvpar/zvion_p, zvphi_p/zvion_p, &
          !     zz*zdx_x_par_d_phi, 1./zphi_d_par, rmjion*vpll/(zr_p(i)*zvpar)
       else
          !
          ! --- point is outside the r,z grid ---
          !
          zxeq_arg(i) = sign(1._fp,zb_arg(i))    ! cleans up map_larmor_plot
          zveq_arg(i) = 1._fp  ! these should not be set to zero, see map_larmor_newton
          zvelocity   = 0._fp
       end if
       zvelocity_arg(1,i) = zvelocity(1)
       zvelocity_arg(2,i) = zvelocity(2)
       zvelocity_arg(3,i) = zvelocity(3)
       !write(81,'(10(es14.6,1x))') za_arg(i), zb_arg(i), zxeq_arg(i), zveq_arg(i), zvelocity_arg(:,i), zxp_arg(:,i)
    end do
    ! write(65,*) ' '
  end subroutine map_larmor_points
 
  !
  ! --- zzdot ---
  ! replaces sum(a(:)*b(:))
  !
  function zzdot(a,b) result(c)
    real(fp) :: a(*),b(*),c
    c = a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
  end function zzdot
 
  !
  ! ============================== map_larmor_bffget ========================================
  !
  subroutine map_larmor_bffget(ivec, zz_xx, zz_yy, zz_zz, ztol, zbvec, ilist, ifcns, &
       iivecd, zfvals, jregion, ierr)
 
    use nbi_com
 
    !
    ! this subroutine (almost) matches the functionality of xplasma eqxyz_bffget
    ! but runs faster by taking advantage of the r,z representation of the b field
    ! and functions.  the region return code is not the same as eqxyz_bffget.  the
    ! region code returned here is 1 if no errors and 4 if there was an error.
    ! the module option ibffget=.true. will force the xplasma version to be called
    ! in this subroutine.
    !
 
    integer, parameter   :: jfcnd = ifcnd+3  ! maximum number of functions and b field coordinates
 
    integer, intent(in)  :: ivec           ! number of coordinates
    integer, intent(in)  :: ilist          ! number of functions to evaluate
    integer, intent(in)  :: ifcns(ilist)   ! list of functions
    integer, intent(in)  :: iivecd         ! dimension of output arrays
    integer, intent(out) :: jregion(ivec)  ! crippled region code for each point
                                           ! 1-> point is inside r,z grid (unless ibffget=.true.)
                                           ! 4-> point is outside r,z grid
    integer, intent(out) :: ierr           ! nonzero on error
 
    real(fp), intent(in)  :: zz_xx(ivec), zz_yy(ivec), zz_zz(ivec)  ! coordinates in meters
    real(fp), intent(in)  :: ztol                   ! tolerance for xplasma
    real(fp), intent(out) :: zbvec(3,ivec)          ! x,y,z coordinates of b field
    real(fp), intent(out) :: zfvals(iivecd, ilist)  ! function values returned
 
    integer :: jf                          ! number of functions and b field coordinates
    integer :: jfcns(jfcnd)                ! list of functions and br,bz,bphi
    integer :: i,j                         ! temporary loop variables
    integer :: j_br, j_bz, j_bphi          ! index into function array for b field
    integer :: k                           ! index into function array for point
    integer :: jerr                        ! vector error return code
    integer :: kerr                        ! per point return code
 
    integer :: ilista                      ! adjusted list length (dmc nov 12 2001)
 
    real(fp) :: zr(ivecd)                      ! r coordinate of points
    real(fp) :: zans(ivecd,jfcnd)              ! result of xplasma call
    real(fp) :: zcos, zsin                     ! radial direction vector
 
    ! -------------------------------------------------------------------------------------------
    ! dmc:  trailer items in ifcns can be 0, i.e. omitted, return 0.0
 
    do i=ilist,1,-1
       if(ifcns(i).ne.0) then
          ilista=i
          exit
       endif
    enddo
 
    zfvals=0.0_fp
 
    if (ivec>ivecd .or. ilist>ifcnd) then
       call nbimsg_exit('?map_larmor_bffget: inconsistent size of ivec>ivecd or ilist>ifcnd')
    end if
 
    if (ibffget) then
       ! -- do it the slow way --
       call eqxyz_bffget(ivec, zz_xx, zz_yy, zz_zz, ztol, zbvec, ilista, ifcns, &
            iivecd, zfvals, jregion, ierr)
    else
       ! --- do it the fast way ---
       ! -- fill points --
       do i = 1, ivec
          zr(i) = sqrt(zz_xx(i)*zz_xx(i) + zz_yy(i)*zz_yy(i))
       end do
 
       ! -- fill functions --
       j_br   = ilista+1                  ! b field functions start after regular ones
       j_bz   = ilista+2
       j_bphi = ilista+3
 
       jfcns(1:ilista) = ifcns(1:ilista)   ! copy over standard functions
       jfcns(j_br)   = id_br             ! add in b fields to list of functions
       jfcns(j_bz)   = id_bz
       jfcns(j_bphi) = id_bphi
       jf = j_bphi
 
       ! -- call xplasma for all points --
       ! if this does not work, xplasma will be called for each
       ! point individually, this should rarely if ever happen
       !
       call eq_frz(ivec, zr, zz_zz, jf, jfcns, ivecd, zans, jerr)  ! interpolate in r,z plane
 
       do i = 1, ivec         ! loop over x,y,z points
          if (jerr==0) then
             k = i            ! index into function array for point i
             kerr = jerr
          else                ! vector call failed, evaluate each point individually
             call eq_frz(1, zr(i:i), zz_zz(i:i), jf, jfcns, ivecd, zans, kerr)  ! interpolate in r,z plane
             k = 1
          end if
 
          if (kerr==0) then
             ! -- this point is ok --
             jregion(i) = 1      ! could actually be 1,2,3
 
             ! -- set function values --
             do j=1,ilista
                zfvals(i,j) = zans(k,j)
             end do
 
             ! -- set b field direction --
             if (zr(i) == 0._fp) then
                zcos = 1._fp
                zsin = 0._fp
             else
                zcos = zz_xx(i)/zr(i)
                zsin = zz_yy(i)/zr(i)
             end if
             zbvec(1,i) = (zcos*zans(k,j_br) - zsin*zans(k,j_bphi))  ! r,z -> x,y,z
             zbvec(2,i) = (zsin*zans(k,j_br) + zcos*zans(k,j_bphi))
             zbvec(3,i) = zans(k,j_bz)
          else
             ! -- this is a bad point --
             jregion(i) = 4
             do j=1,ilista
                zfvals(i,j) = 0._fp
             end do
             zbvec(1,i) = 0._fp
             zbvec(2,i) = 0._fp
             zbvec(3,i) = 0._fp
          end if
       end do               ! end loop over points
       ierr = 0             ! let region code supply error information
    end if                  ! end if of using eqxyz_bffget
  end subroutine map_larmor_bffget
 
  !
  ! ============================== map_larmor_compare ========================================
  !
  subroutine map_larmor_compare(istep)
    !
    ! junk routine for comparing the  result of nbiflr with depsub
    ! when map_larmor was used
    ! the vpar,vperp returned by map_larmor is relative to the guiding center
    ! whereas depsub requires the velocity vector at the deposition point -- the
    ! variables zvr_ml, zvy_ml, zvtor_ml can be used to feed depsub.
    ! differences on the order of .2-.5 cm for a 70 kev particle have been observed
    ! due to differences in the magnetic field returned by nblocb in depsub.
    !
    use nbi_com
 
    integer, intent(in) :: istep          ! header integer
    real(fp) :: zz_rmjion, zz_yion, zz_vion   ! changed in depsub
    real(fp) :: zz_xksid, zz_xion, zz_th
    real(fp) :: zphi0, zphigc, ze0, zegc, zbt0, zbtgc
    real(fp) :: zb, zomega, zbpar(3), ziomega
    integer :: jregion(ivecd)
    real(fp) :: zvpar, zv2, zvperp, zrho
 
    if (.not. nlbgflr) then
       print *, '? map_larmor_compare: does not make sense to call this'
       print *, '  with nlbgflr=.false.'
    else
       zz_rmjion = rmjion
       zz_yion   = yion
       zz_vion   = vion
       zz_xksid  = xksid
       zz_xion   = xion
       zz_th     = th
 
       ifcns(1) = id_psi
       ifcns(2) = id_phi   ! electrostatic potential
 
       zz_xx(1) = rmjflr/100._fp
       zz_yy(1) = 0._fp
       zz_zz(1) = yflr/100._fp
       if(ifcns(2).ne.0) then
          call map_larmor_bffget(1, zz_xx, zz_yy, zz_zz, ztol, zbvec, ifcnd, ifcns, &
               ivecd, zfvals, jregion, ierr)
       else
          call map_larmor_bffget(1, zz_xx, zz_yy, zz_zz, ztol, zbvec, 1, ifcns(1:1), &
               ivecd, zfvals(:,1), jregion, ierr)
          zfvals(:,2)=0.0_fp
       endif

       if (ierr /= 0 .or. jregion(1)>3) then
          write(nout,'(a)') '? map_larmor_compare: error in call to eq_bffget'
          call nbi_exit
       end if
       zbvec(:,1) = 1.e4_fp*zbvec(:,1)             ! to gauss
       zb    = sqrt(sum(zbvec(:,1)*zbvec(:,1))) ! compare with zbmod in depsub
       zbpar = zbvec(:,1)/zb                    ! compare with zbdir in depsub
       zomega = (xzbeami*zz_e_c*zb)/rmassb
       ziomega = 1._fp/zomega                      ! compare with zrlov in depsub
 
       zvpar = zbpar(1)*zvr_ml+zbpar(2)*zvtor_ml+zbpar(3)*zvy_ml
       zv2   = zvr_ml**2+zvtor_ml**2+zvy_ml**2
       zvperp= sqrt(max(0._fp,zv2-zvpar*zvpar))
       zrho  = zvperp*ziomega
 
       call depsub(rmjflr, yflr, zvtor_ml, zvr_ml, zvy_ml, 2, &
            zphi0,zphigc,ze0,zegc,zbt0,zbtgc)
 
       print '(/a,i5,a)', '----------- ', istep, ' --------------'
       print '(a,f12.6)', 'rmjflr = ', rmjflr
       print '(a,f12.6)', 'yflr   = ', yflr
       print '(a,f12.6)', 'pho    = ', (180._fp/xpi)*atan2(sinpho,cospho)
       print '(a,f12.6,1x,f12.6)', 'rmjion = ', zz_rmjion, rmjion
       print '(a,f12.6,1x,f12.6)', 'yion   = ', zz_yion, yion
       print '(a,f12.6,1x,f12.6)', 'xksid  = ', zz_xksid, xksid
       print '(a,e17.6,1x,e17.6)', 'vion   = ', zz_vion, vion
 
       print '(a,f12.6,a)', 'b    = ', zb*1.e-4_fp, ' teslas'
       print '(a,f12.6,a)', 'rho  = ', zrho*1.e-2_fp, ' m'
       print '(a,f12.6)',   'vperp/v = ', zvperp/sqrt(zv2)
       print '(a,f12.6)',   'vpar/v =  ', zvpar/sqrt(zv2)
 
       rmjion = zz_rmjion
       yion   = zz_yion
       vion   = zz_vion
       xksid  = zz_xksid
       xion   = zz_xion
       th     = zz_th
    end if
  end subroutine map_larmor_compare
end module map_larmor_mod
! 10sep2004 fgtok -s r8_precision.sub nubeam_test_fp.sub "r8con.csh conversion"
! 10sep2004 fgtok
