904 lines
36 KiB
Fortran
904 lines
36 KiB
Fortran
module beams
|
||
use const_and_precisions, only : wp_, one
|
||
implicit none
|
||
|
||
contains
|
||
|
||
subroutine read_beam0(params, err, unit)
|
||
! 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 utils, only : get_free_unit
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(antenna_parameters), intent(inout) :: params
|
||
integer, intent(out) :: err
|
||
integer, intent(in), optional :: unit
|
||
|
||
! local variables
|
||
integer :: u
|
||
real(wp_) :: k0, w0(2), d0(2), z_R(2), phi
|
||
|
||
u = get_free_unit(unit)
|
||
|
||
open(unit=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, unit)
|
||
! 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 : get_free_unit,locate
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(antenna_parameters), intent(inout) :: params
|
||
integer, intent(out) :: err
|
||
integer, optional, intent(in) :: unit
|
||
|
||
! 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
|
||
|
||
u = get_free_unit(unit)
|
||
|
||
open(unit=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, nisteer, 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, unit)
|
||
! 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 contains the number N of beams defined in the file
|
||
! 2. The following N lines contain the following data for each of the
|
||
! N beams:
|
||
! ID,iox,f,na,nb
|
||
! 3. The meaning of the data is
|
||
! - ID is a label to identify the beam
|
||
! - iox=1,2 is a flag to select O-mode (1) or X-mode (2) polarization
|
||
! - f is the wave frequency (GHz)
|
||
! - n*m is the number of rows of the corresponding table below with
|
||
! the beam parameters.
|
||
! * If the beam steering is fixed: na=nb=1
|
||
! * For a steering around a single axis:
|
||
! na>1 and nb=1, or na=1 and nb>1.
|
||
! * For a steering with two independent axes:
|
||
! n>1 and m>1.
|
||
! The row number l of the following table is mapped to two indexes
|
||
! i,j (1≤i≤na, 1≤j≤nb) via l=i+na*(j-1), i.e., the index i "runs
|
||
! faster".
|
||
! Index i is assumed to be associated to a steering mainly in the
|
||
! poloidal direction and the poloidal launch angle α(i,j) must be
|
||
! monotonous along its first dimension.
|
||
! Index j is assumed to be associated to a steering mainly in the
|
||
! toroidal direction and the toroidal launch angle β(i,j) must be
|
||
! monotonous along its second dimension.
|
||
! 4. The rest of file is a sequence of N tables with the following
|
||
! columns:
|
||
! α, β, x₀, y₀, z₀, w₁, w₂, k₁, k₂, φ_w, φ_R
|
||
! 5. The meaning of the columns is
|
||
! - α, β 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 utils, only : get_free_unit, intlin, locate
|
||
use reflections, only : inside
|
||
use dierckx, only : curfit, splev, surfit, bispev
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(antenna_parameters), intent(inout) :: params
|
||
integer, intent(in) :: beamid
|
||
integer, intent(out), optional :: err
|
||
integer, intent(in), optional :: unit
|
||
|
||
! local variables
|
||
character(len=20) :: beamname
|
||
integer :: u
|
||
integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta
|
||
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, &
|
||
xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, &
|
||
ypolygC, xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC
|
||
real(wp_), DIMENSION(4) :: xvert, yvert
|
||
real(wp_), dimension(1) :: fi
|
||
integer, parameter :: kspl=1
|
||
real(wp_), parameter :: sspl=0.01_wp_
|
||
|
||
u = get_free_unit(unit)
|
||
|
||
open(unit=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, params%iox, params%fghz, nalpha, nbeta
|
||
jumprow = jumprow+nalpha*nbeta
|
||
end do
|
||
! c====================================================================================
|
||
!
|
||
! beam of interest
|
||
read(u,*) beamname, params%iox, params%fghz, nalpha, nbeta
|
||
!
|
||
! 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), &
|
||
xpolyg(npolyg), ypolyg(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
|
||
xpolyg(1:nxcoord) = xcoord(1:nxcoord)
|
||
ypolyg(1:nxcoord) = ycoord(1:nxcoord)
|
||
!
|
||
! c====================================================================================
|
||
do i=1,nycoord-2
|
||
xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord)
|
||
xpolyg(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1)
|
||
ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord)
|
||
ypolyg(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1)
|
||
end do
|
||
! c====================================================================================
|
||
do i=1,nxcoord
|
||
xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1)
|
||
ypolyg(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(inside(xpolyg,ypolyg,npolyg,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(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), &
|
||
ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), &
|
||
xpolygD(nycoord), ypolygD(nycoord), &
|
||
xoutA(nxcoord+3), youtA(nxcoord+3), xoutB(nycoord+3), &
|
||
youtB(nycoord+3), xoutC(nxcoord+3), youtC(nxcoord+3))
|
||
! coordinates of vertices v1,v2,v3,v4
|
||
xvert(1) = xpolyg(1)
|
||
xvert(2) = xpolyg(nxcoord)
|
||
xvert(3) = xpolyg(nxcoord+nycoord-1)
|
||
xvert(4) = xpolyg(2*nxcoord+nycoord-2)
|
||
yvert(1) = ypolyg(1)
|
||
yvert(2) = ypolyg(nxcoord)
|
||
yvert(3) = ypolyg(nxcoord+nycoord-1)
|
||
yvert(4) = ypolyg(2*nxcoord+nycoord-2)
|
||
! coordinates of side A,B,C,D
|
||
xpolygA = xpolyg(1:nxcoord)
|
||
ypolygA = ypolyg(1:nxcoord)
|
||
xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1)
|
||
ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1)
|
||
xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
|
||
ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
|
||
xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg)
|
||
xpolygD(nycoord) = xpolyg(1)
|
||
ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg)
|
||
ypolygD(nycoord) = ypolyg(1)
|
||
! contour of outside regions A (1,2), B(3,4), C(5,6)
|
||
xoutA = huge(one)
|
||
xoutA(1:nxcoord) = xpolygA
|
||
xoutA(nxcoord+3) = xvert(1)
|
||
youtA = -huge(one)
|
||
youtA(1:nxcoord) = ypolygA
|
||
youtA(nxcoord+1) = yvert(2)
|
||
xoutB = huge(one)
|
||
xoutB(1:nycoord) = xpolygB
|
||
xoutB(nycoord+1) = xvert(3)
|
||
youtB = huge(one)
|
||
youtB(1:nycoord) = ypolygB
|
||
youtB(nycoord+3) = yvert(2)
|
||
xoutC = -huge(one)
|
||
xoutC(1:nxcoord) = xpolygC
|
||
xoutC(nxcoord+3) = xvert(3)
|
||
youtC = huge(one)
|
||
youtC(1:nxcoord) = ypolygC
|
||
youtC(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(inside(xoutA,youtA,nxcoord+3,xcoord0,ycoord0)) then
|
||
in = 1
|
||
if(xcoord0.GT.xvert(2)) then
|
||
in = 2
|
||
end if
|
||
else if(inside(xoutB,youtB,nycoord+3,xcoord0,ycoord0)) then
|
||
in = 3
|
||
if(ycoord0.GT.yvert(3)) then
|
||
in = 4
|
||
end if
|
||
else if(inside(xoutC,youtC,nxcoord+3,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(xpolygA,nxcoord,xcoord0,ii)
|
||
! find corresponding y value on side A for xcoord position
|
||
call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0)
|
||
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(ypolygB,nycoord,ycoord0,ii)
|
||
call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0)
|
||
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(xpolygC,nxcoord,xcoord0,ii)
|
||
call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0)
|
||
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(ypolygD,nycoord,ycoord0,ii)
|
||
call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0)
|
||
incheck = 1
|
||
CASE (8)
|
||
! params%alpha and params%beta outside table range
|
||
xcoord0 = xvert(1)
|
||
ycoord0 = yvert(1)
|
||
ii = 1
|
||
END SELECT
|
||
! c----------------------------------------------------------------------------------
|
||
!
|
||
deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, &
|
||
xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC)
|
||
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, xpolyg, ypolyg, 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
|
||
|
||
implicit none
|
||
|
||
! 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_
|
||
implicit none
|
||
|
||
! 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
|
||
|
||
end module beams
|