! check the orbit crossing of last flux surface in (r,z) coordinate
subroutine orbrz_check_bdy(iorb, istep)

  use orbrzv
  use nbi_com,only:xbmbnd
  use xplasma_obj_instance

  implicit none
  integer, intent(in) :: iorb, istep ! orbit index and time step from orbit.for
  integer :: iwarn, ierr
  real(fp) :: zpsi

!mg  call xplasma_ctrans(s, iwarn, r_in = orbrzv_zz(1)*1.0e-2_fp, z_in = orbrzv_zz(2)*1.0e-2_fp, &
!mg                    & maptype=1, tol=1.0e-6_fp, i2pi=2, rho_out = orbrzv_rho, theta_out = orbrzv_th)

  call xs_fxget1_r8(0.01_fp, 0.0_fp, orbrzv_zz(1), orbrzv_zz(2), & 
       & orbrzv_rho, orbrzv_th, iwarn)

  !call ezspline_interp2(orbrzv_psi_spo,orbrzv_zz(1),orbrzv_zz(2),zpsi,ierr)
  !if(zpsi < orbrzv_psix) then
  !  call ezspline_interp1(orbrzv_rhot_spo, zpsi, orbrzv_rho, ierr)
  !endif

  if(orbrzv_zz(1) > orbrzv_rmin .and. orbrzv_zz(1) < orbrzv_rmax .and. &
   & orbrzv_zz(2) > orbrzv_zmin .and. orbrzv_zz(2) < orbrzv_zmax) then
    if(orbrzv_rho < xbmbnd) then
      orbrzv_pflag = 0
    else
      orbrzv_pflag = 1
    end if
  else
    orbrzv_pflag = -1
  end if

end subroutine orbrz_check_bdy


! restore initial status (i.e. before time advancement)
subroutine orbrz_backto0

  use orbrzv

  orbrzv_zz(1:4)    = orbrzv_zz0(1:4)
  orbrzv_dzzdt(1:3) = orbrzv_dzzdt0(1:3)
  orbrzv_psiv(1:6)  = orbrzv_psiv0(1:6)
  orbrzv_gv(1:3)    = orbrzv_gv0(1:3)
  orbrzv_phiv(1:3)  = orbrzv_phiv0(1:3)    
  orbrzv_b          = orbrzv_b0
  ! bugfix : jul.23.2007
  orbrzv_teng       = orbrzv_teng0
  orbrzv_tcmt       = orbrzv_tcmt0

end subroutine orbrz_backto0


! check rf resonance condition before(isub=0) and after(isub=1)
! the orbit time stepping
subroutine orbrz_check_rfres(iorb, istep)

  use orbrzv
  use nbi_com,only:bitime,twopi,goose,lsbeam,zc

  implicit none
  integer, intent(in) :: iorb, istep
  integer :: iant, inphi, il, ik
  real(fp) :: zqb0_mc, zqb_mc, zkl(1:2), zdsf0, zdsf, zlam0, zlam, &
            & ztime0, ztime, zttmp, zomg, zq_mc, zw0, zw1, zww

  if(orbrzv_pflag < 0) call bad_exit('?orbrz_chk_rfres : invalid pflag')

  orbrzv_rf_nkick = 0 ! count total number of kick
  orbrzv_rf_nbadc = 0 ! count total number of hfs -> lfs crossing

  ! gyro-frequency
  zq_mc = orbrzv_qs/(orbrzv_ms*zc)
  zqb0_mc = zq_mc*orbrzv_b0
  zqb_mc  = zq_mc*orbrzv_b

  ! exact time of rf-interaction is subtle.
  ! need to improve later...
  ztime  = bitime(1) + orbrzv_dt*goose
  ztime0 = orbrzv_time0 + orbrzv_dt

  ! after time stepping ----------------------------------------------------
  do iant = 1, orbrzv_rf_nant
    do inphi = 1, orbrzv_rf_nnphi(iant)
      if(orbrzv_rfw(iant,inphi)%freq.le.0.0_fp) cycle !mg, no data for this antenna
      zomg = orbrzv_rfw(iant,inphi)%freq
      zomg = twopi*zomg
      zkl(1) = orbrzv_rfw(iant,inphi)%kl*orbrzv_bpsign
      zkl(2) = -zkl(1)

      do ik = 1, orbrzv_rfw(iant,inphi)%ksym
        ! doffler shifted frequency
        zdsf0 = zomg - zkl(ik)*orbrzv_zz0(3)/orbrzv_zz0(1) 
        zdsf = zomg - zkl(ik)*orbrzv_zz(3)/orbrzv_zz(1)
          
        do il = orbrzv_lmin(lsbeam), orbrzv_lmax(lsbeam)
          ! resonance condition after the timestep
          zlam0 = zdsf0 - real(il,fp)*zqb0_mc !mg this is rf - resonance condition
          zlam = zdsf - real(il,fp)*zqb_mc


          ! check crossing
          if(zlam*zlam0 <= 0.0_fp) then

            ! increase total counter
            orbrzv_rf_nkick = orbrzv_rf_nkick + 1

            zw0 = abs(zlam)
            zw1 = abs(zlam0)
            zww = zw0 + zw1

            if(zww > 0.0_fp) then
              zw0 = zw0/zww
              orbrzv_rfw(iant,inphi)%lam0(il,ik) = zw0
              orbrzv_rfw(iant,inphi)%lam(il,ik) = 1.0_fp - zw0
            else
              orbrzv_rfw(iant,inphi)%lam0(il,ik) = 0.5_fp
              orbrzv_rfw(iant,inphi)%lam(il,ik) = 0.5_fp
            endif

            ! record time 
            zttmp = orbrzv_rfw(iant,inphi)%time(il,ik)
            orbrzv_rfw(iant,inphi)%dt(il,ik) = ztime - zttmp
            orbrzv_rfw(iant,inphi)%time(il,ik) = ztime

            ! record non-goosed time 
            zttmp = orbrzv_rfw(iant,inphi)%time0(il,ik)
            orbrzv_rfw(iant,inphi)%dt0(il,ik) = ztime0 - zttmp
            orbrzv_rfw(iant,inphi)%time0(il,ik) = ztime0

            ! check direction (bad corssing)
            if(zlam0 <= 0.0_fp) then
              orbrzv_rf_nbadc = orbrzv_rf_nbadc + 1
              ! bad crossing
              orbrzv_rfw(iant,inphi)%dsgn(il,ik) = -1
            else
              ! normal crossing
              orbrzv_rfw(iant,inphi)%dsgn(il,ik) = 1
            endif
          else
            ! no crossing
            orbrzv_rfw(iant,inphi)%dsgn(il,ik) = 0
            ! count number of steps between resonance crossing
            ! individually for each resonances
            orbrzv_rfw(iant,inphi)%nstep(il,ik) = 1.0_fp + &
               &  orbrzv_rfw(iant,inphi)%nstep(il,ik)
          endif

        enddo ! il
      enddo ! ik
    enddo ! inphi
  enddo ! iant

end subroutine orbrz_check_rfres

! calculate time step control factor for precise rf-resonance control
subroutine orbrz_rf_timft(ztf_min, zrfd_min, irff)

  use orbrzv
  use nbi_com,only:epslon,twopi,lsbeam,orbrzv_zzerr_con,zc

  implicit none
  integer, intent(in) :: irff
  real(fp), intent(out) :: ztf_min, zrfd_min

  real(fp) :: zq_mc, zqb_mc, zdr1, zdr2, zdsf, ztf, zrfd, zdlam, zdrd
  real(fp) :: zkl(1:2), zdvdr, zdkvdr, zomg
  integer :: iant, inphi, ik, il

  integer :: ierr
  real(fp) :: zlrf_inv

  ztf_min  = 1.0_fp
  zrfd_min = 1.0_fp

  zq_mc = orbrzv_qs/(orbrzv_ms*zc)
  zqb_mc = zq_mc*orbrzv_b

  ! dvll / dr
  if(orbrzv_zz(3) .ne. 0.0_fp) then
    zdvdr = orbrzv_zz(4)*orbrzv_dbdr + orbrzv_qs*orbrzv_phiv(1)
    zdvdr = -zdvdr/(orbrzv_ms*orbrzv_zz(3))
  else
    zdvdr = 0.0_fp
  endif

  ! estimated variation of r in the next time step
  zdr1 = orbrzv_tdt*orbrzv_dzzdt(1)
  zdr1 = zdr1 + sign(epslon, zdr1)

  ! referential distance fixed
  if(orbrzv_nmidc > 1) then
    zdrd = orbrzv_bnw
  else
    zdrd = 0.1_fp*(orbrzv_prmax - orbrzv_rc)
  endif

  do iant = 1, orbrzv_rf_nant
    do inphi = 1, orbrzv_rf_nnphi(iant)
      if(orbrzv_rfw(iant,inphi)%freq.le.0.0_fp) cycle !mg, no data for this antenna
      zomg = orbrzv_rfw(iant,inphi)%freq
      zomg = twopi*zomg
      zkl(1) = orbrzv_rfw(iant,inphi)%kl
      zkl(1) = zkl(1)/orbrzv_zz(1)
      zkl(2) = -zkl(1)

      do ik = 1, orbrzv_rfw(iant,inphi)%ksym
        ! doffler shifted frequency
        zdsf = zomg - zkl(ik)*orbrzv_zz(3)
        ! d(kll*vll)/dr
        zdkvdr = zkl(ik)*zdvdr - zkl(ik)/orbrzv_zz(1)

        do il = orbrzv_lmin(lsbeam), orbrzv_lmax(lsbeam)

          if(irff == 1 .and. orbrzv_rfw(iant,inphi)%dsgn(il,ik) < 0) then
          ! re-step for bad crossing
            ztf = 1.0_fp
            !ztf = orbrzv_rf_dist
          else
          ! normal time step
            zdlam = real(il,fp)*zq_mc*orbrzv_dbdr + zdkvdr
            zdlam = zdlam + dsign(epslon, zdlam)
            zdr2 = -(real(il,fp)*zqb_mc - zdsf)/zdlam

            ztf  = dabs(zdr2/zdr1)
            zrfd = dabs(zdr2/zdrd)
            !ztf = zdr2/zdr1
            !if(ztf < 0.0_fp) then
            !  ztf = orbrzv_rf_dist ! resonance is moving away !
            !else

            ! version-1
            !ztf = dmax1(1.0e-1_fp, dmin1(ztf, 1.0_fp))

            ! version-2
            !if(ztf >= 1.0_fp .and. ztf < orbrzv_rf_dist) then
            !  ztf = 1.0_fp - 1.0_fp/ztf
            !endif

            ! version-3
            !if(zdr2 <= 0.1_fp*(orbrzv_prmax-orbrzv_prmin)) then
            !  ztf = dmax1(0.1_fp*orbrzv_tdt, orbrzv_dt_min*100.0_fp)
            !else
            !  ztf = orbrzv_tdt
            !endif

            ! version-4
            !if(ztf < orbrzv_rf_dist) then
            !  ztf = min(0.2_fp, ztf)
            !endif

            ! version-5
            !ztf = ztf/orbrzv_rf_dist
            !ztf = max(ztf, 1.0/orbrzv_rf_dist)

            !ztf  = max(0.05_fp, 0.3_fp*ztf)
            !ztf  = 1.0_fp 
            ztf  = max(0.1_fp, 0.5_fp*ztf)
            zrfd = max(0.1_fp, zrfd)

          endif

          orbrzv_rfw(iant,inphi)%ztf(il,ik) = ztf
          ztf_min  = min(ztf, ztf_min)
          zrfd_min = min(zrfd, zrfd_min)
        enddo
      enddo
    enddo
  enddo

  orbrzv_ztf_min  = ztf_min
  orbrzv_zrfd_min = min(orbrzv_zrfd_min, zrfd_min)
  !orbrzv_zrfd_min = zrfd_min

  if(ztf_min < 1.0_fp) then
    orbrzv_zzerr_con0 = orbrzv_zzerr_con
    orbrzv_zzerr_con = min(orbrzv_zzerr_con,orbrzv_rferr_con)
  else
    orbrzv_zzerr_con = orbrzv_zzerr_con0
  endif

  ! keep track of maximum bnw/lrf along the orbit path
  if(zrfd_min < 1.0_fp) then
    call ezspline_interp2(orbrzv_lrfi_spo, orbrzv_zz(1), orbrzv_zz(2), &
         zlrf_inv, ierr)
    orbrzv_bnw_lrf = max(orbrzv_bnw_lrf, zlrf_inv*orbrzv_bnw)
  endif

end subroutine orbrz_rf_timft
