gray/src/beams.f90

module beams
  use const_and_precisions, only : wp_, one

  implicit none

contains

  subroutine read_beam0(params, err)
    ! Reads the wave launcher parameters for the simple case
    ! where w(z) and 1/R(z) are fixed.
    !
    ! The file should be formatted as follows:
    !
    !   1 f
    !   2 x₀ y₀ z₀
    !   3 w₀₁ w₀₂ d₀₁ d₀₂ φ
    !
    ! where:
    !   - f is the frequency (GHz)
    !   - x₀, y₀, z₀ are the launcher position (cm)
    !   - w₀₁,w₀₂ are the beam waists in the two principal directions (cm)
    !   - d₀₁,d₀₂ are the distances of the beam waists from the launch
    !     point (cm)
    !   - φ is the rotation angle from the horizontal direction to the
    !     first principal direction (deg)
    !
    ! Note: this case implies simple astigmatism, i.e. the
    ! amplitude and phase ellipses in the transverse plane
    ! are aligned.

    use const_and_precisions, only : pi, vc=>ccgs_
    use gray_params,          only : antenna_parameters
    use logger,               only : log_error

    ! subroutine arguments
    type(antenna_parameters), intent(inout) :: params
    integer,                  intent(out)   :: err

    ! local variables
    integer :: u
    real(wp_) :: k0, w0(2), d0(2), z_R(2), phi

    open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
    if (err /= 0) then
      call log_error('opening beams file ('//trim(params%filenm)//') failed!', &
                     mod='beams', proc="read_beam0")
      return
    end if

    read(u, *) params%fghz  ! Wave frequency (GHz)
    read(u, *) params%pos   ! Launcher position (Cartesian coordinates)
    read(u, *) w0, &        ! Beam waists (in the ξ,η axes)
               d0, &        ! Waists (along ξ,η axes) distance from launcher
               phi          ! rotation angle from horizontal direction to
                            ! ξ axis in the transverse plane
    close(u)

    ! Wavevector k₀
    k0 = 2.0e9_wp_* pi * params%fghz / vc

    ! Rayleigh ranges in the ξ, η directions:
    ! z_R = ½k₀w₀²
    z_R = k0/2 * w0**2

    ! Beam widths in the ξ, η directions at z=0:
    ! w² = w₀² (1 + z₀²/z_R²)
    params%w = w0 * sqrt(1 + (d0/z_R)**2)

    ! Curvature in the ξ, η directions at z=0:
    ! K = 1/R = (z-z₀)/[(z-z₀)² + z_R²]
    params%ri = -d0 / (d0**2 + z_R**2)

    ! Ellipse axes angle
    params%phi = phi
  end subroutine read_beam0


  subroutine read_beam1(params, err)
    ! Reads the wave launcher parameters for the case
    ! where w(z, α) and 1/R(z, α) depend on the launcher angle α.
    !
    ! Format notes:
    !   1. The first two lines contain the frequency in GHz
    !      and the number of rows.
    !   2. The rest of file is a table with the following columns:
    !        θ, α, β, x₀, y₀, z₀, w₁, w₂, k₁, k₂, φ_w, φ_R
    !   3 The meaning of the columns is
    !     - θ is the mechanical steering angle (unused)
    !     - α, β are the poloidal and toroidal launch angles (deg)
    !     - x₀, y₀, z₀ are the launcher position (mm)
    !     - w₁,w₂ are the beam widths in the two principal directions (mm)
    !     - k₁,k₂ are the wavefront curvatures in the two principal
    !       directions (mm⁻¹)
    !     - φ_w, φ_R are the rotation angles of the amplitude and phase
    !       ellipses in the transverse plane at the launch point (deg)

    use gray_params,  only : antenna_parameters
    use splines,      only : spline_simple
    use utils,        only : locate
    use logger,       only : log_error

    ! subroutine arguments
    type(antenna_parameters), intent(inout) :: params
    integer,                  intent(out)   :: err

    ! local variables
    integer :: u, nisteer, i, k, ii
    real(wp_) :: steer, dal
    real(wp_), dimension(:), allocatable ::                &
      alphastv, betastv, x00v, y00v,                       &
      z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v
    type(spline_simple) :: beta, waist1, waist2,   &
                           rci1, rci2, phi1, phi2, &
                           x0, y0, z0

    open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
    if (err /= 0) then
      call log_error('opening beams file ('//trim(params%filenm)//') failed!', &
                     mod='beams', proc="read_beam1")
      return
    end if
    read(u,*) params%fghz
    read(u,*) nisteer

    allocate(alphastv(nisteer), betastv(nisteer), waist1v(nisteer), &
             waist2v(nisteer), rci1v(nisteer), rci2v(nisteer),      &
             phi1v(nisteer), phi2v(nisteer), x00v(nisteer),         &
             y00v(nisteer), z00v(nisteer))

    do i=1,nisteer
      read(u, *) steer, alphastv(i), betastv(i), &
                 x00v(i), y00v(i), z00v(i),      &
                 waist1v(i), waist2v(i),         &
                 rci1v(i), rci2v(i),             &
                 phi1v(i), phi2v(i)
    end do
    close(u)

    ! initial beam data measured in mm -> transformed to cm
    x00v    = 0.1_wp_ * x00v
    y00v    = 0.1_wp_ * y00v
    z00v    = 0.1_wp_ * z00v
    waist1v = 0.1_wp_ * waist1v
    waist2v = 0.1_wp_ * waist2v
    rci1v   = 10 * rci1v
    rci2v   = 10 * rci2v

    call beta%init(alphastv, betastv, nisteer)
    call waist1%init(alphastv, waist1v, nisteer)
    call waist2%init(alphastv, waist2v, nisteer)
    call rci1%init(alphastv, rci1v, nisteer)
    call rci2%init(alphastv, rci2v, nisteer)
    call phi1%init(alphastv, phi1v, nisteer)
    call phi2%init(alphastv, phi2v, nisteer)
    call x0%init(alphastv, x00v, nisteer)
    call y0%init(alphastv, y00v, nisteer)
    call z0%init(alphastv, z00v, nisteer)

    if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then
      call locate(alphastv, params%alpha, k)
      dal           = params%alpha - alphastv(k)
      params%beta   = beta%raw_eval(k, dal)
      params%pos(1) = x0%raw_eval(k, dal)
      params%pos(2) = y0%raw_eval(k, dal)
      params%pos(3) = z0%raw_eval(k, dal)
      params%w(1)   = waist1%raw_eval(k, dal)
      params%w(2)   = waist2%raw_eval(k, dal)
      params%ri(1)  = rci1%raw_eval(k, dal)
      params%ri(2)  = rci2%raw_eval(k, dal)
      params%phi(1) = phi1%raw_eval(k, dal)
      params%phi(2) = phi2%raw_eval(k, dal)
    else
      ! params%alpha outside table range
      if(params%alpha >= alphastv(nisteer)) ii=nisteer
      if(params%alpha <= alphastv(1)) ii=1
      params%alpha  = alphastv(ii)
      params%beta   = betastv(ii)
      params%pos(1) = x00v(ii)
      params%pos(2) = y00v(ii)
      params%pos(3) = z00v(ii)
      params%w(1)   = waist1v(ii)
      params%w(2)   = waist2v(ii)
      params%ri(1)  = rci1v(ii)
      params%ri(2)  = rci2v(ii)
      params%phi(1) = phi1v(ii)
      params%phi(2) = phi2v(ii)
    end if

    deallocate(alphastv, betastv, waist1v, waist2v, &
               rci1v, rci2v, phi1v, phi2v, &
               x00v, y00v, z00v)

    call beta%deinit
    call waist1%deinit
    call waist2%deinit
    call rci1%deinit
    call rci2%deinit
    call phi1%deinit
    call phi2%deinit
    call x0%deinit
    call y0%deinit
    call z0%deinit

  end subroutine read_beam1


  subroutine read_beam2(params, beamid, err, set_pol)
    ! Reads the wave launcher parameters for the general case
    ! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β.
    !
    ! Format notes:
    !   1. The first line specifies the number N of beams described by the file.
    !   2. The rest of the files consists of N 2D tables preceded by a header
    !        id, mode, f, na, nb
    !      where
    !       - id is a string identifier of the beam
    !       - mode indicates where the beam has O-mode (1) or X-mode (2)
    !         polarisation
    !       - f is the frequency (GHz)
    !       - nα, nβ are the numbers of rows and columns of the table.
    !   3. The 2D table is stored in row-major order over *nα×nβ* lines, that is, the
    !      i,j-th record is stored on the l-th line, with l = i + nα*j.
    !   4. The poloidal angle *α(i,j)* must be monotonic along *i* and the toroidal angle
    !      *β(i, j)* must be monotonic along *j*.
    !   5. Each line stores one record with the same fields as in the 1D format.

    use gray_params, only : antenna_parameters
    use utils,       only : linear_interp, locate
    use types,       only : contour
    use dierckx,     only : curfit, splev, surfit, bispev
    use logger,      only : log_error

    ! subroutine arguments
    type(antenna_parameters), intent(inout)         :: params      ! beam parameters
    integer,                  intent(in)            :: beamid      ! which beam to load
    integer,                  intent(out), optional :: err
    logical,                  intent(in),  optional :: set_pol     ! whether to set params%iox

    ! local variables
    character(len=20) :: beamname
    logical :: set_pol_
    integer :: u
    integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta, iox
    integer :: iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk
    integer :: nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2
    integer :: nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0
    integer :: nxz0, nyz0, kx, ky, ii, npolyg, nmax, lwrk2, in
    integer :: nxycoord
    integer, DIMENSION(:), ALLOCATABLE :: iwrk
    real(wp_) :: alphast,betast, waist1, waist2, rci1, rci2, phi1, phi2
    real(wp_) ::  fp, minx, maxx, miny, maxy, eps, xcoord0, ycoord0
    real(wp_), DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv, &
       betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord,   &
       ycoord, wrk, txwaist1, tywaist1, txwaist2, tywaist2,    &
       txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2,  &
       txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1,  &
       cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2
    type(contour) :: polyg, polygA, polygB, polygC, polygD, outA, outB, outC
    real(wp_), DIMENSION(4) :: xvert, yvert
    real(wp_), dimension(1) :: fi
    integer, parameter :: kspl=1
    real(wp_), parameter :: sspl=0.01_wp_

    set_pol_ = .true.
    if (present(set_pol)) set_pol_ = set_pol

    open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
    if (err /= 0) then
      call log_error('opening beams file ('//trim(params%filenm)//') failed!', &
                     mod='beams', proc="read_beam1")
      return
    end if

!=======================================================================================
!     # of beams        
        read(u,*) nbeam
!       
!     unused beams' data
        jumprow=0
!  c====================================================================================
        do i=1,beamid-1
          read(u,*) beamname, iox, params%fghz, nalpha, nbeta
          jumprow = jumprow+nalpha*nbeta
        end do
!  c====================================================================================
!       
!     beam of interest
        read(u,*) beamname, iox, params%fghz, nalpha, nbeta
        if (set_pol_) params%iox = iox
!       
!  c====================================================================================
!     unused beams' data grids
        do i=1,(nbeam - beamid)
          read(u,*) beamname
        end do
        do i=1,jumprow
          read(u,*) alphast,betast,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2
        end do
!  c====================================================================================
!       
!      # of elements in beam data grid
        nisteer = nalpha*nbeta
!       
        allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer),   &
         waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), &
         phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer),      &
         xcoord(nisteer),ycoord(nisteer))
!       
!  c====================================================================================
!     beam data grid reading
        do i=1,nisteer
          read(u,*) alphast,betast,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2
!       
!     initial beam data (params%pos(1), params%pos(2), params%pos(3)) are measured in mm -> transformed to cm
          x00v(i)=0.1d0*params%pos(1)
          y00v(i)=0.1d0*params%pos(2)
          z00v(i)=0.1d0*params%pos(3)
          alphastv(i)=alphast
          betastv(i)=betast
          waist1v(i)=0.1d0*waist1
          rci1v(i)=1.0d1*rci1
          waist2v(i)=0.1d0*waist2
          rci2v(i)=1.0d1*rci2
          phi1v(i)=phi1
          phi2v(i)=phi2
        end do
    close(u)
!  c====================================================================================
!       
!     fdeg = 0 alpha, beta free variables
!            1 alpha free variable
!            2 beta free variable
!            3 no free variables
        fdeg = 2*(1/nalpha) + 1/nbeta

!#######################################################################################
!     
!     no free variables
      if(fdeg.eq.3) then
        params%alpha=alphastv(1)
        params%beta=betastv(1)
        params%pos(1)=x00v(1)
        params%pos(2)=y00v(1)
        params%pos(3)=z00v(1)
        params%w(1)=waist1v(1)
        params%w(2)=waist2v(1)
        params%ri(1)=rci1v(1)
        params%ri(2)=rci2v(1)
        params%phi(2)=phi1v(1)
        params%phi(1)=phi2v(1)
        deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, &
        phi2v,x00v,y00v,z00v,xcoord,ycoord)
        return
      end if
!#######################################################################################
!
!
!#######################################################################################
      if(fdeg.eq.2) then
!     beta = independent variable
!     alpha = dependent variable
        xcoord = betastv
        ycoord = alphastv
        xcoord0 = params%beta
        ycoord0 = params%alpha
        kx=min(nbeta-1,kspl)
!  c====================================================================================
      else
!  c====================================================================================
!     alpha = independent variable
!     beta = dependent/independent (fdeg = 1/0)
        xcoord = alphastv
        ycoord = betastv
        xcoord0 = params%alpha
        ycoord0 = params%beta
        nxcoord = nalpha
        nycoord = nbeta
        kx=min(nalpha-1,kspl)
        ky=min(nbeta-1,kspl)
      end if
!#######################################################################################
!
      iopt = 0
      incheck = 0
!     
!#######################################################################################
      if(fdeg.ne.0) then
        nxest = kx + nxcoord + 1
        lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx))
        kwrk = nxest
        allocate(cycoord(nxest), txycoord(nxest), cwaist1(nxest),    &
        txwaist1(nxest), cwaist2(nxest), txwaist2(nxest),            &
        crci1(nxest), txrci1(nxest), crci2(nxest), txrci2(nxest),    &
        cphi1(nxest), txphi1(nxest), cphi2(nxest), txphi2(nxest),    &
        cx0(nxest), txx0(nxest), cy0(nxest), txy0(nxest),            &
        cz0(nxest), txz0(nxest), w(nxcoord), wrk(lwrk), iwrk(kwrk))
!       
        w = 1.d0
!       
!     2D interpolation
        call curfit(iopt,nxcoord,xcoord,ycoord,w,  &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord,  &
        txycoord,cycoord,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,waist1v,w, &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1,  &
        txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,waist2v,w, &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2,  &
        txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,rci1v,w,   &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1,    &
        txrci1,crci1,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,rci2v,w,   &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci2,    &
        txrci2,crci2,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,phi1v,w,   &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1,    &
        txphi1,cphi1,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,phi2v,w,   &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2,    &
        txphi2,cphi2,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,x00v,w,    &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0,      &
        txx0,cx0,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,y00v,w,    &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0,      &
        txy0,cy0,fp,wrk,lwrk,iwrk,ier)
!      
        call curfit(iopt,nxcoord,xcoord,z00v,w,    &
        xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0,      &
        txz0,cz0,fp,wrk,lwrk,iwrk,ier)
!     
!     check if xcoord0 is out of table range
!       incheck = 1 inside / 0 outside
        if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer)) incheck=1
!  c====================================================================================
      else
!  c====================================================================================
        npolyg = 2*(nxcoord+nycoord-2)
        minx = minval(xcoord)
        maxx = maxval(xcoord)
        miny = minval(ycoord)
        maxy = maxval(ycoord)
        nxest = 2*(kx + 1)
        nyest = 2*(ky + 1)
        nmax = max(nxest,nyest)+max(kx,ky)
        eps = 10.**(-8)
        lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19
        lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2
        kwrk = nisteer+(nmax-3)*(nmax-3)
        allocate(cwaist1(nxest*nyest), txwaist1(nmax), tywaist1(nmax), &
        cwaist2(nxest*nyest), txwaist2(nmax), tywaist2(nmax),          &
        crci1(nxest*nyest), txrci1(nmax), tyrci1(nmax),                &
        crci2(nxest*nyest), txrci2(nmax), tyrci2(nmax),                &
        cphi1(nxest*nyest), txphi1(nmax), typhi1(nmax),                &
        cphi2(nxest*nyest), txphi2(nmax), typhi2(nmax),                &
        cx0(nxest*nyest), txx0(nmax), tyx0(nmax),                      &
        cy0(nxest*nyest), txy0(nmax), tyy0(nmax),                      &
        cz0(nxest*nyest), txz0(nmax), tyz0(nmax),                      &
        wrk(lwrk), wrk2(lwrk2), iwrk(kwrk),                            &
        polyg%R(npolyg), polyg%z(npolyg), w(nisteer))
!       
        w = 1.d0
!       
!     3D interpolation
        call surfit(iopt,nisteer,xcoord,ycoord,waist1v,w,      &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxwaist1,txwaist1,nywaist1,tywaist1,cwaist1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
        call surfit(iopt,nisteer,xcoord,ycoord,waist2v,w,      &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxwaist2,txwaist2,nywaist2,tywaist2,cwaist2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
        call surfit(iopt,nisteer,xcoord,ycoord,rci1v,w,        &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxrci1,txrci1,nyrci1,tyrci1,crci1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
        call surfit(iopt,nisteer,xcoord,ycoord,rci2v,w,        &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxrci2,txrci2,nyrci2,tyrci2,crci2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
        call surfit(iopt,nisteer,xcoord,ycoord,phi1v,w,        &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxphi1,txphi1,nyphi1,typhi1,cphi1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
        call surfit(iopt,nisteer,xcoord,ycoord,phi2v,w,        &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxphi2,txphi2,nyphi2,typhi2,cphi2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
        call surfit(iopt,nisteer,xcoord,ycoord,x00v,w,         &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxx0,txx0,nyx0,tyx0,cx0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!    
       call surfit(iopt,nisteer,xcoord,ycoord,y00v,w,          &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxy0,txy0,nyy0,tyy0,cy0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!     
        call surfit(iopt,nisteer,xcoord,ycoord,z00v,w,         &
        minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,           &
        nxz0,txz0,nyz0,tyz0,cz0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!     data range polygon
        polyg%R(1:nxcoord) = xcoord(1:nxcoord)
        polyg%z(1:nxcoord) = ycoord(1:nxcoord)
!        
!  c====================================================================================
        do i=1,nycoord-2
          polyg%R(nxcoord+i) = xcoord((i+1)*nxcoord)
          polyg%R(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1)
          polyg%z(nxcoord+i) = ycoord((i+1)*nxcoord)
          polyg%z(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1)
        end do
!  c====================================================================================
        do i=1,nxcoord
          polyg%R(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1)
          polyg%z(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1)
        end do
!  c====================================================================================
!        
!     check if (xcoord0, ycoord0) is out of table range
!       incheck = 1 inside / 0 outside
        if (polyg%contains(xcoord0, ycoord0)) incheck = 1
      end if
      deallocate(wrk,iwrk)
!#######################################################################################
!     
!     
!#######################################################################################
      if(fdeg.ne.0) then
!  c====================================================================================
        if(incheck.eq.1) then
          call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier)
          ycoord0=fi(1)
          call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier)
          params%w(1)=fi(1)
          call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier)
          params%w(2)=fi(1)
          call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier)
          params%ri(1)=fi(1)
          call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier)
          params%ri(2)=fi(1)
          call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier)
          params%phi(2)=fi(1)
          call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier)
          params%phi(1)=fi(1)
          call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier)
          params%pos(1)=fi(1)
          call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier)
          params%pos(2)=fi(1)
          call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier)
          params%pos(3)=fi(1)
!    c----------------------------------------------------------------------------------
        else
!    c----------------------------------------------------------------------------------
          if(xcoord0.ge.xcoord(nisteer)) ii=nisteer
          if(xcoord0.le.xcoord(1)) ii=1
!         
          xcoord0=xcoord(ii)
          ycoord0=ycoord(ii)
          params%pos(1)=x00v(ii)
          params%pos(2)=y00v(ii)
          params%pos(3)=z00v(ii)
          params%w(1)=waist1v(ii)
          params%w(2)=waist2v(ii)
          params%ri(1)=rci1v(ii)
          params%ri(2)=rci2v(ii)
          params%phi(2)=phi1v(ii)
          params%phi(1)=phi2v(ii)
        end if
!  c====================================================================================
      else
!  c====================================================================================
        if(incheck.eq.0) then
          allocate(polygA%R(nxcoord), polygA%z(nxcoord))
          allocate(polygC%R(nxcoord), polygC%z(nxcoord))
          allocate(polygB%R(nycoord), polygB%z(nycoord))
          allocate(polygD%R(nycoord), polygD%z(nycoord))
          allocate(outA%R(nxcoord+3), outA%z(nxcoord+3))
          allocate(outB%R(nycoord+3), outB%z(nycoord+3))
          allocate(outC%R(nxcoord+3), outC%z(nxcoord+3))
!     coordinates of vertices v1,v2,v3,v4
          xvert(1) = polyg%R(1)
          xvert(2) = polyg%R(nxcoord)
          xvert(3) = polyg%R(nxcoord+nycoord-1)
          xvert(4) = polyg%R(2*nxcoord+nycoord-2)
          yvert(1) = polyg%z(1)
          yvert(2) = polyg%z(nxcoord)
          yvert(3) = polyg%z(nxcoord+nycoord-1)
          yvert(4) = polyg%z(2*nxcoord+nycoord-2)
!     coordinates of side A,B,C,D 
          polygA%R = polyg%R(1:nxcoord)
          polygA%z = polyg%z(1:nxcoord)
          polygB%R = polyg%R(nxcoord:nxcoord+nycoord-1)
          polygB%z = polyg%z(nxcoord:nxcoord+nycoord-1)
          polygC%R = polyg%R(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
          polygC%z = polyg%z(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
          polygD%R(1:nycoord-1) = polyg%R(2*nxcoord+nycoord-2:npolyg)
          polygD%R(nycoord) = polyg%R(1)
          polygD%z(1:nycoord-1) = polyg%z(2*nxcoord+nycoord-2:npolyg)
          polygD%z(nycoord) = polyg%z(1)
!     contour of outside regions A (1,2), B(3,4), C(5,6)
          outA%R = huge(one)
          outA%R(1:nxcoord) = polygA%R
          outA%R(nxcoord+3) = xvert(1)
          outA%z = -huge(one)
          outA%z(1:nxcoord) = polygA%z
          outA%z(nxcoord+1) = yvert(2)
          outB%R = huge(one)
          outB%R(1:nycoord) = polygB%R
          outB%R(nycoord+1) = xvert(3)
          outB%z = huge(one)
          outB%z(1:nycoord) = polygB%z
          outB%z(nycoord+3) = yvert(2)
          outC%R = -huge(one)
          outC%R(1:nxcoord) = polygC%R
          outC%R(nxcoord+3) = xvert(3)
          outC%z = huge(one)
          outC%z(1:nxcoord) = polygC%z
          outC%z(nxcoord+1) = yvert(4)
!     c----------------------------------------------------------------------------------
!      search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid
!      
!      (6)   |        (5)       |        (4)    
!      _ _ _v4__________________v3 _ _ _ _ _  
!             \        C        \
!              \                 \              (1)->(8) outside regions
!      (7)    D \                 \ B    (3)    v1->v4 grid vertices
!                \                 \            A-D grid sides   
!      _ _ _ _ _ _\_________________\_ _ _ _
!                v1         A       v2         
!      (8)        |      (1)        |    (2)     
!     
          if (outA%contains(xcoord0, ycoord0)) then
            in = 1
            if(xcoord0.GT.xvert(2)) then
              in = 2
            end if
          else if (outB%contains(xcoord0, ycoord0)) then
            in = 3
            if(ycoord0.GT.yvert(3)) then
              in = 4
            end if
          else if (outC%contains(xcoord0, ycoord0)) then
            in = 5
            if(xcoord0.LT.xvert(4)) then
              in = 6
            end if
          else
            in = 7
            if(ycoord0.LT.yvert(1)) then
              in = 8
            end if
          end if
!    c----------------------------------------------------------------------------------
!     (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border
!     depending on the region
!       1: xcoord0 unchanged, ycoord0 moved on side A
!       3: xcoord0 moved on side B, ycoord0 unchanged
!       5: xcoord0 unchanged, ycoord0 moved on side C
!       7: xcoord0 moved on side D, ycoord0 unchanged
!       2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates
!     in 1,3,5,7 incheck is set back to 1 to evaluate params%pos(1),params%pos(2),params%pos(3),waist,rci,phi in
!     new (xcoord0,ycoord0)
!     in 2,4,6,8 incheck remains 0 and params%pos(1),params%pos(2),params%pos(3),waist,rci,phi values at the
!     (xcoord0,ycoord0) vertex are used
          params%alpha = xcoord0
          params%beta = ycoord0
          SELECT CASE (in)
            CASE (1)
              !    params%beta outside table range
!             locate position of xcoord0 with respect to x coordinates of side A
              call locate(polygA%R, xcoord0, ii)
!             find corresponding y value on side A for xcoord position
              ycoord0 = linear_interp(polygA%R(ii), polygA%z(ii), polygA%R(ii+1), polygA%z(ii+1), xcoord0)
              incheck = 1
            CASE (2)
              !    params%alpha and params%beta outside table range
!             xcoord0, ycoord0 set 
              xcoord0 = xvert(2)
              ycoord0 = yvert(2)
              ii = nxcoord !indice per assegnare valori waist, rci, phi
            CASE (3)
              !    params%alpha outside table range
              call locate(polygB%z, ycoord0, ii)
              xcoord0 = linear_interp(polygB%z(ii), polygB%R(ii), polygB%z(ii+1), polygB%R(ii+1), ycoord0)
              incheck = 1
            CASE (4)
              !    params%alpha and params%beta outside table range
              xcoord0 = xvert(3)
              ycoord0 = yvert(3)
              ii = nxcoord+nycoord-1
            CASE (5)
              !    params%beta outside table range
              call locate(polygC%R, xcoord0, ii)
              ycoord0 = linear_interp(polygC%R(ii+1), polygC%z(ii+1), polygC%R(ii), polygC%z(ii), xcoord0)
              incheck = 1
            CASE (6)
              !    params%alpha and params%beta outside table range
              xcoord0 = xvert(4)
              ycoord0 = yvert(4)
              ii = 2*nxcoord+nycoord-2
            CASE (7)
              !    params%alpha outside table range
              call locate(polygD%z, ycoord0, ii)
              xcoord0 = linear_interp(polygD%z(ii), polygD%R(ii), polygD%z(ii+1), polygD%R(ii+1), ycoord0)
              incheck = 1
            CASE (8)
              !    params%alpha and params%beta outside table range
              xcoord0 = xvert(1)
              ycoord0 = yvert(1)
              ii = 1
          END SELECT
!    c----------------------------------------------------------------------------------
!     
          deallocate(polygA%R, polygA%z, polygC%R, polygC%z, polygB%R, polygB%z, &
                     polygD%R, polygD%z, outA%R, outA%z, outB%R, outB%z, outC%R, outC%z)
        end if
!  c====================================================================================  
!
!  c====================================================================================
        if(incheck.eq.1) then
          lwrk = 2*(kx+ky+2)
          kwrk = 4
          allocate(wrk(lwrk),iwrk(kwrk))
          call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%w(1)=fi(1)
          call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%w(2)=fi(1)
          call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1,           &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%ri(1)=fi(1)
          call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2,           &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%ri(2)=fi(1)
          call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1,           &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%phi(2)=fi(1)
          call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2,           &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%phi(1)=fi(1)
          call bispev(txx0,nxx0,tyx0,nyx0,cx0,                     &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%pos(1)=fi(1)
          call bispev(txy0,nxy0,tyy0,nyy0,cy0,                     &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%pos(2)=fi(1)
          call bispev(txz0,nxz0,tyz0,nyz0,cz0,                     &
            kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
          params%pos(3)=fi(1)
          deallocate(wrk,iwrk)
!    c----------------------------------------------------------------------------------
        else
!    c----------------------------------------------------------------------------------
          params%pos(1)=x00v(ii)
          params%pos(2)=y00v(ii)
          params%pos(3)=z00v(ii)
          params%w(1)=waist1v(ii)
          params%w(2)=waist2v(ii)
          params%ri(1)=rci1v(ii)
          params%ri(2)=rci2v(ii)
          params%phi(2)=phi1v(ii)
          params%phi(1)=phi2v(ii)
        end if
!  c====================================================================================
      end if
!#######################################################################################
!
      if(fdeg.ne.0) then
        deallocate(cycoord, txycoord, cwaist1, txwaist1, cwaist2, &
        txwaist2,  crci1, txrci1, crci2, txrci2, cphi1, txphi1,   &
        cphi2, txphi2, cx0, txx0, cy0, txy0, cz0, txz0, w)
      else
        deallocate(cwaist1, txwaist1, tywaist1, cwaist2, txwaist2, tywaist2, &
        crci1, txrci1, tyrci1, crci2, txrci2, tyrci2,       &
        cphi1, txphi1, typhi1, cphi2, txphi2, typhi2,       &
        cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0,  &
        wrk2, polyg%R, polyg%z, w)
      end if
!
!#######################################################################################
!     set correct values for alpha, beta
      if(fdeg.eq.2) then
        params%alpha = ycoord0
        params%beta = xcoord0
      else
        params%alpha = xcoord0
        params%beta = ycoord0
      end if
!#######################################################################################
      deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, &
        phi2v,x00v,y00v,z00v,xcoord,ycoord)
!
  end subroutine read_beam2


  subroutine launchangles2n(params, N)
    ! Given the wave launcher `params` computes the initial
    ! wavevector N̅ = ck̅/ω.
    !
    ! Notes:
    !  1. the result is in Cartesian coordinates;
    !  2. the beam is assumed to start in a vacuum, so N = 1.

    use const_and_precisions, only : degree
    use gray_params,          only : antenna_parameters

    ! subroutine arguments
    type(antenna_parameters), intent(in) :: params
    real(wp_), intent(out) :: N(3)

    ! local variables
    real(wp_) :: r, Nr, Nphi, a, b

    R = hypot(params%pos(1), params%pos(2))
    a = params%alpha * degree
    b = params%beta  * degree

    ! N̅ in cylindrical coordinates using the poloidal
    ! and toroidal launch angles.
    Nr   = -cos(a)*cos(b)
    Nphi = sin(b)
 
    ! Convert to Cartesian coordinates
    ! Notes:
    !  1. The unit vectors in Cartesian coordinates are:
    !       R^ = x^ cosφ + y^ sinφ
    !       φ^ = x^ sinφ - y^ sinφ
    !  2. tanφ = y/x ⇒ { cosφ = x/√(x²+y²) = x/R
    !                  { sinφ = y/√(x²+y²) = x/R
    N(1) = (Nr*params%pos(1) - Nphi*params%pos(2)) / R
    N(2) = (Nr*params%pos(2) + Nphi*params%pos(1)) / R
    N(3) = -cos(b)*sin(a)
  end subroutine launchangles2n


  subroutine xgygcoeff(fghz, k0, Bres, xgcn)
    ! Given the EC wave frequency computes:
    !
    !   1. vacuum wavevector `k0` (k₀ = ω/c),
    !   2. resonant magnetic field `Bres` (qB/m = ω),
    !   3. adimensional `xgcn` parameter (X = ω_p²/ω² = nq²/ε₀mω²).
    use const_and_precisions, only : qe=>ecgs_, me=>mecgs_, &
                                     vc=>ccgs_, pi, wce1_

    ! subroutine arguments
    real(wp_), intent(in)  :: fghz
    real(wp_), intent(out) :: k0, Bres, xgcn

    ! local variables
    real(wp_) :: omega

    omega = 2.0e9_wp_*pi*fghz ! [rad/s]
    k0 = omega/vc             ! [rad/cm]

    ! yg = Btot/Bres
    Bres = omega/wce1_ ! [T]

    ! xg = xgcn*dens19
    xgcn = 4.0e13_wp_ * pi * qe**2/(me * omega**2) ! [10^-19 m^3]
  end subroutine xgygcoeff


  pure subroutine gaussian_eikonal(Q, r, z, S, grad_S, hess_S)
    ! Computes the complex eikonal, gradient and Hessian of a Gaussian beam

    ! subroutine arguments
    complex(wp_), intent(in)  :: Q(2,2)      ! complex curvature tensor Q(z)
    real(wp_),    intent(in)  :: r(2), z     ! coordinates
    complex(wp_), intent(out) :: S           ! complex eikonal
    complex(wp_), intent(out) :: grad_S(3)   ! gradient of imaginary eikonal
    complex(wp_), intent(out) :: hess_S(3,3) ! Hessian of imaginary eikonal

    ! local variables
    complex(wp_) :: Q2(2,2), Q3(2,2)
    Q2 = matmul(Q, Q)
    Q3 = matmul(Q, Q2)

    ! The most general Gaussian beam is given by
    !
    !   E̅(r̅,z) = e̅(z) exp[-ik₀z - ik₀r̅⋅Q(z)⋅r̅/2 + iη(z)]
    !
    ! where:
    !   - e̅(z) is the amplitude
    !   - Q(z) = Q₀ (I + zQ₀)⁻¹ is the complex curvature tensor
    !   - η(z) is the Gouy phase
    !   - k₀ = ω/c
    !   - r̅ = [x, y]
    !
    ! See https://doi.org/10.1364/AO.52.006030 for more details.
    !
    ! Since the eikonal ansatz in GRAY is E̅(x̅,t) = e̅(r̅) exp[-ik₀S(x̅) + iωt],
    ! the eikonal for a Gaussian beam is:
    !
    !   S(x̅) = z + ½ r̅⋅Q(z)⋅r̅ - η(z)/k₀
    !
    S = z + dot_product(r, matmul(Q, r))/2

    ! Splitting up x̅ into (r̅, z), the gradient of S can be written as
    !
    !   ∇S(r̅,z) = [Q(z)⋅r̅, 1 + ½ r̅⋅dQ/dz⋅r̅]
    !
    ! For dQ/dz, we use the matrix identity dA/dz = - A⁻¹ dA/dz A⁻¹:
    !
    !   dQ/dz = Q₀ d/dz (I + zQ₀)⁻¹
    !         = - Q₀ (I + zQ₀)⁻¹ Q₀ (I + zQ₀)⁻¹
    !         = - Q(z)²
    !
    ! So ∇S(r̅,z) = [Q(z)⋅r̅, 1 - ½ r̅⋅Q(z)²⋅r̅]
    !
    grad_S = [matmul(Q, r), 1 - dot_product(r, matmul(Q2, r))/2]

    ! Carrying on, the Hessian splits into 4 blocks (2x2, 2x1, 1x2, 1x1):
    !
    !  HS(r̅,z) = k₀ [ Q(z),          dQ/dz⋅r̅] = k₀ [ Q(z),    -Q(z)²⋅r̅]
    !               [-Q(z)²⋅r̅, -½ r̅⋅dQ²/dz⋅r̅]      [-Q(z)²⋅r̅,   r̅⋅Q³⋅r̅]
    !
    ! having used dQ²/dz = 2Q dQ/dz = -2Q³.
    !
    hess_S(1:2,:) = reshape([Q, -matmul(Q2, r)], [2,3])
    hess_S(3,  :) = [-matmul(Q2, r), dot_product(r, matmul(Q3, r))]

    ! Note: this ignores the Gouy phase, but for the computation of the
    ! initial conditions (x̅, N̅_R=∇S_R, N̅_I=∇S_I for each ray) it is
    ! inconsequential because:
    !
    !  1. the initial conditions are given at z=const, so η=const as well
    !  2. N̅_I = ∇S_I does not depend on η(z)
    !  3. N̅_R does depend on η(z), but is obtained indirectly as the
    !     solution of the vacuum quasioptical dispersion relation:
    !
    !      { N̅_R ⋅ N̅_I   = 0
    !      { N_R² - N_I² = N₀²
    !
    ! The latter implies N̅_I and N̅_R are not exactly consistent with
    ! a Gaussian beam, but this is a necessary trade-off.
  end subroutine gaussian_eikonal

end module beams