编辑代码

program B_Spline_Basis_Functions
  implicit none

  integer, parameter :: N = 10  ! 控制点数
  integer, parameter :: K = 3   ! B样条阶数
  real(8), parameter :: U(N+K+1) = (/0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0/)  ! 节点矢量
  real(8), parameter :: t = 4.5  ! 要计算基函数的位置

  real(8) :: N_k(N, K)  ! 存储所有K次基函数的值

  ! 计算K次B样条基函数
  call Compute_B_Spline_Basis_Functions(N, K, U, t, N_k)

  ! 打印结果
  write(*,*) "B样条基函数："
  do i = 1, N
    write(*,*) "N_", i, "(", K, ")(", t, ") = ", N_k(i, K)
  end do

contains

  subroutine Compute_B_Spline_Basis_Functions(N, K, U, t, N_k)
    integer, intent(in) :: N, K
    real(8), intent(in) :: U(N+K+1)
    real(8), intent(in) :: t
    real(8), intent(out) :: N_k(N, K)

    integer :: i, j
    real(8) :: alpha, beta

    ! 计算0次B样条基函数
    do i = 1, N
      if (U(i) <= t .and. t < U(i + 1)) then
        N_k(i, 0) = 1.0
      else
        N_k(i, 0) = 0.0
      end if
    end do

    ! 递推计算K次B样条基函数
    do j = 1, K
      do i = 1, N
        alpha = 0.0
        beta = 0.0
        if (U(i+j) - U(i) > 0.0) then
          alpha = (t - U(i)) / (U(i+j) - U(i))
        end if
        if (U(i+j+1) - U(i+1) > 0.0) then
          beta = (U(i+j+1) - t) / (U(i+j+1) - U(i+1))
        end if
        N_k(i, j) = alpha * N_k(i, j-1) + beta * N_k(i+1, j-1)
      end do
    end do

  end subroutine Compute_B_Spline_Basis_Functions

end program B_Spline_Basis_Functions