revert to simple belicu implementation

src/NESTOR/belicu.f90

@@ -32,29 +32,36 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
 
   ! net toroidal plasma current in A
   current = torcur/mu0
-  
+
   ! initialize target array
   bx = 0
   by = 0
   bz = 0
 
-  ! first point (at index 0) is equal to last point --> closed curve
-  xpts(1, 0) = raxis_nestor(nv)*(cosper(nvper)*cosuv(nv) - sinper(nvper)*sinuv(nv))
-  xpts(2, 0) = raxis_nestor(nv)*(sinper(nvper)*cosuv(nv) + cosper(nvper)*sinuv(nv))
-  xpts(3, 0) = zaxis_nestor(nv)
-
   ! loops over source geometry
-  i = 1
+  i = 0
   DO kper = 1, nvper
      DO kv = 1, nv
+        i = i + 1
 
         ! xpts == xpt of _s_ource (current filament)
         xpts(1, i) = raxis_nestor(kv)*(cosper(kper)*cosuv(kv) - sinper(kper)*sinuv(kv))
         xpts(2, i) = raxis_nestor(kv)*(sinper(kper)*cosuv(kv) + cosper(kper)*sinuv(kv))
         xpts(3, i) = zaxis_nestor(kv)
+      END DO
+  END DO
+
+  ! last point is equal to first point --> closed curve
+  xpts(1, nvp + 1) = xpts(1, 1)
+  xpts(2, nvp + 1) = xpts(2, 1)
+  xpts(3, nvp + 1) = xpts(3, 1)
+
+  ! iterate over all wire segments that make up the axis;
+  ! the number of wire segments is one less than number of points of the closed loop
+  DO i = 1, nvper * nv
 
         ! filament geometry: from current point (R_i == xpts(:,i)) to previous point (R_f == xpts(:,i-1))
-        dvec = xpts(:,i)-xpts(:,i-1)
+        dvec = xpts(:,i+1)-xpts(:,i)
         L = norm2(dvec)
 
         ! loop over evaluation points
@@ -64,9 +71,9 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
            xpt(2) = rp(j) * sin1(j)
            xpt(3) = zp(j)
 
-           Ri_vec = xpt - xpts(:,i-1)
+           Ri_vec = xpt - xpts(:,i)
            Ri = norm2(Ri_vec)
-           Rf = norm2(xpt - xpts(:,i))
+           Rf = norm2(xpt - xpts(:,i + 1))
            Ri_p_Rf = Ri + Rf
 
            ! 1.0e-7 == mu0/4 pi
@@ -78,8 +85,6 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
            bz(j) = bz(j) + Bmag * (dvec(1)*Ri_vec(2) - dvec(2)*Ri_vec(1))
         end do
 
-        i = i + 1
-     END DO
   END DO
 
 END SUBROUTINE belicu

src/NESTOR/bextern.f90

@@ -18,6 +18,7 @@ SUBROUTINE bextern(plascur, wint)
   REAL(rprec), DIMENSION(nuv2), INTENT(in) :: wint
 
   INTEGER :: i
+  logical :: dbgout_active
 
   ! exterior Neumann problem
 
@@ -30,6 +31,19 @@ SUBROUTINE bextern(plascur, wint)
   ! This sets brad, bphi and bz to the interpolated field from the mgrid.
   CALL becoil(r1b, z1b, bvac(1,1), bvac(1,2), bvac(1,3))
 
+  dbgout_active = open_dbg_context("vac1n_bextern", num_eqsolve_retries)
+  if (dbgout_active) then
+
+    ! these are only from the mgrid at this point
+    call add_real_2d("mgrid_brad", nv, nu3, brad)
+    call add_real_2d("mgrid_bphi", nv, nu3, bphi)
+    call add_real_2d("mgrid_bz",   nv, nu3, bz)
+
+    ! This should be in Amperes.
+    call add_real("axis_current", plascur/mu0)
+
+  end if ! dbgout_active
+
   ! COMPUTE B (ON PLASMA BOUNDARY) FROM NET TOROIDAL PLASMA CURRENT
   ! THE NET CURRENT IS MODELLED AS A WIRE AT THE MAGNETIC AXIS, AND THE
   ! BIOT-SAVART LAW IS USED TO COMPUTE THE FIELD AT THE PLASMA SURFACE
@@ -58,8 +72,12 @@ SUBROUTINE bextern(plascur, wint)
   ! NOTE: BEXN == NP*F = -B0 dot [Xu cross Xv] NP        (see PKM, Eq. 2.13)
   bexni(:nuv2) = wint(:nuv2)*bexn(:nuv2)*pi2*pi2
 
-  if (open_dbg_context("vac1n_bextern", num_eqsolve_retries)) then
+  if (dbgout_active) then
+
+    ! axis geometry used in belicu
+    call add_real_2d("xpts_axis", 3, nvper * nv + 1, xpts)
 
+    ! these are now mgrid + axis-current
     call add_real_2d("brad", nv, nu3, brad)
     call add_real_2d("bphi", nv, nu3, bphi)
     call add_real_2d("bz",   nv, nu3, bz)

src/NESTOR/data/vacmod.f90

@@ -229,8 +229,9 @@ subroutine allocate_nestor
   IF (i .ne. 0) STOP 'allocation error in bextern'
 
   ! from tolicu
-  ! need 0:nvp for "virtual" point at index 0 which is equal to last point for closed curve
-  ALLOCATE (xpts(3,0:nvp), stat=i)
+  ! need nvp+1 for "virtual" point at last index,
+  ! which is equal to first point for a closed curve
+  ALLOCATE (xpts(3,nvp+1), stat=i)
   IF (i .ne. 0) STOP ' allocation error in tolicu'
 
   ! from scalpot