build_periodic_zero_mode_height_plan Subroutine

  subroutine build_periodic_zero_mode_height_plan(source_heights, area_xy, plan, status, message)
    real(dp), intent(in) :: source_heights(:, :)
    real(dp), intent(in) :: area_xy
    type(periodic_zero_mode_plan_type), intent(out) :: plan
    integer(i32), intent(out) :: status
    character(len=*), intent(out) :: message
    real(dp), allocatable :: heights(:), unique_heights(:)
    real(dp) :: scale, tolerance
    integer(i32) :: elem, vertex, count, nelem

    status = periodic_zero_mode_invalid
    message = ''
    if (.not. ieee_is_finite(area_xy) .or. area_xy <= 0.0_dp) then
      message = 'periodic zero-mode area_xy must be finite and positive'
      return
    end if
    if (size(source_heights, 1) /= 3 .or. size(source_heights, 2) <= 0) then
      message = 'periodic zero-mode heights must have shape (3, nsource)'
      return
    end if
    if (any(.not. ieee_is_finite(source_heights))) then
      message = 'periodic zero-mode panel heights must be finite'
      return
    end if
    nelem = int(size(source_heights, 2), i32)
    allocate (heights(3*nelem))
    do elem = 1_i32, nelem
      heights(3*elem - 2:3*elem) = source_heights(:, elem)
    end do
    call sort_real(heights)
    scale = max(1.0_dp, maxval(abs(heights)))
    tolerance = 128.0_dp*epsilon(1.0_dp)*scale
    allocate (unique_heights(size(heights)))
    count = 1_i32
    unique_heights(1) = heights(1)
    do vertex = 2, size(heights)
      if (abs(heights(vertex) - unique_heights(count)) > tolerance) then
        count = count + 1_i32
        unique_heights(count) = heights(vertex)
      end if
    end do

    plan%nelem = nelem
    plan%nbreak = count
    plan%area_xy = area_xy
    allocate (plan%break_z(count), plan%panel(nelem))
    plan%break_z = unique_heights(1:count)
    do elem = 1, nelem
      plan%panel(elem)%z = source_heights(:, elem)
      call sort_three(plan%panel(elem)%z)
      plan%panel(elem)%horizontal = &
        abs(plan%panel(elem)%z(3) - plan%panel(elem)%z(1)) <= tolerance
      do vertex = 1, 3
        plan%panel(elem)%break_index(vertex) = nearest_break(plan%break_z, plan%panel(elem)%z(vertex))
      end do
    end do
    status = periodic_zero_mode_ok
  end subroutine build_periodic_zero_mode_height_plan