! quasi-linear heating operator version 1
! simpled goosed rf operator used
! with option to choose 
! 1. n independent q_rf
! 2. one q_rfg = sqrt(n)*q_rf
subroutine orbrz_qlrf(iorb, istep, ierr)

  use orbrzv
  use nbi_com, only:goose,epslon,orbrzv_toric_ufrnm,twopi,wght,orbrzv_prftot,nzones,&
       &xi,lcentr,lep1,lsbeam,orbrzv_prf,epsinv,ntest_max,zel,zmp,zc
  use nbi_random_pll
  use logmod
  !use xplasma_obj_instance

  implicit none
  intrinsic bessel_jn
  integer, intent(in) :: iorb, istep ! orbit index and time step from orbit.for
  integer, intent(out) :: ierr ! status code, 0 = ok

  integer :: iant, inphi, ik, il, inbz, idum, izg, itrap, &
           & iuncon, ixiav, ii, ixi, ixi1, idpt, idpts, iidp
  real(fp) :: zzm(1:3), zbm, zbdbm, ztr1, ztr2, ztr, &
            & zdomg, zddomg, zomgl, zdomgl, zddomgl
  real(fp) :: zvl, zvp2, zkl(1:2), zbb, zrho2, zomg, zdvl, zdvp2, zgdvp2
  real(fp) :: zll, zq_m, zq_mc, zq_c, zqe_m2, zrn1, zddt, zrfg, zrnum
  real(fp) :: zteng0, zdteng, zdwgt, zbfl, zrtmp, zztmp, znrfg, ztr_max, &
            & xjadep, zxa, zxksav, zeiav, zphiav, zxhw, zdx, zwt, &
            & zxdps, zzdek, zra
  real(fp) :: zvv0, zv0, zgvv, zgv, zgdv, zgdx, zdt, zdwdt,zxa_old
  integer :: ilb
  real(fp) :: zrn1_zrfg
  integer :: izrfg,i_zrfg

  ! initialize error flag
  ierr  = 0

  zq_m  = orbrzv_qs/orbrzv_ms
  zq_c  = orbrzv_qs/zc
  zq_mc = zq_m/zc

  ! trap or passing
  if(orbrzv_trpas > 0.0_fp) then
    zrnum = 4.0_fp ! trap
  else if(orbrzv_trpas < 0.0_fp) then
    zrnum = 2.0_fp ! passing
  else 
    zrnum = 0.0_fp ! undetermined
  endif

  ! goosed time (physical time!)
  zdt = orbrzv_dt*goose
  orbrzv_dbg_dt=0.0_fp
! mg i think it is supposed to be here
  ! if there is any single resonance crossing hfs -> lfs,
  ! take back to original position
  if(orbrzv_rf_nbadc > 0) call orbrz_backto0

 ! parallel & perp2 velocity
  zvl  = orbrzv_zz(3)
  zvp2 = 2.0_fp*orbrzv_b*orbrzv_zz(4)/orbrzv_ms

  ! jump to profile diagnostics if non-heating case
  if(orbrzv_rf_nkick .eq. 0) then
    zdteng = 0.0_fp
    orbrzv_where = 2
    goto 9999
  endif

  ! interpolate kp2, e+, pol
  call orbrz_intp_rf(iorb, istep, ierr)

!mg  ! if there is any single resonance crossing hfs -> lfs,
!mg  ! take back to original position
!mg  if(orbrzv_rf_nbadc > 0) call orbrz_backto0

  ! initial kinetic energy and velocity
  zteng0 = orbrzv_teng
  zvv0   = 2.0_fp*(orbrzv_teng-orbrzv_qs*orbrzv_phiv(3))/orbrzv_ms
  zv0    = dsqrt(zvv0)

!#ifdef __DEBUG
  ! debugging purpose initial call for indexing 
  orbrzv_where = 1
  zxa_old = xjadep(min(orbrzv_rho,1.0_fp),orbrzv_th,orbrzv_zz(1),orbrzv_zz(2), &
             & zv0, orbrzv_zz(3), &
             & itrap, iuncon, zxksav, zeiav, zphiav, ixiav, zxhw)
!#endif

  ! compute velocity change by rf...
  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, cycle if no data for this antenna
      zkl(1) = orbrzv_rfw(iant,inphi)%kl*orbrzv_bpsign
      zkl(2) = -zkl(1)

      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.and.il.ne.0) then
            
            zzm(1:3) = orbrzv_rfw(iant,inphi)%zzm(il,ik,1:3)
            zbm = orbrzv_rfw(iant,inphi)%zzm(il,ik,4)

            call ezspline_interp2(orbrzv_bdb_spo,zzm(1),zzm(2),zbdbm,ierr)
            zbdbm = abs(zbdbm) + epslon
            zbdbm = max(zbdbm, orbrzv_bdb_min)

            ! gyro-frequency and its derivatives
            zomg   = zq_mc*zbm
            zdomg  = zzm(3)*zbdbm*zq_mc
            zddomg = orbrzv_zz(4)*zq_c*(zbdbm/orbrzv_ms)**2
            zdomg  = abs(zdomg) + epslon
            zddomg = abs(zddomg) + epslon

            zrho2 = 2.0_fp*zbm*orbrzv_zz(4)/orbrzv_ms ! vp2m
            zrho2 = zrho2/(zomg*zomg) ! gyro-radius ** 2
  
            ! compute bessel function factor (fortran 2008 intrinsic)
            ! orbrzv_zbesl(l) = j_(l-1)
            zbb = orbrzv_rfw(iant,inphi)%kp2(il,ik)*zrho2
            zbb = sqrt(zbb)
            if(il.gt.0) then
               ilb=il
            else
               ilb=abs(il-1)+1 ! as orbrzv_zbesl(l) = j_(l-1)
            endif
            orbrzv_zbesl(1:ilb) = bessel_jn(0, ilb-1, zbb)

            zqe_m2 = zq_m*orbrzv_rfw(iant,inphi)%ep(il,ik)
            zqe_m2 = zqe_m2*orbrzv_toric_ufrnm(lsbeam)
            zqe_m2 = zqe_m2*zqe_m2

            zll = real(il, fp)
            ztr_max = orbrzv_trmax(lsbeam)/sqrt(abs(zll))
            zomgl   = zll*zomg
            zdomgl  = zll*zdomg
            zddomgl = zll*zddomg
            !mg 
            if(zvl .ne. 0.0_fp) then 
!mg, use this formula for nstxlike rf heating
!mg (many resonances)                
!mg               
               if(ntest_max.ne.0) then !option for nicola, nstx test  
                  ztr= orbrzv_b/abs(zzm(3)*zbdbm*zll) !mg 04/27/2017  see phys plasmas
                                                                        ! 12, 072524 (2005) page 7 on the paragraph between 
                                                                        ! formular (20) and (21)
!mg 04/27/2017               ztr = sqrt(twopi/abs(zdomg))
!mg    
               else
                  ztr = sqrt(twopi/abs(zdomgl)) 
                                        !mg the resonance duration time for particles, 
                                        !which have their turning points far away 
                                        !from cyclotron resonance 
                                        ! from eq(7) g.d. kerbel, m.g.mccoy 
                                        !computer phusics comm 40 (1986) 105-113
                                        ! also see
                                        !eq(19) kovanen,nucl.fusion(1992), vol.32, p787
                                        ! the difference of our formula from theirs is  (c/n)^(-1/2)
                                        ! where n -is harmonic number
                                        ! wave -particle interaction (resonance number)
               endif
            else 
               ztr = ztr_max
            end if
                                        ! in kovanen formulas eq(19)and eq(20) we can see (q/m) factor
                                        ! supposed to be (l*q/(m*c)) 
                                        ! common formula is given in eq(7) and (8) in 
                                        ! also eq (3) for  'nu' has to be taken into account
                                        ! g.d. kerbel and m.g. mccoy, 
                                        ! computer physics communications 40(1986)105-113
                                       

            if(abs(zvv0-zvp2).le.epslon) then 
                                        !mg the resonance duration time for particles, 
                                        !which have their turning points on the
                                        !cyclotron resonance layer
                                        ! from eq(8) g.d. kerbel, m.g.mccoy 
                                        !computer phusics comm 40 (1986) 105-113
                                        ! also see
                                        !eq(20) kovanen,nucl.fusion(1992), vol.32, p787
                                        ! the difference our formula from theirs is  (c/n)^(-1/3)
                                        ! where n -is harmonic number
                                        ! wave -particle interaction (resonance number)
               ztr2 = (0.5_fp*abs(zddomgl))**(1.0_fp/3.0_fp)
 !mg              ztr2 = (0.5_fp*abs(zddomg))**(1.0_fp/3.0_fp)
              ztr = twopi*0.355_fp/ztr2 !mg 0.355 came from airy function ai(0)=1/(3^(2/3)gamma(2/3))
            else
!              ztr = ztr_max
            end if
!mg  as the turning point approaches the cyclotron resonance and | v_|| | goes to zero
!    time ztr1 is increased towards infinity. for particles with turning points close 
!    to the resonance layer, the  variation of v_|| has to be taken into account (ztr2).
!            ztr = min(ztr1, ztr2)
            !ztr = min(ztr, ztr_max)
	    
            orbrzv_dbg_dt=orbrzv_dbg_dt+ztr
            zvv0 = zvl*zvl + zvp2
            zv0  = sqrt(zvv0)
            zrn1 = nbrandm_pll(idum) - 0.5_fp
            ! new rf goosing factor : may.29.2007
            zrfg = orbrzv_rfw(iant,inphi)%dt(il,ik)/orbrzv_rfw(iant,inphi)%dt0(il,ik)
            znrfg = orbrzv_dt*zrfg
            zddt = zqe_m2*orbrzv_zbesl(ilb)*orbrzv_zbesl(ilb)
            zddt = zddt*znrfg ! goosed rf : bugfix apr.27.2007
            zddt=1.0_fp/4.0_fp*ztr*zddt*(abs(zomgl)/orbrzv_rfw(iant,inphi)%freq/twopi)**2
            zgdvp2 =sqrt(24.0_fp*zvp2*abs(zddt)) !mg sqrt(12) =  2*sqrt(3)
!mg the expression below has been derived using eq(18) t.h.stix nucler fusion 15(1975), 737-754
!   and eqs.(a28), (a33-a35) x.q. xu, m.n. rosenbluth phys. fluids b, 3 (1991), 627-643
            zdvp2= zddt + zrn1*zgdvp2
            
            zdvp2=max(-zvp2,zdvp2) !mg 05/08/2013 to avoid zvp2<0 
            zvp2 = zvp2 + zdvp2   
         
            zvp2 = max(0.0_fp, zvp2)
            ! dvl = kl vp dvp / l omg_i    !mg formular (4.9) kennel, engelmann, phys.fluid (1966),vol.9,p2377
            ! vp dvp = 0.5 dvp2
            zdvl = zkl(ik)*0.5_fp/(zomgl*zzm(1))
!!!mg11.16 
!!!mg11.16            zvl = zvl + zddt*zdvl + zrn1*zgdvp2*abs(zdvl)
            zvl = zvl + sign(1.0_fp,orbrzv_zz(3))*zddt*abs(zdvl) + zrn1*zgdvp2*abs(zdvl)
            zgvv = abs(zvp2 - zdvp2 + zddt + 0.5_fp*zgdvp2 + zvl*zvl)
            zgv  = sqrt(zgvv)

            zgdv  = abs(zgv - zv0)/zv0
            zgdx  = abs(zgdvp2)/zvv0

            ! record rf-induced changes for goose update
            orbrzv_rfdv(0) = orbrzv_rfdv(0) + 1.0_fp
            orbrzv_rfdv(1) = orbrzv_rfdv(1) + zgdv
            orbrzv_rfdv(2) = orbrzv_rfdv(2) + zgdx

!#ifdef __DEBUG
            ! monitor resoance crossing ------------------------------------
            rfdbg_pdiag(1,ixiav) = rfdbg_pdiag(1,ixiav) + 1.0_fp
            rfdbg_pdiag(2,ixiav) = rfdbg_pdiag(2,ixiav) + zddt
            rfdbg_pdiag(3,ixiav) = rfdbg_pdiag(3,ixiav) + zrfg
            rfdbg_pdiag(4,ixiav) = rfdbg_pdiag(4,ixiav) + zgdv
            rfdbg_pdiag(5,ixiav) = rfdbg_pdiag(5,ixiav) + zgdx
            rfdbg_pdiag(6,ixiav) = rfdbg_pdiag(6,ixiav) + goose
            rfdbg_pdiag(7,ixiav) = rfdbg_pdiag(7,ixiav) + zdt*zrnum
            rfdbg_pdiag(8,ixiav) = rfdbg_pdiag(8,ixiav) + orbrzv_bnw_lrf
            rfdbg_pdiag(9,ixiav) = rfdbg_pdiag(9,ixiav) + orbrzv_bnw

            orbrzv_rfw(iant,inphi)%nstep(il,ik) = 0.0_fp
            ! monitor resoance crossing ------------------------------------
!#endif


          endif
        enddo
      enddo

    enddo
  enddo

  orbrzv_zz(3) = zvl
  orbrzv_zz(4) = (0.5_fp*orbrzv_ms*zvp2)/orbrzv_b

  ! re-evaluate time derivative for next step
  call orbrz_derivs(4, orbrzv_zz, orbrzv_ms, orbrzv_qs, orbrzv_dzzdt)

  ! energy change by rf
  zdteng = orbrzv_teng - zteng0

  ! multiply particle weight
  zdteng = zdteng*wght

  ! recalculate time integration
  if((orbrzv_rf_nbadc > 0)) then  
    ! ensure to pass though the bad resonance in this re-stepping
     orbrzv_tdt = orbrzv_dt
     call orbrz4(iorb, istep, 1, ierr) ! rk4
    !call orbrz23(iorb, istep, 1, ierr) ! rk2 & rk3

    ! after orbit integration
     if(orbrzv_pflag.ne.0) then
        call warnlog(' orbrz_qlrf: orbit lost after re-advance in qlrf')
        goto 9999
     endif
  endif

  ! add rf-power deposition
  orbrzv_prftot(lsbeam) = orbrzv_prftot(lsbeam) + zdteng



  orbrzv_where = 2
  zxa = xjadep(min(orbrzv_rho,1.0_fp), orbrzv_th, &
             & orbrzv_zz(1),orbrzv_zz(2),dsqrt(zvl*zvl+zvp2),zvl,&
             & itrap, iuncon, zxksav, zeiav, zphiav, ixiav, zxhw)
  !if(iuncon == 0) then
  !  orbrzv_prf(lsbeam,ixiav) = orbrzv_prf(lsbeam,ixiav) + zdteng
  !endif

  !zwt = (zxa-xi(lcentr,1))/zdx
  !ixi = idnint(zwt) + lcentr; ixi = max(ixi, lcentr)
  !ixi1 = ixi + 1; ixi1 = min(ixi1, lep1)
  !zwt = zwt - real(ixi-lcentr,fp)
  !if(iuncon == 0) then
  !  orbrzv_prf(lsbeam,ixi)  = orbrzv_prf(lsbeam,ixi)  + (1.0_fp-zwt)*zdteng
  !  orbrzv_prf(lsbeam,ixi1) = orbrzv_prf(lsbeam,ixi1) + zwt*zdteng
  !endif

  if(iuncon.eq.0) then
    zdx = 1.0_fp/real(nzones,fp)
    if(zxa.eq.epsinv) then !if so it will case a crash
       call warnlog(' orbrz_qlrf: could be so called bethatron particles')
       call warnlog(' orbrz_qlrf: zxa=',zxa_old )
       if(zxa_old.eq.epsinv) then 
          zwt=0.0_fp
          ixi=lcentr
          ixi1 = ixi + 1
       else
          zxa=zxa_old !mg 06/03/2013
          !could be so called bethatron particles - passing 
          !with  vperp = 0
          !they are orbiting near magnetic axes, never crosses midplane
          zwt = (zxa-xi(lcentr,1))/zdx
          ixi = int(zwt) + lcentr; ixi = max(ixi, lcentr)
          ixi1 = ixi + 1; ixi1 = min(ixi1, lep1)
          zwt = zwt - real(ixi-lcentr,fp)
          if(zwt.gt.1.0_fp.or.zwt.lt.0.0_fp) then
             call errorlog(' orbrz_qlrf: zwt not in [0,1], zwt=',zwt)
          endif
       endif
       orbrzv_prf(ixi,lsbeam)  = orbrzv_prf(ixi,lsbeam)  + (1.0_fp-zwt)*zdteng
       orbrzv_prf(ixi1,lsbeam) = orbrzv_prf(ixi1,lsbeam) + zwt*zdteng
    else

       idpt = idnint(zxhw/zdx) + 1
       idpts = 2*idpt - 1
       zzdek = zdteng/real(idpts, fp)
       zdwgt = wght/real(idpts, fp)
       zdwdt = zdwgt*zdt
       
       do iidp = 1, idpts 
          zxdps = zxa + real(iidp-idpt,fp)*zdx
          zwt = (zxdps-xi(lcentr,1))/zdx
          ixi  = int(zwt) + lcentr
          ixi1 = ixi + 1
          ixi  = min(max(lcentr,ixi),lep1)
          ixi1 = min(max(lcentr,ixi1),lep1)
          zwt = zwt - real(ixi-lcentr,fp)
          zwt=max(0.0_fp,zwt)
          zwt=min(1.0_fp,zwt)
          if(zwt.gt.1.0_fp.or.zwt.lt.0.0_fp) then
             call errorlog(' orbrz_qlrf: zwt not in  the [0,1], zwt=',zwt)
          endif
          orbrzv_prf(ixi,lsbeam)  = orbrzv_prf(ixi,lsbeam)  + (1.0_fp-zwt)*zzdek
          orbrzv_prf(ixi1,lsbeam) = orbrzv_prf(ixi1,lsbeam) + zwt*zzdek
          
!#ifdef __DEBUG
          orbrzv_nrf(lsbeam,1,ixi)=orbrzv_nrf(lsbeam,1,ixi)+(1.0_fp-zwt)*zdwdt
          orbrzv_nrf(lsbeam,1,ixi1)=orbrzv_nrf(lsbeam,1,ixi1)+zwt*zdwdt
          
          if(orbrzv_trpas > 0.0_fp) then
             orbrzv_nrf(lsbeam,2,ixi)=orbrzv_nrf(lsbeam,2,ixi)+(1.0_fp-zwt)*zdwdt
             orbrzv_nrf(lsbeam,2,ixi1)=orbrzv_nrf(lsbeam,2,ixi1)+zwt*zdwdt
          else if(orbrzv_trpas < 0.0_fp) then
             orbrzv_nrf(lsbeam,3,ixi)=orbrzv_nrf(lsbeam,3,ixi)+(1.0_fp-zwt)*zdwdt
             orbrzv_nrf(lsbeam,3,ixi1)=orbrzv_nrf(lsbeam,3,ixi1)+zwt*zdwdt
          endif
!#endif

       enddo
    endif
 endif

9999 continue   ! jump to here for lost and for non-heating case
end subroutine orbrz_qlrf


