直方体「2022年 早稲田本庄高等学院 数学[3]」をsympyとFreeCADでやってみた。

AB=3,AD=6,AE=6の直方体ABCD-EFGH....

オリジナル ますいしい様

~~https://ameblo.jp/mathisii/entry-12734081753.html~~

https://ameblo.jp/mathisii/entry-12767795972.html



# sympyで(ますいしい様の手順で)
［3］,［3］(別解) を勉強中

# sympyで(計算だけ)
・「同じ平面上」は四面体の体積計算(Matrix)で判別する方程式(solve)です。
・ヘロンの公式で空間三角形の面積計算です。外積の勉強中。
・座標系の向きは、FreeCAD用で別解と異なります。


```python
from sympy import *
var('a DJ')
def SimentaiTaiseki(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 Herons(a,b,c):
    s=(a+b+c)/2
    return sqrt(s*(s-a)*(s-b)*(s-c))
def Herons3(p,q,r):
    return Herons(p.distance(q), q.distance(r), r.distance(p))
def Herons4(p,q,r,s):
    return Herons3(p,q,s)+Herons3(s,q,r)
AB=3
AD=6
AE=6
A=Point( 0,  0,  0)
B=Point( 0,-AB,  0)
C=Point(AD,-AB,  0)
D=Point(AD,  0,  0)
E=Point( 0,  0,-AE)
F=Point( 0,-AB,-AE)
G=Point(AD,-AB,-AE)
H=Point(AD,  0,-AE)
I=B+Point(0, 0, -a)
J=D+Point(0, 0,-DJ)
DJ=solve(SimentaiTaiseki(A,I,G,J),DJ)[0]
print("#問1",DJ)
I=I.subs({a:3})
J=D-Point(0,0,DJ.subs({a:3}))
print("#問2",Herons4(A,I,G,J).simplify())
K=Rational(1,2)*(F+G)
L=Rational(1,2)*(G+H)
Plain_KL=Plane(K,L,K+Point(0,0,1))
P=Plain_KL.intersection(Line(G,I))[0]
Q=Plain_KL.intersection(Line(G,J))[0]
print("#問3", Rational(1,2)*AB*AD*AE \
             -Rational(1,3)*Herons4(K,P,Q,L).simplify()*Plain_KL.distance(G))
#問1 6 - a
#問2 27
#問3 207/4
```

# FreeCADのマクロで
上のsympy計算に、FreeCAD作図を追加しました。
```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('a DJ')
def SimentaiTaiseki(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 Herons(a,b,c):
    s=(a+b+c)/2
    return sqrt(s*(s-a)*(s-b)*(s-c))
def Herons3(p,q,r):
    return Herons(p.distance(q), q.distance(r), r.distance(p))
def Herons4(p,q,r,s):
    return Herons3(p,q,s)+Herons3(s,q,r)
AB=3
AD=6
AE=6
A=Point( 0,  0,  0)
B=Point( 0,-AB,  0)
C=Point(AD,-AB,  0)
D=Point(AD,  0,  0)
E=Point( 0,  0,-AE)
F=Point( 0,-AB,-AE)
G=Point(AD,-AB,-AE)
H=Point(AD,  0,-AE)
I=B+Point(0, 0, -a)
J=D+Point(0, 0,-DJ)
DJ=solve(SimentaiTaiseki(A,I,G,J),DJ)[0]
print("#問1",DJ)
I=I.subs({a:3})
J=D-Point(0,0,DJ.subs({a:3}))
print("#問2",Herons4(A,I,G,J).simplify())
K=Rational(1,2)*(F+G)
L=Rational(1,2)*(G+H)
Plain_KL=Plane(K,L,K+Point(0,0,1))
P=Plain_KL.intersection(Line(G,I))[0]
Q=Plain_KL.intersection(Line(G,J))[0]
print("#問3", Rational(1,2)*AB*AD*AE \
             -Rational(1,3)*Herons4(K,P,Q,L).simplify()*Plain_KL.distance(G))
#問1 6 - a
#問2 27
#問3 207/4
#########################################################################################################
# 3D作図
def myXYZ2Txt(A):
    return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxt_XYZ(A,myWedgei):
    P5x=float(A.x)
    P5y=float(A.y)
    P5z=float(A.z)
    p5 = FreeCAD.Vector(P5x, P5y, P5z)
    myText = Draft.makeText(myWedgei, p5)
    myText.Label = myWedgei
    FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.5 mm'
def myLine(A,B):
    Ax,Ay,Az=float(A.x),float(A.y),float(A.z)
    Bx,By,Bz=float(B.x),float(B.y),float(B.z)
    pl = FreeCAD.Placement()
    pl.Rotation.Q = (0.4247081540122249, 0.17592004639554645, 0.33985110062924484, 0.8204732460821097)
    pl.Base = FreeCAD.Vector(-3.9166066876399563, -2.1670824762243774, 1.7495260956243028)
    points = [FreeCAD.Vector(Ax,Ay,Az), FreeCAD.Vector(Bx,By,Bz)]
    line = Draft.make_wire(points, placement=pl, closed=False, face=True, support=None)
    Draft.autogroup(line)
myTxt_XYZ(A,"A"+myXYZ2Txt(A))
myTxt_XYZ(B,"B"+myXYZ2Txt(B))
myTxt_XYZ(C,"C"+myXYZ2Txt(C))
myTxt_XYZ(D,"D"+myXYZ2Txt(D))
myTxt_XYZ(E,"E"+myXYZ2Txt(E))
myTxt_XYZ(F,"F"+myXYZ2Txt(F))
myTxt_XYZ(G,"G"+myXYZ2Txt(G))
myTxt_XYZ(H,"H"+myXYZ2Txt(H))
myTxt_XYZ(I,"I"+myXYZ2Txt(I))
myTxt_XYZ(J,"J"+myXYZ2Txt(J))
myTxt_XYZ(K,"K"+myXYZ2Txt(K))
myTxt_XYZ(L,"L"+myXYZ2Txt(L))
myTxt_XYZ(P,"P"+myXYZ2Txt(P))
myTxt_XYZ(Q,"Q"+myXYZ2Txt(Q))
doc = App.activeDocument()
myLine(A,B)
myLine(B,C)
myLine(C,D)
myLine(D,A)
#
myLine(E,F)
myLine(F,G)
myLine(G,H)
myLine(H,E)
#
myLine(A,E)
myLine(B,F)
myLine(C,G)
myLine(D,H)
#
myLine(A,I)
myLine(I,G)
myLine(G,J)
myLine(J,A)
#
myLine(K,P)
myLine(P,Q)
myLine(Q,L)
#
App.ActiveDocument.addObject("App::Origin", "Origin")
App.ActiveDocument.getObject('Origin').Visibility = True
App.ActiveDocument.recompute()
Gui.activeDocument().activeView().viewAxonometric()
Gui.SendMsgToActiveView("ViewFit")
```
現在ワイヤーフレームで申し訳ありません。
FreeCADの操作をすれば、体積(float)Vを求める事ができる？はずです。

![png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/46ac4758-c643-b33d-021a-7292b549cd2b.png)


# 参考
四面体の体積計算(Matrix)＜

https://qiita.com/mrrclb48z/items/1dbdb218308190528cac