fix simple limiter in analytic equilibrium

In the case of analytic equilibrium without a limiter contour, the
simple limiter was built incorrectly due to an unnecessary conversion
from cm (the equilibrium data are already in metre).

src/const_and_precisions.f90

@@ -57,5 +57,6 @@ module const_and_precisions
   real(wp_), parameter :: mc2_   = me_*c_**2/kev_ ! Electron rest energy, keV
   real(wp_), parameter :: mu0inv = 1._wp_/mu0_    ! Inverse magnetic permeability of vacuum, m/H
   real(wp_), parameter :: wce1_  = e_/me_         ! ECR (angular) frequency / magnetic field, rad/s/T
+  real(wp_), parameter :: cm = 0.01_wp_           ! Converts cm to m
 
 end module const_and_precisions

src/coreprofiles.f90

@@ -451,9 +451,9 @@ contains
     ! invalid as they are given assuming a vacuum (N=1)
     if (present(launch_pos)) then
       block
-        use equilibrium, only : equinum_psi
+        use equilibrium,          only : equinum_psi
+        use const_and_precisions, only : cm
         real(wp_) :: R, Z, psi
-        real(wp_), parameter :: cm = 1.0e-2_wp_
 
         ! Convert from cartesian to cylindrical coordinates
         R = hypot(launch_pos(1), launch_pos(2)) * cm  ! R = √(x²+y²)

src/gray_jetto1beam.f90

@@ -60,15 +60,21 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
       call set_globals(params)
     end block init_params
 
-    ! Set a simple limiter
+    ! Set a simple limiter following the boundary of the data grid
     simple_limiter: block
-      use reflections, only : range2rect
+      use reflections,          only : range2rect
-      use limiter,     only : limiter_set_globals=>set_globals
+      use limiter,              only : limiter_set_globals=>set_globals
+      use const_and_precisions, only : cm
+
+      ! Avoid clipping out the launcher
+      real(wp_) :: R_launch, R_max
+      R_launch = norm2(params%antenna%pos(1:2)) * cm
+      R_max = max(R_launch, R(mr))
+
+      ! Convert to a list of R,z
+      call range2rect(xmin=params%misc%rwall, xmax=R_max, ymin=z(1), ymax=z(mz), &
+                      xv=data%equilibrium%rlim, yv=data%equilibrium%zlim)
 
-      real(wp_) :: r0m
-      r0m = norm2(params%antenna%pos(1:2)) * 0.01_wp_
-      call range2rect(params%misc%rwall, max(r0m, r(mr)), z(1), z(mz), &
-                      data%equilibrium%rlim, data%equilibrium%zlim)
       ! Set the global variables of `limiter`
       call limiter_set_globals(data%equilibrium)
     end block simple_limiter

src/main.f90

@@ -384,8 +384,9 @@ contains
 
   subroutine init_misc(params, data)
     ! Performs miscellanous initial tasks, before the gray_main subroutine.
-    use reflections, only : range2rect
+    use reflections,          only : range2rect
-    use limiter,     only : limiter_set_globals=>set_globals
+    use limiter,              only : limiter_set_globals=>set_globals
+    use const_and_precisions, only : cm
 
     implicit none
 
@@ -393,27 +394,49 @@ contains
     type(gray_parameters), intent(inout) :: params
     type(gray_data),       intent(inout) :: data
 
-    ! Build a basic limiter when one is not provided by the EQDSK
+    ! Build a basic limiter if one is not provided by equilibrium.txt
     if (.not. allocated(data%equilibrium%rlim) &
         .or. params%raytracing%ipass < 0) then
-      block
-        real(wp_) :: rmxm, r0m, z0m, dzmx
-        r0m  = norm2(params%antenna%pos(1:2)) * 0.01_wp_
-        dzmx = abs(params%raytracing%ipass) * &
-               params%raytracing%dst * params%raytracing%nstep * 0.01_wp_
-        z0m  = params%antenna%pos(3) * 0.01_wp_
 
-        allocate(data%equilibrium%rlim(5))
+      ! Restore sign of ipass
-        allocate(data%equilibrium%zlim(5))
+      params%raytracing%ipass = abs(params%raytracing%ipass)
-        params%raytracing%ipass = abs(params%raytracing%ipass)
+
+      allocate(data%equilibrium%rlim(5))
+      allocate(data%equilibrium%zlim(5))
+      block
+        real(wp_) :: R_launch, z_launch, R_max, delta_z
+
+        ! Launcher coordinates
+        R_launch = norm2(params%antenna%pos(1:2)) * cm
+        z_launch = params%antenna%pos(3) * cm
+
+        ! Max vertical distance a ray can travel
+        delta_z = abs(params%raytracing%ipass) * &
+                  params%raytracing%dst * params%raytracing%nstep * cm
+
+        ! Max radius, either due to the plasma extent or equilibrium grid
         if (params%equilibrium%iequil < 2) then
-          rmxm = (data%equilibrium%rv(1) + data%equilibrium%rv(2)) * 0.01_wp_
+          ! analytic equilibrium, use R₀+a
+          R_max = data%equilibrium%rv(1) + data%equilibrium%rv(2)
         else
-          rmxm = data%equilibrium%rv(size(data%equilibrium%rv))
+          ! numeric equilibrium, use max R of the grid
+          R_max = data%equilibrium%rv(size(data%equilibrium%rv))
         end if
-        call range2rect(params%misc%rwall, max(r0m, rmxm), &
+
-                        z0m - dzmx, z0m + dzmx,            &
+        ! Avoid clipping out the launcher
-                        data%equilibrium%rlim, data%equilibrium%zlim)
+        R_max = max(R_launch, R_max)
+
+        ! Convert to a list of R,z:
+        !
+        !  (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz)
+        !                       ║4   3║
+        !                       ║     ║
+        !                       ║1   2║
+        !  (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz)
+        !
+        call range2rect(xmin=params%misc%rwall, xmax=R_max,               &
+                        ymin=z_launch - delta_z, ymax=z_launch + delta_z, &
+                        xv=data%equilibrium%rlim, yv=data%equilibrium%zlim)
       end block
     end if