用户登陆
站点日历
站点统计
最新评论
空间三点画园 [ 日期：2007-04-06 ] [ 来自：空间三点画园
算法：
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
'空间三点绘制圆***********************************************************空间三点绘制圆**********
'
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
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
[本日志由 田草 于 2007-04-13 09:19 PM 编辑]
|
|
暂时没有评论
| 发表评论 - 不要忘了输入验证码哦！ |