804 lines
33 KiB
Fortran
804 lines
33 KiB
Fortran
module beams
|
||
use const_and_precisions, only : wp_, one
|
||
implicit none
|
||
|
||
contains
|
||
|
||
subroutine read_beam0(params, unit)
|
||
! Reads the wave launcher parameters for the simple case
|
||
! where w(z) and 1/R(z) are fixed.
|
||
|
||
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(in), optional :: unit
|
||
|
||
! local variables
|
||
integer :: u
|
||
integer :: err
|
||
real(wp_) :: ak0,zrcsi,zreta
|
||
|
||
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")
|
||
call exit(1)
|
||
end if
|
||
|
||
read(u, *) params%fghz
|
||
read(u, *) params%pos
|
||
read(u, *) params%w, params%ri, params%phi(1)
|
||
close(u)
|
||
|
||
ak0 = 2.0e9_wp_* pi * params%fghz / vc
|
||
zrcsi = 0.5_wp_ * ak0 * params%w(1)**2
|
||
zreta = 0.5_wp_ * ak0 * params%w(2)**2
|
||
|
||
params%w(1) = params%w(1) * sqrt(1.0_wp_ + (params%ri(1)/zrcsi)**2)
|
||
params%w(2) = params%w(2) * sqrt(1.0_wp_ + (params%ri(2)/zreta)**2)
|
||
params%ri(1) = -params%ri(1) / (params%ri(1)**2 + zrcsi**2)
|
||
params%ri(2) = -params%ri(2) / (params%ri(2)**2 + zreta**2)
|
||
params%phi(2) = params%phi(1)
|
||
end subroutine read_beam0
|
||
|
||
|
||
subroutine read_beam1(params, unit)
|
||
! Reads the wave launcher parameters for the case
|
||
! where w(z, α) and 1/R(z, α) depend on the launcher angle α.
|
||
|
||
use gray_params, only : antenna_parameters
|
||
use simplespline, only : spli, difcs
|
||
use utils, only : get_free_unit,locate
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(antenna_parameters), intent(inout) :: params
|
||
integer, intent(in), optional :: unit
|
||
|
||
! local variables
|
||
integer :: u, iopt, ier, nisteer, i, k, ii
|
||
real(wp_) :: steer, dal
|
||
real(wp_), dimension(:), allocatable :: &
|
||
alphastv, betastv, x00v, y00v, &
|
||
z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, &
|
||
cbeta, cx0, cy0, cz0, cwaist1, cwaist2, &
|
||
crci1, crci2, cphi1, cphi2
|
||
integer :: err
|
||
|
||
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")
|
||
call exit(1)
|
||
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), cbeta(4*nisteer), &
|
||
cx0(4*nisteer), cy0(4*nisteer), cz0(4*nisteer), &
|
||
cwaist1(4*nisteer), cwaist2(4*nisteer), crci1(4*nisteer), &
|
||
crci2(4*nisteer), cphi1(4*nisteer), cphi2(4*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._wp_ * rci1v
|
||
rci2v = 10._wp_ * rci2v
|
||
|
||
iopt = 0
|
||
call difcs(alphastv, betastv, nisteer, iopt, cbeta, ier)
|
||
call difcs(alphastv, waist1v, nisteer, iopt, cwaist1, ier)
|
||
call difcs(alphastv, rci1v, nisteer, iopt, crci1, ier)
|
||
call difcs(alphastv, waist2v, nisteer, iopt, cwaist2, ier)
|
||
call difcs(alphastv, rci2v, nisteer, iopt, crci2, ier)
|
||
call difcs(alphastv, phi1v, nisteer, iopt, cphi1, ier)
|
||
call difcs(alphastv, phi2v, nisteer, iopt, cphi2, ier)
|
||
call difcs(alphastv, x00v, nisteer, iopt, cx0, ier)
|
||
call difcs(alphastv, y00v, nisteer, iopt, cy0, ier)
|
||
call difcs(alphastv, z00v, nisteer, iopt, cz0, ier)
|
||
|
||
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 = spli(cbeta, nisteer, k, dal)
|
||
params%pos(1) = spli(cx0, nisteer, k, dal)
|
||
params%pos(2) = spli(cy0, nisteer, k, dal)
|
||
params%pos(3) = spli(cz0, nisteer, k, dal)
|
||
params%w(1) = spli(cwaist1, nisteer, k, dal)
|
||
params%w(2) = spli(cwaist2, nisteer, k, dal)
|
||
params%ri(1) = spli(crci1, nisteer, k, dal)
|
||
params%ri(2) = spli(crci2, nisteer, k, dal)
|
||
params%phi(1) = spli(cphi1, nisteer, k, dal)
|
||
params%phi(2) = spli(cphi2, nisteer, 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, cbeta, &
|
||
cx0, cy0, cz0, cwaist1, cwaist2, &
|
||
crci1, crci2, cphi1, cphi2)
|
||
|
||
end subroutine read_beam1
|
||
|
||
|
||
subroutine read_beam2(params, beamid, unit)
|
||
! Reads the wave launcher parameters for the general case
|
||
! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β.
|
||
|
||
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(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_
|
||
integer :: err
|
||
|
||
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")
|
||
call exit(1)
|
||
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, anv)
|
||
! Given the wave launcher `params` computes the initial
|
||
! wavevector `anv`, defined by n̅ = ck̅/ω, in cartesian coordinates.
|
||
|
||
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) :: anv(3)
|
||
|
||
! local variables
|
||
real(wp_) :: r, anr, anphi, a, b
|
||
|
||
r = sqrt(params%pos(1)**2 + params%pos(2)**2)
|
||
a = degree*params%alpha
|
||
b = degree*params%beta
|
||
|
||
! Angles α, β in a local reference system
|
||
! as proposed by Gribov et al.
|
||
anr = -cos(b)*cos(a)
|
||
anphi = sin(b)
|
||
|
||
anv(1) = (anr*params%pos(1) - anphi*params%pos(2))/r ! = anx
|
||
anv(2) = (anr*params%pos(2) + anphi*params%pos(1))/r ! = any
|
||
anv(3) = -cos(b)*sin(a) ! = anz
|
||
end subroutine launchangles2n
|
||
|
||
|
||
subroutine xgygcoeff(fghz, ak0, 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) :: ak0, bres, xgcn
|
||
|
||
! local variables
|
||
real(wp_) :: omega
|
||
|
||
omega = 2.0e9_wp_*pi*fghz ! [rad/s]
|
||
ak0 = 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
|