空间三点画园


空间三点画园


算法: 

      1、在空间取得三点  P1、P2、P3
      2、判断三点是否共线,是继续,否则直接提示三点共线,退出。
      3、空间三点求平面方程,(Ax+By+Cz+D=0)
      4、计算线段P1P2和线段P2P3的中点坐标,P12 和P23.
      4、求出直线P1P2 的方程{x=(x1-x2)t+x1;y=(y1-y2)t+y1;z=(z1-z2)t+z1)} 
      5、 线段P1P2的中垂线( 在P1P2P3平面内)向量和 直线P1P2 的向量乘积为零。
           设中垂线上的一点为P4
           (x1-x2)*(x12-x4)+(y1-y2)*(y12-y4)+(z1-z2)*(z12-z4)=0
           这是P1P2的中垂面方程。
           他和平面P1P2P3的交线就是中垂线方程。
       6、同样的方面计算出P2P3的中垂线方程。
       7、两个中垂线的交点就是三点圆的圆心。两个中垂线方程组成三元一次方程组。及得到圆的圆心坐标。
       8、圆心到三点中的任意一点的距离及是圆的半径。
       9、在当前视图平面内绘制等半径的一个圆,圆心就是计算出的圆心。
       10、将圆旋转到P1P2P3的平面内。
              a 圆的旋转轴确定,两个平面的交线既是其旋转轴。
              c 计算旋转角度,通过两个平面的法向量的成绩计算两个平面的夹角。
                
       理论上不错 但是 三元一次方程 问题 还是没有解决 ,也不知道 哪儿弄错了。就是没有实现。现在想通过 圆心到三点的距离相等求圆心的坐标。
   如果谁会三元一次方程解决,请帮帮我。

  圆心到三点的距离相等。
  |P0P1= |P0P2|=|P0P3|

还是通过一开始的方法 我还是实现了,但是在解方程得时候 还有点问题。方程的一般试的某个系数为零的情况没有解决我想我会慢慢完善它的。

'********************************************************************************************

'空间三点绘制圆***********************************************************空间三点绘制圆**********
'
Sub Circle3D3P()
    Dim P1 As Variant
    Dim P2 As Variant
    Dim P3 As Variant
    ThisDrawing.Utility.Prompt "空间三点绘制圆:"
    P1 = ThisDrawing.Utility.GetPoint(, "第一点:")
    ThisDrawing.Utility.InitializeUserInput 1, ""
    P2 = ThisDrawing.Utility.GetPoint(, "第二点:")
    ThisDrawing.Utility.InitializeUserInput 1, ""
    P3 = ThisDrawing.Utility.GetPoint(, "第三点:")
    
    ThisDrawing.ModelSpace.AddLine P1, P2
    ThisDrawing.ModelSpace.AddLine P2, P3
    ThisDrawing.ModelSpace.AddLine P3, P1
    ThisDrawing.ModelSpace.AddPoint P1
    ThisDrawing.ModelSpace.AddPoint P2
    ThisDrawing.ModelSpace.AddPoint P3
    '判断三点不共线
    If ThreeP_IsOnline(P1, P2, P3) = True Then
       ThisDrawing.Utility.Prompt "所选三点共线" & vbCrLf
       Exit Sub
    End If
    '计算三点空间平面方程
    Dim A1  As Double, B1 As Double, C1 As Double, D1 As Double
    KJPMFC P1, P2, P3, A1, B1, C1, D1
    'P1P2的中点
    Dim P_12 As Variant
    P_12 = CenterPoint(P1, P2)
    ThisDrawing.ModelSpace.AddPoint P_12
    'P1P2 的向量
    Dim P12(0 To 2) As Double
    P12(0) = P1(0) - P2(0):      P12(1) = P1(1) - P2(1):            P12(2) = P1(2) - P2(2)
    '中垂面的法向量 {P12(0),P12(1),P12(2)}
    '中垂线向量和P1P2向量积为零(P1(0)-P2(0))*(x-P_12(0))+(P1(1)-P2(1))*(y-P_12(1))+(P1(2)-P2(2))*(z-P_12(2))=0
    '中垂面的系数
    Dim A2  As Double, B2 As Double, C2 As Double, D2 As Double
    A2 = P12(0)
    B2 = P12(1)
    C2 = P12(2)
    D2 = -P12(0) * P_12(0) - P12(1) * P_12(1) - P12(2) * P_12(2)
    'P2P3的中点
    Dim P_23 As Variant
    P_23 = CenterPoint(P2, P3)
    ThisDrawing.ModelSpace.AddPoint P_23
    'P2P3 的向量
    Dim P23(0 To 2) As Double
    P23(0) = P2(0) - P3(0):      P23(1) = P2(1) - P3(1):            P23(2) = P2(2) - P3(2)
    Dim A3  As Double, B3 As Double, C3 As Double, d3 As Double
    A3 = P23(0)
    B3 = P23(1)
    C3 = P23(2)
    d3 = -P23(0) * P_23(0) - P23(1) * P_23(1) - P23(2) * P_23(2)
    '两个中垂面和三点所在平面的分别的交线正好组成一个三元一次方程组。
    '    高斯消元法解方程,这个不能出现方程的系数为零的情况。(以后再完善)。
    Dim Xishu(1 To 3, 1 To 3) As Double, Changshu(1 To 3) As Double
    Xishu(1, 1) = A1: Xishu(1, 2) = B1: Xishu(1, 3) = C1:  Changshu(1) = -D1
    Xishu(2, 1) = A2: Xishu(2, 2) = B2: Xishu(2, 3) = C2:  Changshu(2) = -D2
    Xishu(3, 1) = A3: Xishu(3, 2) = B3: Xishu(3, 3) = C3:  Changshu(3) = -d3
    Dim x0 As Double, y0 As Double, z0 As Double
    XXFCZ Xishu, Changshu
    x0 = Changshu(1)
    y0 = Changshu(2)
    z0 = Changshu(3)
    'Point3d(x0,y0,z0) 即为圆心坐标
    ThisDrawing.ModelSpace.AddPoint Point3D(x0, y0, z0)
    '验证圆心坐标是否正确(到三点的距离相等)。
    Dim R1 As Double, R2 As Double, R3 As Double
    R1 = P2PDistance(P1, Point3D(x0, y0, z0))
    R2 = P2PDistance(P2, Point3D(x0, y0, z0))
    R3 = P2PDistance(P3, Point3D(x0, y0, z0))
    MsgBox "R1=" & R1
    MsgBox "R2=" & R2
    MsgBox "R3=" & R3
    '圆心是否满足平面方程?
    MsgBox A1 * x0 + B1 * y0 + C1 * z0 + D1
    Dim C3P As AcadCircle
    '过圆心在当前视图平面(x,y平面)内画圆。
    Set C3P = ThisDrawing.ModelSpace.AddCircle(Point3D(x0, y0, z0), R1)
    '计算圆的旋转轴
    '   计算圆和平面的交线方程 上的另一点(已经存在一点就是圆心)
    '      首先计算圆的方程
    '      圆的方程应该是x=R*cos(angle)+x0,y=R*sin(angle)+y0,z=z0
     Dim M1(0 To 2) As Double
     '该圆垂直于z轴的一个平面,他和三点所在平面的交线当然也垂直于z轴,z坐标不变,x和y坐标可以连续变化。
     '   令这个点的x坐标为x0+100,z坐标不变即满足了圆所在的平面方程,通过P1P2P3平面方程计算y坐标
     M1(0) = x0 + 100
     M1(1) = -(A1 * (x0 + 100) + C1 * z0 + D1) / B1
     M1(2) = z0
     ThisDrawing.ModelSpace.AddPoint M1
     '画出旋转轴
     ThisDrawing.ModelSpace.AddLine M1, Point3D(x0, y0, z0)
     
    '计算圆和平面的夹角
    '    还是要圆计算平面的方程的一般试
    '        在圆上找两个点
    Dim M2(0 To 2) As Double, M3(0 To 2) As Double
    M2(0) = R1 * Cos(Atn(1)) + x0
    M2(1) = R1 * Sin(Atn(1)) + y0
    M2(2) = z0
    M3(0) = R1 * Cos(Atn(1) * 2) + x0
    M3(1) = R1 * Sin(Atn(1) * 2) + y0
    M3(2) = z0
    Dim A4 As Double, B4 As Double, C4 As Double, D4 As Double
    KJPMFC Point3D(x0, y0, z0), M2, M3, A4, B4, C4, D4
    '平面的夹角
    Dim CosA As Double
    CosA = (A1 * A4 + B1 * B4 + C1 * C4) / (Sqr(A1 ^ 2 + B1 ^ 2 + C1 ^ 2) * Sqr(A4 ^ 2 + B4 ^ 2 + C4 ^ 2))
    Dim J As Double
    J = Arccos(CosA)
    
    C3P.Rotate3d Point3D(x0, y0, z0), M1, -J
    
End Sub
'**********************************************************************************************

'线性方程组的解法***********************************************************线性方程组得解法********
'
Sub XXFCZ(ByRef A() As Double, ByRef B() As Double)
    '高斯消元法
    'A(i,j)是系数
    'B(i)是右端项
    Dim n As Long
    Dim k As Long
    Dim l As Double
    Dim J As Long
    Dim sum As Double
    n = UBound(B)
    For k = 1 To n - 1
        For i = k + 1 To n
            l = A(i, k) / A(k, k)
            For J = k + 1 To n
                A(i, J) = A(i, J) - l * A(k, J)
            Next J
            B(i) = B(i) - l * B(k)
        Next i
    Next k '以上是消元过程
    B(n) = B(n) / A(n, n)
    For i = n - 1 To 1 Step -1
        sum = 0
        For J = i + 1 To n
            sum = sum + A(i, J) * B(J)
        Next J
        B(i) = (B(i) - sum) / A(i, i)
    Next i '以上是回代过程
End Sub


下面的 高斯消元法就可以解决 方程系数为零的情况。它可以求解任意 线性方程。A() 为方程系数,N为方程的未知数,B为方程右侧的常熟。
Function GaoSi(A(), N, B())
    Dim Ipiv(50), INdxr(50), Indxc(50)
    Dim J
    Dim K
    Dim L
    Dim LL
    Dim Dum
    Dim Pivinv
    For J = 1 To N
        Ipiv(J) = 0
    Next J
    Dim i
    Dim Big
    Dim Irow
    Dim Icol

    For i = 1 To N
        Big = 0#
        For J = 1 To N
            If Ipiv(J) <> 1 Then
                For K = 1 To N
                    If Ipiv(K) = 0 Then
                        If Abs(A(J, K)) >= Big Then
                            Big = Abs(A(J, K))
                            Irow = J
                            Icol = K
                        End If
                    ElseIf Ipiv(K) > 1 Then
                        MsgBox "异常矩阵"
                    End If
                Next K
            End If
        Next J
        Ipiv(Icol) = Ipiv(Icol) + 1
        If Irow <> Icol Then
            For L = 1 To N
                Dum = A(Irow, L)
                A(Irow, L) = A(Icol, L)
                A(Icol, L) = Dum
            Next L
            Dum = B(Irow)
            B(Irow) = B(Icol)
            B(Icol) = Dum
        End If
        INdxr(i) = Irow
        Indxc(i) = Icol
        If A(Icol, Icol) = 0# Then MsgBox "异常矩阵"
        Pivinv = 1# / A(Icol, Icol)
        A(Icol, Icol) = 1#
        For L = 1 To N
            A(Icol, L) = A(Icol, L) * Pivinv
        Next L
        B(Icol) = B(Icol) * Pivinv
        For LL = 1 To N
            If LL <> Icol Then
                Dum = A(LL, Icol)
                A(LL, Icol) = 0#
                For L = 1 To N
                    A(LL, L) = A(LL, L) - A(Icol, L) * Dum
                Next L
                B(LL) = B(LL) - B(Icol) * Dum
            End If
        Next LL
    Next i
    For L = N To 1 Step -1
        If INdxr(L) <> Indxc(L) Then
            For K = 1 To N
                Dum = A(K, INdxr(L))
                A(K, INdxr(L)) = A(K, Indxc(L))
                A(K, Indxc(L)) = Dum
            Next K
        End If
    Next L
End Function




vb 表达试计算

系统优化批处理

欢迎关注微信公众账号ByCAD