编辑代码

program query_lines_within_sphere
    implicit none

    integer, parameter :: dp = selected_real_kind(15) ! 双精度实数
    real(dp) :: x0 = 1.0_dp, y0 = 1.0_dp, z0 = 1.0_dp ! 球体中心点坐标
    real(dp) :: r = 10.0_dp ! 球体半径

    ! 线段数组，每行存储一条线段的两个端点坐标
    real(dp), dimension(3, 2, 3) :: lines = reshape([ &
        0.0_dp, 0.0_dp, 0.0_dp, &
        1.0_dp, 1.0_dp, 1.0_dp, &
        1.0_dp, 1.0_dp, 1.0_dp, &
        2.0_dp, 2.0_dp, 2.0_dp, &
        0.0_dp, 1.0_dp, 0.0_dp, &
        1.0_dp, 0.0_dp, 1.0_dp &
    ], [3, 2, 3])

    integer :: indices(3) ! 存储在球体内的线段的下标
    integer :: i, j, k
    real(dp) :: d1, d2

    j = 0
    do i = 1, size(lines, 1)
        k = 0
        d1 = sqrt((lines(i, 1, 1)-x0)**2 + (lines(i, 1, 2)-y0)**2 + (lines(i, 1, 3)-z0)**2)
        d2 = sqrt((lines(i, 2, 1)-x0)**2 + (lines(i, 2, 2)-y0)**2 + (lines(i, 2, 3)-z0)**2)
        if (d1 <= r .and. d2 <= r) then
            j = j + 1
            indices(j) = i
        end if
    end do

    write(*, *) "Indices of lines within sphere:", indices(1:j)

end program query_lines_within_sphere