panel_potential_sum_translated Subroutine

  subroutine panel_potential_sum_translated( &
    count, vx, vy, vz, nx, ny, nz, edge, ox, oy, oz, area, index, charge, geometry, target, phi &
    )
    integer(i32), intent(in) :: count
    real(dp), intent(in) :: vx(3, *), vy(3, *), vz(3, *), nx(*), ny(*), nz(*), edge(3, *)
    real(dp), intent(in) :: ox(3, *), oy(3, *), oz(3, *), area(*)
    integer(i32), intent(in) :: index(*)
    real(dp), intent(in) :: charge(:)
    type(panel_geometry_type), intent(in) :: geometry(:)
    real(dp), intent(in) :: target(3)
    real(dp), intent(inout) :: phi
    real(dp) :: out(panel_chunk)
    logical :: active(panel_chunk), on_plane(panel_chunk)
    real(dp) :: ax, ay, az, bx, by, bz, cx, cy, cz, ra, rb, rc, height, numerator, denominator, omega
    real(dp) :: px, py, pz, l1, l2, l3, d1, d2, d3, integral, scale, potential, field(3)
    logical :: edge_hit, evaluate_lane
    integer(i32) :: chunk_start, chunk_size, k, m, idx

    do chunk_start = 1_i32, count, panel_chunk
      chunk_size = min(panel_chunk, count - chunk_start + 1_i32)
      edge_hit = .false.
      !$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) reduction(.or.:edge_hit)
      do k = 1_i32, chunk_size
        m = chunk_start + k - 1_i32
        active(k) = abs(charge(index(m))) > tiny(1.0_dp)
        scale = k_coulomb*charge(index(m))/area(m)
        height = (target(1) - vx(1, m))*nx(m) + (target(2) - vy(1, m))*ny(m) + (target(3) - vz(1, m))*nz(m)
        px = target(1) - height*nx(m)
        py = target(2) - height*ny(m)
        pz = target(3) - height*nz(m)
        on_plane(k) = abs(height) <= 64.0_dp*epsilon(1.0_dp)*max(edge(1, m), edge(2, m), edge(3, m))
        evaluate_lane = active(k) .and. .not. on_plane(k)
        ax = vx(1, m) - target(1)
        ay = vy(1, m) - target(2)
        az = vz(1, m) - target(3)
        bx = vx(2, m) - target(1)
        by = vy(2, m) - target(2)
        bz = vz(2, m) - target(3)
        cx = vx(3, m) - target(1)
        cy = vy(3, m) - target(2)
        cz = vz(3, m) - target(3)
        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
        omega = -2.0_dp*atan2(merge(numerator, 0.0_dp, evaluate_lane), merge(denominator, 1.0_dp, evaluate_lane))
        l1 = edge(1, m)/(ra + rb)
        l2 = edge(2, m)/(rb + rc)
        l3 = edge(3, m)/(rc + ra)
        edge_hit = edge_hit .or. (evaluate_lane .and. max(l1, l2, l3) >= 1.0_dp)
        l1 = merge(l1, 0.0_dp, evaluate_lane .and. l1 < 1.0_dp)
        l2 = merge(l2, 0.0_dp, evaluate_lane .and. l2 < 1.0_dp)
        l3 = merge(l3, 0.0_dp, evaluate_lane .and. l3 < 1.0_dp)
        l1 = log((1.0_dp + l1)/(1.0_dp - l1))
        l2 = log((1.0_dp + l2)/(1.0_dp - l2))
        l3 = log((1.0_dp + l3)/(1.0_dp - 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
        integral = integral + d1*l1
        integral = integral + d2*l2
        integral = integral + d3*l3
        out(k) = scale*integral
      end do
      if (edge_hit) error stop 'panel target lies on an edge or vertex.'
      do k = 1_i32, chunk_size
        if (.not. active(k)) cycle
        if (on_plane(k)) then
          m = chunk_start + k - 1_i32
          idx = index(m)
          call panel_potential_field_on_vertices( &
            geometry(idx), vx(:, m), vy(:, m), vz(:, m), charge(idx), target, potential, field &
            )
          out(k) = potential
        end if
        phi = phi + out(k)/k_coulomb
      end do
    end do
  end subroutine panel_potential_sum_translated