panel_field_sum_translated Subroutine

  subroutine panel_field_sum_translated( &
    count, vx, vy, vz, nx, ny, nz, edge, ox, oy, oz, area, index, charge, geometry, target, ex, ey, ez &
    )
    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) :: ex, ey, ez
    real(dp) :: outx(panel_chunk), outy(panel_chunk), outz(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) :: l1, l2, l3, fx, fy, fz, 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, l1, l2, l3, fx, fy, fz, 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)
        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
        ! Keep unused lanes non-singular without branching around the vector math calls.
        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))
        fx = nx(m)*omega
        fy = ny(m)*omega
        fz = nz(m)*omega
        fx = fx + ox(1, m)*l1
        fy = fy + oy(1, m)*l1
        fz = fz + oz(1, m)*l1
        fx = fx + ox(2, m)*l2
        fy = fy + oy(2, m)*l2
        fz = fz + oz(2, m)*l2
        fx = fx + ox(3, m)*l3
        fy = fy + oy(3, m)*l3
        fz = fz + oz(3, m)*l3
        outx(k) = scale*fx
        outy(k) = scale*fy
        outz(k) = scale*fz
      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 &
            )
          outx(k) = field(1)
          outy(k) = field(2)
          outz(k) = field(3)
        end if
        ex = ex + outx(k)/k_coulomb
        ey = ey + outy(k)/k_coulomb
        ez = ez + outz(k)/k_coulomb
      end do
    end do
  end subroutine panel_field_sum_translated