Merge pull request #2 from jonathanschilling/no_abscab_in_belicu

No abscab in belicu

src/NESTOR/belicu.f90

@@ -28,6 +28,10 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
 
   real(rprec), dimension(3, nuv2) :: eval_pos, magnetic_field
 
+  ! If .true., use ABSCAB for computing the magnetic field contribution
+  ! due to the net toridal current modeled as a filament along the magnetic axis.
+  logical, parameter :: use_abscab_for_axis_current = .false.
+
   ! B_External due to LIne CUrrent
 
   ! net toroidal plasma current in A
@@ -36,37 +40,79 @@ SUBROUTINE belicu(torcur, bx, by, bz, cos1, sin1, rp, zp)
   ! loops over source geometry
   i = 0
   DO kper = 1, nvper
-     DO kv = 1, nv
-        i = i + 1
+      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
+          ! 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)
-
-  DO j = 1, nuv2
-    ! evaluation positions
-    eval_pos(1, j) = rp(j) * cos1(j)
-    eval_pos(2, j) = rp(j) * sin1(j)
-    eval_pos(3, j) = zp(j)
-  end do
-
-  ! initialize target array
-  magnetic_field = 0.0_dp
-
-  ! use ABSCAB to compute the line-current-along-axis magnetic field contribution
-  call magneticFieldPolygonFilament(nvper * nv + 1, xpts, current, &
-                                    nuv2, eval_pos, magnetic_field)
-
-  bx(:) = magnetic_field(1,:)
-  by(:) = magnetic_field(2,:)
-  bz(:) = magnetic_field(3,:)
+  xpts(1, nvp + 1) = xpts(1, 1)
+  xpts(2, nvp + 1) = xpts(2, 1)
+  xpts(3, nvp + 1) = xpts(3, 1)
+
+  if (use_abscab_for_axis_current) then
+
+    DO j = 1, nuv2
+        ! evaluation positions
+        eval_pos(1, j) = rp(j) * cos1(j)
+        eval_pos(2, j) = rp(j) * sin1(j)
+        eval_pos(3, j) = zp(j)
+    end do
+
+    ! initialize target array
+    magnetic_field = 0.0_dp
+
+    ! use ABSCAB to compute the line-current-along-axis magnetic field contribution
+    call magneticFieldPolygonFilament(nvper * nv + 1, xpts, current,  &
+                                      nuv2, eval_pos, magnetic_field)
+
+    bx(:) = magnetic_field(1,:)
+    by(:) = magnetic_field(2,:)
+    bz(:) = magnetic_field(3,:)
+
+  else ! use_abscab_for_axis_current
+
+    ! initialize target array
+    bx = 0
+    by = 0
+    bz = 0
+
+    ! 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+1)-xpts(:,i)
+        L = norm2(dvec)
+
+        ! loop over evaluation points
+        DO j = 1,nuv2
+            ! xpt == evaluation position
+            xpt(1) = rp(j) * cos1(j)
+            xpt(2) = rp(j) * sin1(j)
+            xpt(3) = zp(j)
+
+            Ri_vec = xpt - xpts(:,i)
+            Ri = norm2(Ri_vec)
+            Rf = norm2(xpt - xpts(:,i + 1))
+            Ri_p_Rf = Ri + Rf
+
+            ! 1.0e-7 == mu0/4 pi
+            Bmag = 1.0E-7_dp * current * 2.0_dp * Ri_p_Rf / ( Ri * Rf * (Ri_p_Rf*Ri_p_Rf - L*L) )
+
+            ! cross product of L*hat(eps)==dvec with Ri_vec, scaled by Bmag
+            bx(j) = bx(j) + Bmag * (dvec(2)*Ri_vec(3) - dvec(3)*Ri_vec(2))
+            by(j) = by(j) + Bmag * (dvec(3)*Ri_vec(1) - dvec(1)*Ri_vec(3))
+            bz(j) = bz(j) + Bmag * (dvec(1)*Ri_vec(2) - dvec(2)*Ri_vec(1))
+        end do
+
+    END DO
+
+  end if ! use_abscab_for_axis_current
 
 END SUBROUTINE belicu