panel_field_sum_translated_sp Subroutine

  subroutine panel_field_sum_translated_sp( &
    count, vx, vy, vz, nx, ny, nz, edge, ox, oy, oz, area, index, charge, geometry, axes, shift1, shift2, &
    target_local, target, ex, ey, ez &
    )
    integer(i32), intent(in) :: count
    real(sp), intent(in) :: vx(3, *), vy(3, *), vz(3, *), nx(*), ny(*), nz(*), edge(3, *)
    real(sp), intent(in) :: ox(3, *), oy(3, *), oz(3, *), area(*)
    integer(i32), intent(in) :: index(*), axes(2)
    real(dp), intent(in) :: charge(:), shift1(*), shift2(*)
    type(panel_geometry_type), intent(in) :: geometry(:)
    real(dp), intent(in) :: target_local(3), target(3)
    real(dp), intent(inout) :: ex, ey, ez
    real(sp) :: outx(panel_chunk), outy(panel_chunk), outz(panel_chunk), tx, ty, tz
    logical :: active(panel_chunk), fallback(panel_chunk)
    real(sp) :: ax, ay, az, bx, by, bz, cx, cy, cz, ra, rb, rc, height, numerator, denominator, omega
    real(sp) :: l1, l2, l3, fx, fy, fz, scale
    real(dp) :: potential, field(3), shift(3)
    logical :: evaluate_lane
    integer(i32) :: chunk_start, chunk_size, k, m, idx

    tx = real(target_local(1), sp)
    ty = real(target_local(2), sp)
    tz = real(target_local(3), sp)
    do chunk_start = 1_i32, count, panel_chunk
      chunk_size = min(panel_chunk, count - chunk_start + 1_i32)
      !$omp simd private(m, ax, ay, az, bx, by, bz, cx, cy, cz, ra, rb, rc, height, numerator, denominator) &
      !$omp   private(omega, l1, l2, l3, fx, fy, fz, scale, evaluate_lane)
      do k = 1_i32, chunk_size
        m = chunk_start + k - 1_i32
        active(k) = abs(charge(index(m))) > tiny(1.0_dp)
        scale = real(k_coulomb*charge(index(m)), sp)/area(m)
        height = (tx - vx(1, m))*nx(m) + (ty - vy(1, m))*ny(m) + (tz - vz(1, m))*nz(m)
        ax = vx(1, m) - tx
        ay = vy(1, m) - ty
        az = vz(1, m) - tz
        bx = vx(2, m) - tx
        by = vy(2, m) - ty
        bz = vz(2, m) - tz
        cx = vx(3, m) - tx
        cy = vy(3, m) - ty
        cz = vz(3, m) - tz
        ra = sqrt(ax*ax + ay*ay + az*az)
        rb = sqrt(bx*bx + by*by + bz*bz)
        rc = sqrt(cx*cx + cy*cy + cz*cz)
        numerator = ax*(by*cz - bz*cy) + ay*(bz*cx - bx*cz) + az*(bx*cy - by*cx)
        denominator = ra*rb*rc + (ax*bx + ay*by + az*bz)*rc + (bx*cx + by*cy + bz*cz)*ra + &
                      (cx*ax + cy*ay + cz*az)*rb
        l1 = edge(1, m)/(ra + rb)
        l2 = edge(2, m)/(rb + rc)
        l3 = edge(3, m)/(rc + ra)
        ! Points on the panel plane or within about one percent of an edge lose single-precision
        ! digits in log((1+r)/(1-r)); those lanes are redone with the double kernel below.
        fallback(k) = active(k) .and. (abs(height) <= 64.0_sp*epsilon(1.0_sp)*max(edge(1, m), edge(2, m), edge(3, m)) &
                                       .or. max(l1, l2, l3) >= near_edge_single_limit)
        evaluate_lane = active(k) .and. .not. fallback(k)
        omega = -2.0_sp*atan2(merge(numerator, 0.0_sp, evaluate_lane), merge(denominator, 1.0_sp, evaluate_lane))
        l1 = merge(l1, 0.0_sp, evaluate_lane)
        l2 = merge(l2, 0.0_sp, evaluate_lane)
        l3 = merge(l3, 0.0_sp, evaluate_lane)
        l1 = log((1.0_sp + l1)/(1.0_sp - l1))
        l2 = log((1.0_sp + l2)/(1.0_sp - l2))
        l3 = log((1.0_sp + l3)/(1.0_sp - l3))
        fx = nx(m)*omega + ox(1, m)*l1 + ox(2, m)*l2 + ox(3, m)*l3
        fy = ny(m)*omega + oy(1, m)*l1 + oy(2, m)*l2 + oy(3, m)*l3
        fz = nz(m)*omega + oz(1, m)*l1 + oz(2, m)*l2 + oz(3, m)*l3
        outx(k) = scale*fx
        outy(k) = scale*fy
        outz(k) = scale*fz
      end do
      do k = 1_i32, chunk_size
        if (.not. active(k)) cycle
        if (fallback(k)) then
          m = chunk_start + k - 1_i32
          idx = index(m)
          shift = 0.0_dp
          if (axes(1) > 0_i32) shift(axes(1)) = shift(axes(1)) + shift1(m)
          if (axes(2) > 0_i32) shift(axes(2)) = shift(axes(2)) + shift2(m)
          call panel_potential_field_translated(geometry(idx), shift, charge(idx), target, potential, field)
          ex = ex + field(1)/k_coulomb
          ey = ey + field(2)/k_coulomb
          ez = ez + field(3)/k_coulomb
          cycle
        end if
        ex = ex + real(outx(k), dp)/k_coulomb
        ey = ey + real(outy(k), dp)/k_coulomb
        ez = ez + real(outz(k), dp)/k_coulomb
      end do
    end do
  end subroutine panel_field_sum_translated_sp