panel_potential_sum_translated_sp Subroutine

  subroutine panel_potential_sum_translated_sp( &
    count, vx, vy, vz, nx, ny, nz, edge, ox, oy, oz, area, index, charge, geometry, axes, shift1, shift2, &
    target_local, target, phi &
    )
    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) :: phi
    real(sp) :: out(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) :: px, py, pz, l1, l2, l3, d1, d2, d3, integral, 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, px, py, pz, l1, l2, l3, d1, d2, d3, integral, 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)
        px = tx - height*nx(m)
        py = ty - height*ny(m)
        pz = tz - height*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)
        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))
        d1 = (vx(1, m) - px)*ox(1, m) + (vy(1, m) - py)*oy(1, m) + (vz(1, m) - pz)*oz(1, m)
        d2 = (vx(2, m) - px)*ox(2, m) + (vy(2, m) - py)*oy(2, m) + (vz(2, m) - pz)*oz(2, m)
        d3 = (vx(3, m) - px)*ox(3, m) + (vy(3, m) - py)*oy(3, m) + (vz(3, m) - pz)*oz(3, m)
        integral = -height*omega + d1*l1 + d2*l2 + d3*l3
        out(k) = scale*integral
      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)
          phi = phi + potential/k_coulomb
          cycle
        end if
        phi = phi + real(out(k), dp)/k_coulomb
      end do
    end do
  end subroutine panel_potential_sum_translated_sp