src/beams.f90: document read_beam{0,1} formats

This commit is contained in:
Michele Guerini Rocco 2023-03-26 16:03:57 +02:00
parent 038864a84f
commit 69308901ee
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
2 changed files with 86 additions and 36 deletions

View File

@ -7,6 +7,24 @@ contains
subroutine read_beam0(params, 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 z z φ
!
! where:
! - f is the frequency (GHz)
! - x, y, z are the launcher position (mm)
! - w,w are the beam waists in the two principal directions (mm)
! - z,z are the positions of the foci in the two principal
! directions (mm)
! - φ is the of the ellipse (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
@ -22,7 +40,7 @@ contains
! local variables
integer :: u
integer :: err
real(wp_) :: ak0,zrcsi,zreta
real(wp_) :: k0, w0(2), z0(2), z_R(2), phi
u = get_free_unit(unit)
@ -33,26 +51,50 @@ contains
call exit(1)
end if
read(u, *) params%fghz
read(u, *) params%pos
read(u, *) params%w, params%ri, params%phi(1)
read(u, *) params%fghz ! Wave frequency (GHz)
read(u, *) params%pos ! Launcher position (Cartesian coordinates)
read(u, *) w0, & ! Beam waists (in the ξ,η axes)
z0, & ! Foci positions (in the ξ,η axes)
phi ! Ellipse axes angle
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
! Wavevector k
k0 = 2.0e9_wp_* pi * params%fghz / vc
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)
! Rayleigh ranges in the ξ, η directions:
! z_R = ½kw²
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 + (z0/z_R)**2)
! Curvature in the ξ, η directions at z=0:
! K = 1/R = (z-z)/[(z-z)² + z_R²]
params%ri = -z0 / (z0**2 + z_R**2)
! Ellipse axes angle
params%phi = phi
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 α.
!
! 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 (useless)
! - α, β are the poloidal and toroidal launch angles (deg)
! - x, y, z are the launcher position (mm)
! - w,w are the beam waists in the two principal directions (mm)
! - k,k are the wavefront curvatures in the two principal
! directions (mm¹)
! - φ_w, φ_R are the angles of the amplitude and phase ellipses (deg)
use gray_params, only : antenna_parameters
use simplespline, only : spli, difcs
@ -742,9 +784,13 @@ contains
end subroutine read_beam2
subroutine launchangles2n(params, anv)
subroutine launchangles2n(params, N)
! Given the wave launcher `params` computes the initial
! wavevector `anv`, defined by = ck̅/ω, in cartesian coordinates.
! wavevector = 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
@ -753,31 +799,38 @@ contains
! subroutine arguments
type(antenna_parameters), intent(in) :: params
real(wp_), intent(out) :: anv(3)
real(wp_), intent(out) :: N(3)
! local variables
real(wp_) :: r, anr, anphi, a, b
real(wp_) :: r, Nr, Nphi, a, b
r = sqrt(params%pos(1)**2 + params%pos(2)**2)
a = degree*params%alpha
b = degree*params%beta
R = hypot(params%pos(1), params%pos(2))
a = params%alpha * degree
b = params%beta * degree
! Angles α, β in a local reference system
! as proposed by Gribov et al.
anr = -cos(b)*cos(a)
anphi = sin(b)
! in cylindrical coordinates using the poloidal
! and toroidal launch angles.
Nr = -cos(a)*cos(b)
Nphi = 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
! 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, ak0, bres, xgcn)
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 = ω),
! 2. resonant magnetic field `Bres` (qB/m = ω),
! 3. adimensional `xgcn` parameter (X = ω_p²/ω² = nq²/ε²).
use const_and_precisions, only : qe=>ecgs_, me=>mecgs_, &
vc=>ccgs_, pi, wce1_
@ -785,16 +838,16 @@ contains
! subroutine arguments
real(wp_), intent(in) :: fghz
real(wp_), intent(out) :: ak0, bres, xgcn
real(wp_), intent(out) :: k0, Bres, xgcn
! local variables
real(wp_) :: omega
omega = 2.0e9_wp_*pi*fghz ! [rad/s]
ak0 = omega/vc ! [rad/cm]
k0 = omega/vc ! [rad/cm]
! yg = btot/bres
bres = omega/wce1_ ! [T]
! yg = Btot/Bres
Bres = omega/wce1_ ! [T]
! xg = xgcn*dens19
xgcn = 4.0e13_wp_ * pi * qe**2/(me * omega**2) ! [10^-19 m^3]

View File

@ -92,10 +92,7 @@ contains
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
! Compute the initial cartesian wavevector (anv0)
! from launch angles α,β and the position x:
! NR(α, β, x)
! (α, β, x)
! Nz(α, β, x)
! from launch angles α,β and the position
call launchangles2n(params%antenna, anv0)
! Initialise the ray variables (beamtracing)