triangular pyramid, trigonal pyramid offset-3d-Plane-Intersection 通过 vbasolver
triangular pyramid, trigonal pyramid offset-3d-Plane-Intersection By vbasolver
triangular pyramid, trigonal pyramid offset-3d-Plane-Intersection By vba求解器
x=0.250 y=0.250 z=0.250‖内切球‖中心
每个平面偏移0.25 我挑战vba求解器
为什么全行 0.250 0.250 0.250
我要全行 0.250 0.250 0.250
为什么6行8行错误
结果很奇怪
请告诉我如何解决它。
Const x0 = 0
Const y0 = 0
Const z0 = 2
Const x1 = 0
Const y1 = 0
Const z1 = 0
Const x2 = 1
Const y2 = 0
Const z2 = 0
Const x3 = 0
Const y3 = 1
Const z3 = 0
Const r0 = 0.25
Const r1 = 0.25
Const r2 = 0.25
Const r3 = 0.25
Function myR1C1toA1(i, j)
myR1C1toA1 = Application.ConvertFormula("R" & i & "C" & j, xlR1C1, xlA1)
End Function
Function myPlane(Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz)
a = (By - Ay) * (Cz - Az) - (Cy - Ay) * (Bz - Az)
b = (Bz - Az) * (Cx - Ax) - (Cz - Az) * (Bx - Ax)
c = (Bx - Ax) * (Cy - Ay) - (Cx - Ax) * (By - Ay)
d = -(a * Ax + b * Ay + c * Az)
myPlane = Array(a, b, c, d)
End Function
Function myOffset3PlaneIntersection(irow, x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3, r1, r2, r3)
Dim v1 As Variant
Dim v2 As Variant
Dim v3 As Variant
'
v1 = myPlane(x0, y0, z0, x2, y2, z2, x3, y3, z3)
v2 = myPlane(x0, y0, z0, x3, y3, z3, x1, y1, z1)
v3 = myPlane(x0, y0, z0, x1, y1, z1, x2, y2, z2)
'
my1 = "(Abs(" & v1(0) & " * " & myR1C1toA1(irow, 1) & "+ " & v1(1) & " * " & myR1C1toA1(irow, 2) & "+ " & v1(2) & " * " & myR1C1toA1(irow, 3) & "+ " & v1(3) & ") / Sqrt(" & v1(0) & " ^ 2 + " & v1(1) & " ^ 2 +" & v1(2) & " ^ 2) - " & myR1C1toA1(irow, 1) & ")"
my2 = "(Abs(" & v2(0) & " * " & myR1C1toA1(irow, 1) & "+ " & v2(1) & " * " & myR1C1toA1(irow, 2) & "+ " & v2(2) & " * " & myR1C1toA1(irow, 3) & "+ " & v2(3) & ") / Sqrt(" & v2(0) & " ^ 2 + " & v2(1) & " ^ 2 +" & v2(2) & " ^ 2) - " & myR1C1toA1(irow, 2) & ")"
my3 = "(Abs(" & v3(0) & " * " & myR1C1toA1(irow, 1) & "+ " & v3(1) & " * " & myR1C1toA1(irow, 2) & "+ " & v3(2) & " * " & myR1C1toA1(irow, 3) & "+ " & v3(3) & ") / Sqrt(" & v3(0) & " ^ 2 + " & v3(1) & " ^ 2 +" & v3(2) & " ^ 2) - " & myR1C1toA1(irow, 3) & ")"
Range(myR1C1toA1(irow, 5)).Formula = "=" & my1 & "^ 2 +" & my2 & "^ 2 +" & my3 & "^ 2"
'
Dim ws As Worksheet: Set ws = ActiveSheet
SolverReset
SolverOk setCell:=ws.Range(myR1C1toA1(irow, 5)), _
MaxMinVal:=3, _
ByChange:=ws.Range(myR1C1toA1(irow, 1) & ":" & myR1C1toA1(irow, 3)), _
EngineDesc:="GRG Nonlinear"
SolverSolve UserFinish:=True
End Function
Sub myFormat()
Columns("A:E").Select
Selection.NumberFormatLocal = "0.000_ "
Range("A5").Select
End Sub
Sub aaa_Main()
Dim myXYZ As Variant
ActiveSheet.Cells.Clear
MsgBox "Solver Start"
Cells(1, 1) = x0
Cells(1, 2) = y0
Cells(1, 3) = z0
Cells(2, 1) = x1
Cells(2, 2) = y1
Cells(2, 3) = z1
Cells(3, 1) = x2
Cells(3, 2) = y2
Cells(3, 3) = z2
Cells(4, 1) = x3
Cells(4, 2) = y3
Cells(4, 3) = z3
'-----------------------------------------------------------------------------------------------------
irow = 5
Cells(irow + 0, 4) = r0
Cells(irow + 1, 4) = r1
Cells(irow + 2, 4) = r2
Cells(irow + 3, 4) = r3
myXYZ = myOffset3PlaneIntersection(irow + 0, x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3, r1, r2, r3)
myXYZ = myOffset3PlaneIntersection(irow + 1, x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, z0, r2, r3, r0)
myXYZ = myOffset3PlaneIntersection(irow + 2, x2, y2, z2, x3, y3, z3, x0, y0, z0, x1, y1, z1, r3, r0, r1)
myXYZ = myOffset3PlaneIntersection(irow + 3, x3, y3, z3, x0, y0, z0, x1, y1, z1, x2, y2, z2, r0, r1, r2)
myFormat
End Sub
'1 0.000 0.000 2.000
'2 0.000 0.000 0.000
'3 1.000 0.000 0.000
'4 0.000 1.000 0.000
'5 0.250 0.250 0.250 0.250 0.000
'6 0.000 0.000 0.000 0.250 0.000 ???? I want 0.250 0.250 0.250
'7 0.250 0.250 0.250 0.250 0.000
'8 0.102 0.339 0.102 0.250 0.000 ???? I want 0.250 0.250 0.250
(2021-11-15)我试试sympy
from sympy import *
def myVol(PTO,PTA,PTB,PTC):
return Matrix([[PTA.x-PTO.x, PTA.y-PTO.y, PTA.z-PTO.z], [PTB.x-PTO.x, PTB.y-PTO.y, PTB.z-PTO.z], [PTC.x-PTO.x, PTC.y-PTO.y, PTC.z-PTO.z]]).det()/6
def myUnitVector(myPoint3D):
myL=myPoint3D.distance((0, 0))
return Point3D(myPoint3D.x/myL,myPoint3D.y/myL,myPoint3D.z/myL)
def myHtoP(myHairetu):
return Point3D(myHairetu[0],myHairetu[1],myHairetu[2])
def my3PLaneIntersection(PTO,PTA,PTB,PTC,RA,RB,RC):
vA =myUnitVector(myHtoP(Plane(PTO, PTB, PTC).normal_vector))
PLA=Plane(PTO + RA * vA, normal_vector=vA)
vB =myUnitVector(myHtoP(Plane(PTO, PTC, PTA).normal_vector))
PLB=Plane(PTO + RB * vB, normal_vector=vB)
vC =myUnitVector(myHtoP(Plane(PTO, PTA, PTB).normal_vector))
PLC=Plane(PTO + RC * vC, normal_vector=vC)
return PLC.intersection(PLB.intersection(PLA)[0])
PTO,PTA,PTB,PTC=Point3D(0,0,0),Point3D(1,0,0),Point3D(0,1,0),Point3D(0,0,2)
print("#",myVol(PTO,PTA,PTB,PTC))
myRO,myRA,myRB,myRC=0,0,0,0
print("#",myVol(
my3PLaneIntersection(PTO,PTA,PTB,PTC,myRO,myRO,myRO)[0],
my3PLaneIntersection(PTA,PTC,PTB,PTO,myRA,myRA,myRA)[0],
my3PLaneIntersection(PTB,PTC,PTO,PTA,myRB,myRB,myRB)[0],
my3PLaneIntersection(PTC,PTB,PTA,PTO,myRC,myRC,myRC)[0]))
myRO,myRA,myRB,myRC=1/4,1/4,1/4,1/4
print("#",myVol(
my3PLaneIntersection(PTO,PTA,PTB,PTC,myRO,myRO,myRO)[0],
my3PLaneIntersection(PTA,PTC,PTB,PTO,myRA,myRA,myRA)[0],
my3PLaneIntersection(PTB,PTC,PTO,PTA,myRB,myRB,myRB)[0],
my3PLaneIntersection(PTC,PTB,PTA,PTO,myRC,myRC,myRC)[0]))
myRO,myRA,myRB,myRC=-1/4,-1/4,-1/4,-1/4
print("#",myVol(
my3PLaneIntersection(PTO,PTA,PTB,PTC,myRO,myRO,myRO)[0],
my3PLaneIntersection(PTA,PTC,PTB,PTO,myRA,myRA,myRA)[0],
my3PLaneIntersection(PTB,PTC,PTO,PTA,myRB,myRB,myRB)[0],
my3PLaneIntersection(PTC,PTB,PTA,PTO,myRC,myRC,myRC)[0]))
# 1/3
# 1/3
# 0
# 8/3
我检查
from sympy import *
var('m')
def myUnitVector(myPoint3D):
myL=myPoint3D.distance((0, 0))
return Point3D(myPoint3D.x/myL,myPoint3D.y/myL,myPoint3D.z/myL)
def myHtoP(myHairetu):
return Point3D(myHairetu[0],myHairetu[1],myHairetu[2])
def my3PLaneIntersection(PTO,PTA,PTB,PTC,RA,RB,RC):
vA =myUnitVector(myHtoP(Plane(PTO, PTB, PTC).normal_vector))
PLA=Plane(PTO + RA * vA, normal_vector=vA)
vB =myUnitVector(myHtoP(Plane(PTO, PTC, PTA).normal_vector))
PLB=Plane(PTO + RB * vB, normal_vector=vB)
vC =myUnitVector(myHtoP(Plane(PTO, PTA, PTB).normal_vector))
PLC=Plane(PTO + RC * vC, normal_vector=vC)
return PLC.intersection(PLB.intersection(PLA)[0])
def myProjection(PTO,PTA,PTC,PTOO):
v = myUnitVector(myHtoP(Plane(PTO, PTA, PTC).normal_vector))
return solve(PTOO - Plane(PTO, PTA, PTC).projection(PTOO) - m * v, m)[m] * (-1)
PTO,PTA,PTB,PTC=Point3D(0,0,0),Point3D(1,0,0),Point3D(0,1,0),Point3D(0,0,2)
myRO,myRA,myRB,myRC=-1/4,-1/4,-1/4,-1/4
PTOO= my3PLaneIntersection(PTO, PTA,PTB,PTC,myRO,myRO,myRO)[0]
PTAA= my3PLaneIntersection(PTA, PTC,PTB,PTO,myRA,myRA,myRA)[0]
PTBB= my3PLaneIntersection(PTB, PTC,PTO,PTA,myRB,myRB,myRB)[0]
PTCC= my3PLaneIntersection(PTC, PTB,PTA,PTO,myRC,myRC,myRC)[0]
print("#",PTOO,"\n#",PTAA,"\n#",PTBB,"\n#",PTCC,)
print("#",myProjection(PTO,PTB,PTA,PTOO),
myProjection(PTO,PTA,PTC,PTOO),
myProjection(PTO,PTC,PTB,PTOO))
print("#",myProjection(PTA,PTO,PTB,PTAA),
myProjection(PTA,PTC,PTO,PTAA),
myProjection(PTA,PTB,PTC,PTAA))
print("#",myProjection(PTB,PTA,PTO,PTBB),
myProjection(PTB,PTO,PTC,PTBB),
myProjection(PTB,PTC,PTA,PTBB))
print("#",myProjection(PTC,PTO,PTA,PTCC),
myProjection(PTC,PTB,PTO,PTCC),
myProjection(PTC,PTA,PTB,PTCC))
# Point3D(-1/4, -1/4, -1/4)
# Point3D(7/4, -1/4, -1/4)
# Point3D(-1/4, 7/4, -1/4)
# Point3D(-1/4, -1/4, 15/4)
# -1/4 -1/4 -1/4
# -1/4 -1/4 -1/4
# -1/4 -1/4 -1/4
# -1/4 -1/4 -1/4
triangular pyramid, trigonal pyramid offset-3d-Plane-Intersection By vba求解器
x=0.250 y=0.250 z=0.250‖内切球‖中心
每个平面偏移0.25 我挑战vba求解器
为什么全行 0.250 0.250 0.250
我要全行 0.250 0.250 0.250
为什么6行8行错误
结果很奇怪
请告诉我如何解决它。
Const x0 = 0
Const y0 = 0
Const z0 = 2
Const x1 = 0
Const y1 = 0
Const z1 = 0
Const x2 = 1
Const y2 = 0
Const z2 = 0
Const x3 = 0
Const y3 = 1
Const z3 = 0
Const r0 = 0.25
Const r1 = 0.25
Const r2 = 0.25
Const r3 = 0.25
Function myR1C1toA1(i, j)
myR1C1toA1 = Application.ConvertFormula("R" & i & "C" & j, xlR1C1, xlA1)
End Function
Function myPlane(Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz)
a = (By - Ay) * (Cz - Az) - (Cy - Ay) * (Bz - Az)
b = (Bz - Az) * (Cx - Ax) - (Cz - Az) * (Bx - Ax)
c = (Bx - Ax) * (Cy - Ay) - (Cx - Ax) * (By - Ay)
d = -(a * Ax + b * Ay + c * Az)
myPlane = Array(a, b, c, d)
End Function
Function myOffset3PlaneIntersection(irow, x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3, r1, r2, r3)
Dim v1 As Variant
Dim v2 As Variant
Dim v3 As Variant
'
v1 = myPlane(x0, y0, z0, x2, y2, z2, x3, y3, z3)
v2 = myPlane(x0, y0, z0, x3, y3, z3, x1, y1, z1)
v3 = myPlane(x0, y0, z0, x1, y1, z1, x2, y2, z2)
'
my1 = "(Abs(" & v1(0) & " * " & myR1C1toA1(irow, 1) & "+ " & v1(1) & " * " & myR1C1toA1(irow, 2) & "+ " & v1(2) & " * " & myR1C1toA1(irow, 3) & "+ " & v1(3) & ") / Sqrt(" & v1(0) & " ^ 2 + " & v1(1) & " ^ 2 +" & v1(2) & " ^ 2) - " & myR1C1toA1(irow, 1) & ")"
my2 = "(Abs(" & v2(0) & " * " & myR1C1toA1(irow, 1) & "+ " & v2(1) & " * " & myR1C1toA1(irow, 2) & "+ " & v2(2) & " * " & myR1C1toA1(irow, 3) & "+ " & v2(3) & ") / Sqrt(" & v2(0) & " ^ 2 + " & v2(1) & " ^ 2 +" & v2(2) & " ^ 2) - " & myR1C1toA1(irow, 2) & ")"
my3 = "(Abs(" & v3(0) & " * " & myR1C1toA1(irow, 1) & "+ " & v3(1) & " * " & myR1C1toA1(irow, 2) & "+ " & v3(2) & " * " & myR1C1toA1(irow, 3) & "+ " & v3(3) & ") / Sqrt(" & v3(0) & " ^ 2 + " & v3(1) & " ^ 2 +" & v3(2) & " ^ 2) - " & myR1C1toA1(irow, 3) & ")"
Range(myR1C1toA1(irow, 5)).Formula = "=" & my1 & "^ 2 +" & my2 & "^ 2 +" & my3 & "^ 2"
'
Dim ws As Worksheet: Set ws = ActiveSheet
SolverReset
SolverOk setCell:=ws.Range(myR1C1toA1(irow, 5)), _
MaxMinVal:=3, _
ByChange:=ws.Range(myR1C1toA1(irow, 1) & ":" & myR1C1toA1(irow, 3)), _
EngineDesc:="GRG Nonlinear"
SolverSolve UserFinish:=True
End Function
Sub myFormat()
Columns("A:E").Select
Selection.NumberFormatLocal = "0.000_ "
Range("A5").Select
End Sub
Sub aaa_Main()
Dim myXYZ As Variant
ActiveSheet.Cells.Clear
MsgBox "Solver Start"
Cells(1, 1) = x0
Cells(1, 2) = y0
Cells(1, 3) = z0
Cells(2, 1) = x1
Cells(2, 2) = y1
Cells(2, 3) = z1
Cells(3, 1) = x2
Cells(3, 2) = y2
Cells(3, 3) = z2
Cells(4, 1) = x3
Cells(4, 2) = y3
Cells(4, 3) = z3
'-----------------------------------------------------------------------------------------------------
irow = 5
Cells(irow + 0, 4) = r0
Cells(irow + 1, 4) = r1
Cells(irow + 2, 4) = r2
Cells(irow + 3, 4) = r3
myXYZ = myOffset3PlaneIntersection(irow + 0, x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3, r1, r2, r3)
myXYZ = myOffset3PlaneIntersection(irow + 1, x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, z0, r2, r3, r0)
myXYZ = myOffset3PlaneIntersection(irow + 2, x2, y2, z2, x3, y3, z3, x0, y0, z0, x1, y1, z1, r3, r0, r1)
myXYZ = myOffset3PlaneIntersection(irow + 3, x3, y3, z3, x0, y0, z0, x1, y1, z1, x2, y2, z2, r0, r1, r2)
myFormat
End Sub
'1 0.000 0.000 2.000
'2 0.000 0.000 0.000
'3 1.000 0.000 0.000
'4 0.000 1.000 0.000
'5 0.250 0.250 0.250 0.250 0.000
'6 0.000 0.000 0.000 0.250 0.000 ???? I want 0.250 0.250 0.250
'7 0.250 0.250 0.250 0.250 0.000
'8 0.102 0.339 0.102 0.250 0.000 ???? I want 0.250 0.250 0.250
(2021-11-15)我试试sympy
from sympy import *
def myVol(PTO,PTA,PTB,PTC):
return Matrix([[PTA.x-PTO.x, PTA.y-PTO.y, PTA.z-PTO.z], [PTB.x-PTO.x, PTB.y-PTO.y, PTB.z-PTO.z], [PTC.x-PTO.x, PTC.y-PTO.y, PTC.z-PTO.z]]).det()/6
def myUnitVector(myPoint3D):
myL=myPoint3D.distance((0, 0))
return Point3D(myPoint3D.x/myL,myPoint3D.y/myL,myPoint3D.z/myL)
def myHtoP(myHairetu):
return Point3D(myHairetu[0],myHairetu[1],myHairetu[2])
def my3PLaneIntersection(PTO,PTA,PTB,PTC,RA,RB,RC):
vA =myUnitVector(myHtoP(Plane(PTO, PTB, PTC).normal_vector))
PLA=Plane(PTO + RA * vA, normal_vector=vA)
vB =myUnitVector(myHtoP(Plane(PTO, PTC, PTA).normal_vector))
PLB=Plane(PTO + RB * vB, normal_vector=vB)
vC =myUnitVector(myHtoP(Plane(PTO, PTA, PTB).normal_vector))
PLC=Plane(PTO + RC * vC, normal_vector=vC)
return PLC.intersection(PLB.intersection(PLA)[0])
PTO,PTA,PTB,PTC=Point3D(0,0,0),Point3D(1,0,0),Point3D(0,1,0),Point3D(0,0,2)
print("#",myVol(PTO,PTA,PTB,PTC))
myRO,myRA,myRB,myRC=0,0,0,0
print("#",myVol(
my3PLaneIntersection(PTO,PTA,PTB,PTC,myRO,myRO,myRO)[0],
my3PLaneIntersection(PTA,PTC,PTB,PTO,myRA,myRA,myRA)[0],
my3PLaneIntersection(PTB,PTC,PTO,PTA,myRB,myRB,myRB)[0],
my3PLaneIntersection(PTC,PTB,PTA,PTO,myRC,myRC,myRC)[0]))
myRO,myRA,myRB,myRC=1/4,1/4,1/4,1/4
print("#",myVol(
my3PLaneIntersection(PTO,PTA,PTB,PTC,myRO,myRO,myRO)[0],
my3PLaneIntersection(PTA,PTC,PTB,PTO,myRA,myRA,myRA)[0],
my3PLaneIntersection(PTB,PTC,PTO,PTA,myRB,myRB,myRB)[0],
my3PLaneIntersection(PTC,PTB,PTA,PTO,myRC,myRC,myRC)[0]))
myRO,myRA,myRB,myRC=-1/4,-1/4,-1/4,-1/4
print("#",myVol(
my3PLaneIntersection(PTO,PTA,PTB,PTC,myRO,myRO,myRO)[0],
my3PLaneIntersection(PTA,PTC,PTB,PTO,myRA,myRA,myRA)[0],
my3PLaneIntersection(PTB,PTC,PTO,PTA,myRB,myRB,myRB)[0],
my3PLaneIntersection(PTC,PTB,PTA,PTO,myRC,myRC,myRC)[0]))
# 1/3
# 1/3
# 0
# 8/3
我检查
from sympy import *
var('m')
def myUnitVector(myPoint3D):
myL=myPoint3D.distance((0, 0))
return Point3D(myPoint3D.x/myL,myPoint3D.y/myL,myPoint3D.z/myL)
def myHtoP(myHairetu):
return Point3D(myHairetu[0],myHairetu[1],myHairetu[2])
def my3PLaneIntersection(PTO,PTA,PTB,PTC,RA,RB,RC):
vA =myUnitVector(myHtoP(Plane(PTO, PTB, PTC).normal_vector))
PLA=Plane(PTO + RA * vA, normal_vector=vA)
vB =myUnitVector(myHtoP(Plane(PTO, PTC, PTA).normal_vector))
PLB=Plane(PTO + RB * vB, normal_vector=vB)
vC =myUnitVector(myHtoP(Plane(PTO, PTA, PTB).normal_vector))
PLC=Plane(PTO + RC * vC, normal_vector=vC)
return PLC.intersection(PLB.intersection(PLA)[0])
def myProjection(PTO,PTA,PTC,PTOO):
v = myUnitVector(myHtoP(Plane(PTO, PTA, PTC).normal_vector))
return solve(PTOO - Plane(PTO, PTA, PTC).projection(PTOO) - m * v, m)[m] * (-1)
PTO,PTA,PTB,PTC=Point3D(0,0,0),Point3D(1,0,0),Point3D(0,1,0),Point3D(0,0,2)
myRO,myRA,myRB,myRC=-1/4,-1/4,-1/4,-1/4
PTOO= my3PLaneIntersection(PTO, PTA,PTB,PTC,myRO,myRO,myRO)[0]
PTAA= my3PLaneIntersection(PTA, PTC,PTB,PTO,myRA,myRA,myRA)[0]
PTBB= my3PLaneIntersection(PTB, PTC,PTO,PTA,myRB,myRB,myRB)[0]
PTCC= my3PLaneIntersection(PTC, PTB,PTA,PTO,myRC,myRC,myRC)[0]
print("#",PTOO,"\n#",PTAA,"\n#",PTBB,"\n#",PTCC,)
print("#",myProjection(PTO,PTB,PTA,PTOO),
myProjection(PTO,PTA,PTC,PTOO),
myProjection(PTO,PTC,PTB,PTOO))
print("#",myProjection(PTA,PTO,PTB,PTAA),
myProjection(PTA,PTC,PTO,PTAA),
myProjection(PTA,PTB,PTC,PTAA))
print("#",myProjection(PTB,PTA,PTO,PTBB),
myProjection(PTB,PTO,PTC,PTBB),
myProjection(PTB,PTC,PTA,PTBB))
print("#",myProjection(PTC,PTO,PTA,PTCC),
myProjection(PTC,PTB,PTO,PTCC),
myProjection(PTC,PTA,PTB,PTCC))
# Point3D(-1/4, -1/4, -1/4)
# Point3D(7/4, -1/4, -1/4)
# Point3D(-1/4, 7/4, -1/4)
# Point3D(-1/4, -1/4, 15/4)
# -1/4 -1/4 -1/4
# -1/4 -1/4 -1/4
# -1/4 -1/4 -1/4
# -1/4 -1/4 -1/4