subroutine nbi_aimingangle(rmjflr,yflr,zangl,vpll, vperp,zvec,alpha_v)
  !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 
  !also return vr,vphi,vz velocity of the particle
  use nbi_com, only : xbjaco,nonlin,nsjdotb,id_bphi,id_br,id_bz,fp
  use xstraln_calls, only : xs_fxget
  use xplasma_calls
  use logmod
  implicit none
  real(fp),intent(in) ::  rmjflr,yflr !(r,z) ptcl. [cm,cm]
  real(fp),intent(in) ::  zangl !gyrophase angle [rad]
  real(fp),intent(in) ::  vpll,vperp !v_ll, v_perp [cm/s, cm/s]
  real(fp),intent(out) :: alpha_v ! poloidal aiming angle [rad]
  real(fp), intent(out) :: zvec(3)   ! vr,vphi,vz velocity of the particle [cm/s, rad/s, cm/s]
  real(fp) :: cospho,sinpho !cos(zangl), sin(zangl)
  real(fp) :: xiflr,thflr !(rho,theta) ptcl [cm,deg]
  real(fp) :: zphiarg !phi
  integer :: jfcns(3)       ! list of functions and br,bz,bphi
  real(fp) :: zz,zz2,zzjaco(2,2)
  logical :: ixmomry        ! .true. when jacobian is based on xmomry called within
  !  this subroutine
  logical :: im             ! .true. for nonorthogonality messages
  real(fp) :: zx_gc(3)        !  r,phi,z of guiding center
  real(fp) :: zbvec_gc(1,3)     !  guiding center br,bphi,bz
  !  unit vectors in cylindrical 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
  !  note:    e_g and e_ll are orthogonal vectors
  !
  !              ^ e_perp
  !              |   
  !              |     ^ alpha_gyro
  ! 	     |   |    /   
  !           |  |   /^  positive direction
  !            | |  /  ^    
  !             || /    ^
  !             ||/______^__ >e_g
  !             | ^
  !            |   \
  !           |     \
  !          |       e_ll = b/|b|
  
  real(fp) :: zxi_gc(3)       ! unit vector orthogonal to magnetic surface
  real(fp) :: zpar_gc(3)      ! unit vector in direction of b
  real(fp) :: zperp_gc(3)     ! unit vector on magnetic surface orthogonal to b
  real(fp) :: zm_gc(3)        ! direction of ideal larmor radius displacement
  real(fp) :: v_g,vg_perp     !  v = (v_g, v_ll, vg_perp) - velocity components in 
                            !  cylindrical coord.
  real(fp) :: zb_gc
  integer :: kerr,iwarn
!mg  integer :: id_br,id_bphi,id_bz
  
  cospho=cos(zangl)
  sinpho=sin(zangl)
  ixmomry = .false.         ! first try the original jacobian for xi vector
  
  zx_gc = (/rmjflr, 0._fp, yflr/)
!mg  call eq_gfnum('bphi', id_bphi)
!mg  call eq_gfnum('br',   id_br)
!mg  call eq_gfnum('bz',   id_bz)
  if(min(id_bphi,id_br,id_bz).eq.0) then
     write(nonlin,*) '??nbi_aimingangle: xplasma id missing for br bz bphi'
  endif
  jfcns(1)=id_br            ! add in b fields to list of functions
  jfcns(2)=id_bphi 
  jfcns(3)=id_bz

  zphiarg=0.0_fp
  call xs_fxget(0.01_fp,zphiarg,rmjflr,yflr,xiflr,thflr,iwarn)
  
  ! interpolate in r,z plane
  call eq_frz(1, zx_gc(1:1)/100._fp, zx_gc(3:3)/100._fp, 3, jfcns, 1, zbvec_gc(1:1,:), kerr)
  !     zbvec_gc - (br,bphi,bz)
  zb_gc    = sqrt(sum(zbvec_gc(1,:)*zbvec_gc(1,:))) ! total b field at gc in teslas
  if(zb_gc.gt.0.0_fp)then
     zpar_gc(:)  = zbvec_gc(1,:)/zb_gc ! unit vector in parallel direction
  else
     zpar_gc(:)  = zbvec_gc(1,:)
  endif
  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 = (xiflr<1.0_fp) .or. (abs(zz)>.2_fp) ! reduce number of messages
     if (im) then
        call infolog('nbi_aimingangle: b vector and xi vector not orthogonal')
        call infolog('nbi_aimingangle: angle(degrees) = ', 57.296_fp*asin(zz))
     end if
     if (.not. ixmomry) then
        ixmomry = .true.
        if(im) call infolog('nbi_aimingangle: trying to recover with xmomry..')
        call xmomry(xiflr, thflr, 1, zz, zz2, zzjaco) ! don't change original stuff
        goto 10
     end if
     if (im) then
        call infolog('nbi_aimingangle: enforcing orthogonality and continuing')
     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)
  
  ! find  v_g = (v e_g), vg_perp = (v e_perp)  componentes of velocity 
  
  ! form m vector, unit vector in ideal guiding center direction ---
  zm_gc = cospho*zxi_gc + sinpho*zperp_gc
  zvec = nsjdotb*vpll*zpar_gc + vperp*zm_gc
  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)
  
  return
  
contains
  
  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
  
end subroutine nbi_aimingangle
