module antenna
use precision
implicit none

private
public  antenna_initial, antenna_profile, &
        phiext, gradphiext, &
! test fields
        gradphiext2, phiext1, phiext2, g2phi, g2phi2, &
        phitest, phicoupjcore, phicoupjsol

real(lk), dimension(:,:), allocatable :: phiext
real(lk), dimension(:,:,:), allocatable :: gradphiext
! test fields
real(lk),dimension(:),allocatable ::  phitest,phicoupjcore,phicoupjsol
real(lk),dimension(:,:),allocatable :: phiext1,phiext2,g2phi,g2phi2
real(lk),dimension(:,:,:), allocatable :: gradphiext2

save

contains
!=================================================================
subroutine antenna_initial
    use grids_setup, only: mgrid, mpara_sep_sol, mpara_share
    use simulation_domain,only:couplefrc
    use global_parameters, only: mype, gtcout
    if(mype==0) then
        write(gtcout, *) "============================================="
        write(gtcout, *) "Antenna_initial Starts:"
        write(gtcout, *) "============================================="
        write(gtcout, *) "" !this is a new line
        flush(gtcout)
    endif

    allocate(phiext(0:1,mgrid), gradphiext(3,0:1,mgrid))

    phiext=0.0_lk
    gradphiext=0.0_lk

! allocate memory for test fields
    allocate(g2phi(0:1,mgrid),g2phi2(0:1,mgrid),gradphiext2(3,0:1,mgrid),&
             phiext1(0:1,mgrid),phiext2(0:1,mgrid))
  
    if(couplefrc==1)then
        allocate(phitest(mpara_sep_sol),phicoupjsol(mpara_share),phicoupjcore(mpara_share))
        phitest=0.0_lk
        phicoupjsol=0.0_lk
        phicoupjcore=0.0_lk
    endif

    phiext1=0.0_lk
    phiext2=0.0_lk
    g2phi=0.0_lk
    g2phi2=0.0_lk
    gradphiext2=0.0_lk
end subroutine antenna_initial
!========================================================================
subroutine antenna_profile
    use precision
    use global_parameters
    use particle_array
    use grids_fields
    use grids_setup
    use grids_eqdata
    use laplacian
    use read_in_parameters, only: irun, mstep, tstep, aelectron, qion, aion, &
                                  Zsol_domain0, Zsol_domain1, iflux_sol
    use simulation_domain, only: sfcore, sfsol, couplefrc, R1sol
    implicit none
  
    integer :: i,j,ij,n,m,ii
    real(lk) :: wt_phi,wa,zdum,envelope,temp,rdum,omega_temp_star,omega_ele_star,&
                rhos2, kperp2_rhos2, rhoi
    real(8) dbesjn,dbesyn,gam,gamr,emfactor,canon,kpara
  
    temp=(etemp0/eden0)*7.43e2*7.43e2/(r0*r0)
    emfactor=3.0*1e10*utime/ulength
    canon=rho0*rho0/temp/aelectron/emfactor/emfactor
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    phiext=0.0_lk
    n=nmode   ! zeta mode number  (determined by toroidaln in gtc.in)
    m=1   ! Z mode number
  
    ii=mgrid_core+igrid_sol(iflux_sol)+int(0.5_lk*real(mpara_sol(iflux_sol)))  ! middle radial point value
    gam=n/Rsol(ii-mgrid_core)
    gamr=pi/(rg1d_sol(mpsi_sol)-rg1d_sol(1))
  !  gam=sqrt(gam**2+gamr**2)
    kpara=2.0_lk*pi*m/totals_sol(mpsi_sol)
    omega_ele_star=-meshte(ii)/bmesh(ii)*kapaner(ii)*gam/qion
    omega_temp_star=meshti(ii)*kapatir(ii)*gam/qion/bmesh(ii)
    rhos2=meshte(ii)/bmesh(ii)/bmesh(ii)*(aion/qion/qion)
    rhoi=sqrt(meshti(ii)/bmesh(ii)/bmesh(ii))*(sqrt(aion)/qion)
    kperp2_rhos2=gam*gam*rhos2
    wa=0.5_lk*(kpara**2*rho0**2*omega_temp_star)**(1.0_lk/3.0_lk)
  !  wa=sqrt(meshne(ii)*1836.0*rho0*rho0/temp)  ! electron plasma frequency
  !  wa=1.414_lk*sqrt(meshni(ii)*rho0*rho0/temp)          ! ion plasma frequency
  !  wa=bmesh(ii)*(kpara/gam)*sqrt(aion/aelectron)      ! electron \omega_H mode frequency
  !  wa=kpara*rho0                   ! ion sound wave frequency
  !  wa=1.28_lk*kpara*rho0/sqrt(1.0_lk + kperp2_rhos2)                   ! ion sound wave frequency with polarization correction
  !  wa=1.033_lk
  !   wa=bmesh(ii)*(kpara/gam)                  ! ion omega-H mode frequency
  !   wa=-gam*meshte(ii)*kapanir(ii)/bmesh(ii) ! drift wave frequency
  !   wa=-gam*meshte(ii)*kapanir(ii)/bmesh(ii)/(1.0_lk + kperp2_rhos2) ! drift wave frequency
    if(mype==0 .and. istep==1 .and. irk==1)then
        write(gtcout,*)'Antenna freq=',wa,'=',wa/(pi2*utime),'Hz'
        call FLUSH(gtcout)
    endif
   
! phase of driving source
   ! wt_phi=cos(wa*tstep*(real(istep+mstep*irun)+0.5*irk)) !irun = # of restart runs
    wt_phi=wa*tstep*(real(istep+mstep*irun)+0.5*irk)
    !!!!!!!!!!!!!!!!!!!!!!!!!!   antenna  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if(sfcore==1)then
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                zdum=2.0*pi*real(j-1)*deltas_core1d(i)/totals_core(i)
                rdum=2.0_lk*pi*(psimesh_core(i)-psimesh_core(1))/(psimesh_core(mpsi_core)-psimesh_core(1))
                envelope=1.0_lk
                phiext(1,ij) = envelope*cos(1.0*rdum)*cos(real(m)*zdum)*cos(real(n)*zeta1)
            enddo
        enddo
    endif
  
    if(sfsol==1)then
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=igrid_sol(i)+j
                zdum=2.0*pi*Ssol(ij)
                rdum=2.0_lk*pi*(psimesh_sol(i)-psimesh_sol(1))/(psimesh_sol(mpsi_sol)-psimesh_sol(1))
       !     rdum=1.0_lk*pi*(rg1d_sol(i)-rg1d_sol(1))/(rg1d_sol(mpsi_sol)-rg1d_sol(1))
      !     rdum=0.0_lk
                envelope=1.0_lk
     !           phiext(1,mgrid_core+ij) = envelope*cos(real(m)*zdum)*cos(real(n)*zeta1)*cos(rdum)
                phiext(1,mgrid_core+ij) = envelope*cos(real(m)*zdum)*cos(real(n)*zeta1-wt_phi)!*sin(rdum)
       !     phiext(1,mgrid_core+ij) = envelope*cos(real(m)*zdum)*cos(-real(n)*zeta1-wt_phi)!*sin(rdum)
       !     phiext(1,mgrid_core+ij) = envelope*cos(real(m)*zdum)*wt_phi
       !     phiext(1,mgrid_core+ij) = envelope*cos(2.0*rdum)*cos(real(m)*zdum)!*wt_phi
        !    phiext(1,mgrid_core+ij)=meshni(mgrid_core+ij)
       !      phiext(1,mgrid_core+ij)=-(max(0.0_lk,Zsol(ij)-Zsol(igrid_sol(i)+0.7*mpara_sol(i))))**2
          enddo
        enddo
    endif
    
!call toroidal_filter(phiext)
    if(istep <300)then
!$omp parallel do private(i)
        do i=1,mgrid
    !   phi(1,i)=phi(1,i)+phiext(1,i)
        enddo
    endif
  
  
    if(sfcore==1)then
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                zdum=2.0_lk*pi*(Zcore(ij)-Zsol_domain0)/(Zsol_domain1-Zsol_domain0)
                rdum=2.0_lk*pi*Rcore(ij)/(R1sol)
                envelope=1.0e-2
                phiext1(1,ij) = envelope*cos(2.0_lk*rdum)*cos(2.0_lk*zdum)
            enddo
        enddo
    endif
  
    if(sfsol==1)then
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=igrid_sol(i)+j
                zdum=2.0_lk*pi*(Zsol(ij)-Zsol_domain0)/(Zsol_domain1-Zsol_domain0)
                rdum=2.0_lk*pi*Rsol(ij)/R1sol
                envelope=1.0e-2
                phiext1(1,mgrid_core+ij) = envelope*cos(2.0_lk*rdum)*cos(2.0_lk*zdum)
          enddo
        enddo
    endif
  
    call laplacian_perp(phiext1(1,:),g2phi(1,:))
    if(couplefrc==1)then
        call convert_irreg2reg(phiext1(1,:))
    endif
  
!=================================================================================
!                INTERPOLATION subroutine TEST
!=================================================================================
    if(couplefrc==1)then
        call convert_irreg2reg(phiext(1,:),phitest)
        phiext2=phiext
        call convert_reg2irreg(phiext2(1,:))
        phicoupjcore=0.0_lk
        phicoupjsol=0.0_lk
        do j=2,mpara_share-1
            phicoupjcore(j)=(1.0_lk-interp_weight_core_share(j))*phiext(1,igrid_core(mpsi_core-1)+interp_index_core_share(j))+interp_weight_core_share(j)*phiext(1,igrid_core(mpsi_core-1)+interp_index_core_share(j)+1)
        enddo
    
        do j=2,mpara_share-1
            phicoupjsol(j)=(1.0_lk-interp_weight_sol_share(j))*phiext(1,mgrid_core+igrid_sol(2)+interp_index_sol_share(j))+interp_weight_sol_share(j)*phiext(1,mgrid_core+igrid_sol(2)+interp_index_sol_share(j)+1)
        enddo
    endif
  
end subroutine antenna_profile
!=================================================================
end module antenna
