subroutine orbrz_intp_rf(iorb, istep, ierr)

  use ezspline, only : ezspline_interp2
  use orbrzv
  use nbi_com,only:twopi,lsbeam

  implicit none

  integer, intent(in) :: iorb, istep
  integer, intent(out) :: ierr  ! status code on exit (0=ok)

  integer :: iant, inphi, ik, il
  real(fp) :: repol, impol, zrm, zzm, zxm, ztm, zw0, zw1, zrtmp, zztmp

  ! ---------------------------------------------------------------------
  ! initialize e+, e-, kp for each antenna
  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

      do ik = 1, orbrzv_rfw(iant,inphi)%ksym
        do il = orbrzv_lmin(lsbeam), orbrzv_lmax(lsbeam)

          if(orbrzv_rfw(iant,inphi)%dsgn(il,ik) .ne. 0) then

            zw0 = orbrzv_rfw(iant,inphi)%lam0(il,ik)
            zw1 = orbrzv_rfw(iant,inphi)%lam(il,ik)

            ! linear interpolate the resonance point
            zxm = zw0*orbrzv_psiv0(1) + zw1*orbrzv_psiv(1)

            !if(iorb > 980) then
            !  print *, 'orbrz_intp_rf : ', iorb, istep
            !  print *, '              : ', zw0, zw1, zxm
            !  print *, '              : ', orbrzv_psiv0(1), orbrzv_psiv(1)
            !endif

            zxm = abs(zxm)
            zxm = sqrt(zxm/orbrzv_psix)
            zxm = max(0.0_fp, min(zxm, 1.0_fp))

            !zrm = orbrzv_zz(1)
            !zzm = orbrzv_zz(2)
            !orbrzv_rfw(iant,inphi)%zzm(il,ik,1) = zrm
            !orbrzv_rfw(iant,inphi)%zzm(il,ik,2) = zzm
            zrm = zw0*orbrzv_zz0(1) + zw1*orbrzv_zz(1)
            zzm = zw0*orbrzv_zz0(2) + zw1*orbrzv_zz(2)
            orbrzv_rfw(iant,inphi)%zzm(il,ik,1) = zrm
            orbrzv_rfw(iant,inphi)%zzm(il,ik,2) = zzm

            !orbrzv_rfw(iant,inphi)%zzm(il,ik,3) = orbrzv_zz(3)
            !orbrzv_rfw(iant,inphi)%zzm(il,ik,4) = orbrzv_b
            orbrzv_rfw(iant,inphi)%zzm(il,ik,3) = &
              & zw0*orbrzv_zz0(3) + zw1*orbrzv_zz(3)
            orbrzv_rfw(iant,inphi)%zzm(il,ik,4) = &
              & zw0*orbrzv_b0 + zw1*orbrzv_b

            zrtmp = zrm - orbrzv_rc
            zztmp = zzm - orbrzv_zc

            if(zrtmp == 0.0_fp .and. zztmp == 0.0_fp) then
              ztm = 0.0_fp
            else
              ztm = atan2(zztmp, zrtmp)
            endif
            if(ztm < 0.0_fp) ztm = ztm + twopi

            ! k_perp^2
            call ezspline_interp2(orbrzv_rfw(iant,inphi)%kp_spo, &
                    & zxm, ztm, orbrzv_rfw(iant,inphi)%kp2(il,ik), ierr)
            if(orbrzv_rfw(iant,inphi)%kp2(il,ik).lt.0.0_fp)then
               orbrzv_rfw(iant,inphi)%kp2(il,ik)=0.0_fp
            endif
            !orbrzv_zkp(iant, inphi) = dabs(orbrzv_zkp(iant, inphi))

            ! e+
            call ezspline_interp2(orbrzv_rfw(iant,inphi)%ep_spo, &
                    & zxm, ztm, orbrzv_rfw(iant,inphi)%ep(il,ik), ierr)
            if(orbrzv_rfw(iant,inphi)%ep(il,ik).lt.0.0_fp)then
               orbrzv_rfw(iant,inphi)%ep(il,ik)=0.0_fp
            endif

            ! e-/e+
            call ezspline_interp2(orbrzv_rfw(iant,inphi)%rp_spo, &
                    & zxm, ztm, repol, ierr)
            call ezspline_interp2(orbrzv_rfw(iant,inphi)%ip_spo, &
                    & zxm, ztm, impol, ierr)

            orbrzv_rfw(iant,inphi)%cpol(il,ik) = dcmplx(repol, impol)
          endif

        enddo
      enddo
    enddo
  enddo

end subroutine orbrz_intp_rf
