编辑代码

program B_Spline_Curve
  implicit none

  integer, parameter :: N = 6  ! 控制点数
  integer, parameter :: K = 3  ! B样条阶数
  real(8), parameter :: U(N+K+1) = (/0.0, 0.0, 0.0, 0.1, 0.3, 0.6, 0.8, 1.0, 1.0, 1.0/)  ! 节点矢量

  real(8), dimension(N) :: ControlPoints  ! 控制点坐标
  real(8), dimension(N) :: Weight        ! 控制点权重
  real(8) :: t  ! 参数值
  real(8) :: X, Y, Z  ! 物理坐标

  ! 初始化控制点坐标和权重
  ControlPoints = (/1.0, 2.0, 3.0, 4.0, 5.0, 6.0/)
  Weight = (/1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)

  ! 要计算B样条曲线上的参数点
  t = 0.5  ! 你可以更改t的值

  ! 计算B样条曲线上的物理坐标
  call Compute_B_Spline_Curve(N, K, U, ControlPoints, Weight, t, X, Y, Z)

  ! 打印结果
  write(*,*) "B样条曲线上的物理坐标 (X, Y, Z) 在 t = ", t, " 处："
  write(*,*) "X = ", X
  write(*,*) "Y = ", Y
  write(*,*) "Z = ", Z

contains

  subroutine Compute_B_Spline_Curve(N, K, U, ControlPoints, Weight, t, X, Y, Z)
    integer, intent(in) :: N, K
    real(8), intent(in) :: U(N+K+1)
    real(8), intent(in) :: ControlPoints(N)
    real(8), intent(in) :: Weight(N)
    real(8), intent(in) :: t
    real(8), intent(out) :: X, Y, Z

    integer :: i, j
    real(8) :: N_k(K)
    real(8) :: Denominator

    X = 0.0
    Y = 0.0
    Z = 0.0

    do i = 1, N
      ! 计算K次B样条基函数
      call Compute_B_Spline_Basis(N, K, U, t, N_k)

      ! 计算权重和基函数的乘积
      Denominator = 0.0
      do j = 1, K
        Denominator = Denominator + N_k(j) * Weight(i)
      end do

      ! 更新物理坐标
      X = X + N_k(K) * ControlPoints(i) * Weight(i) / Denominator
      Y = Y + N_k(K) * ControlPoints(i) * Weight(i) / Denominator
      Z = Z + N_k(K) * ControlPoints(i) * Weight(i) / Denominator
    end do

  end subroutine Compute_B_Spline_Curve

  subroutine Compute_B_Spline_Basis(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(K)

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

    ! 递推计算B样条基函数
    ! 计算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

end program B_Spline_Curve