From 69308901eef52ac753d38566c95d850c2f89fb85 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sun, 26 Mar 2023 16:03:57 +0200 Subject: [PATCH] src/beams.f90: document read_beam{0,1} formats --- src/beams.f90 | 117 +++++++++++++++++++++++++++++++++------------- src/gray_core.f90 | 5 +- 2 files changed, 86 insertions(+), 36 deletions(-) diff --git a/src/beams.f90 b/src/beams.f90 index f3b6d42..b6c8759 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -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 = ½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 + (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 n̅ = ck̅/ω, in cartesian coordinates. + ! 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 @@ -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) + ! N̅ 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²/ε₀mω²). 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] diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 4990456..df6bb4e 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -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₀) - ! Nφ(α, β, x₀) - ! Nz(α, β, x₀) + ! from launch angles α,β and the position call launchangles2n(params%antenna, anv0) ! Initialise the ray variables (beamtracing)