空間図形 三角錐の体積「2024 慶應義塾大学 経済[4]」をsympyとFreeCADでやってみたい。

・３次元FreeCADのワイヤーフレームです。
・勉強中。軸線方向の見え方。文字の回転です。(めぐろ塾 様）
オリジナル
問題
https://www.yomiuri.co.jp/nyushi/sokuho/s_mondaitokaitou/keio/mondai/img/keio_213_sugaku_mon.pdf#page=4
解答
https://www.yomiuri.co.jp/nyushi/sokuho/s_mondaitokaitou/keio/kaitou/img/keio_213_sugaku_kai.pdf#page=2

 ＜ 2024 慶應義塾大学 経済 ＜ 読売新聞オンライン 様　ページ
https://www.yomiuri.co.jp/nyushi/sokuho/s_mondaitokaitou/keio/1376715_5113

上と同じです。大学入試数学問題集成 様＞テキスト
未登録(2024/02/23)

公式ホームページ
未登録(2024/02/23)

# 解答解説
めぐろ塾 様
操作はスクロールして下さい。ジャンプしません。[4]

https://www.megurojyuku.com/2024%e6%85%b6%e6%87%89%e7%b5%8c%e6%b8%88%e3%80%90%e6%95%b0%e5%ad%a6%e3%80%91%e8%a7%a3%e7%ad%94%e9%80%9f%e5%a0%b1/

# sympyで
・ver0.2
・mapを使いました。わかりやすいかどうかです。
・2点からの連立方程式、solve関数です。
https://docs.sympy.org/latest/modules/solvers/solvers.html#sympy.solvers.solvers.solve

2024-02-27 追加
・良い子は、マネをしないで下さい。
　4点の同一平面上を、三角錐の体積0でやっています。辺の長さ0を考えていません。
　三角形２つの面積が、nonnegativeだからいいような気もしますが？
```python
# ver0.2
from sympy import *
var('x,y,z',real=True)
var('p,q'  ,real=True,nonnegative=True)
def myMensekiVector3D(P,Q,R):
    return Rational(1,2)*sqrt(P.distance(R)**2*Q.distance(R)** 2 \
           -(P-R).dot(Q-R)**2)
def myTaiseki(PTO,PTA,PTB,PTC):
    return abs(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 myP1P2P3(P1,P2):
    return P1+(P2-P1)*2 
O,A,B,C=map(Point,[(0,0,0),(3,-sqrt(3),0),(3,sqrt(3),0),(p,0,q)])                
L      =2*sqrt(3)
sol    =solve( [Eq(C.distance(O),L),Eq(C.distance(A),L)],[p,q])    
rep    ={p:sol[0][0],q:sol[0][1]}      
T      =Point(3/2,0,q/2)
C,T    =C.subs(rep),T.subs(rep)
D,E,F,G=myP1P2P3(O,T),myP1P2P3(A,T),myP1P2P3(B,T),myP1P2P3(C,T)                
H=Plane(A,B,C).intersection(Line(D,G))[0]  
I=Plane(D,E,G).intersection(Line(C,B))[0]  
J=Plane(D,F,G).intersection(Line(C,A))[0] 
print("#(1)",p,q)
print("#(2)",myMensekiVector3D(C,J,H)  +myMensekiVector3D(C,H,I)  ,myTaiseki(C, J,H,I))
print("#(3)",myTaiseki        (C,J,H,G)+myTaiseki        (C,H,I,G))
#(1) p q
#(2) sqrt(3)/6 0
#(3) 2*sqrt(6)/27
```
・ver0.1 上と同じです。
・３つの球の交点計算で、できませんでした。
・外心から三平方の定理です。(手計算の場合は、外心より重心がラクです。)
property circumcenter
https://docs.sympy.org/latest/modules/geometry/polygons.html#sympy.geometry.polygon.RegularPolygon.circumcenter
```python
# ver0.1
from sympy import *
var('x,y,z',real=True)
def myMensekiVector3D(P,Q,R):
    return Rational(1, 2) * sqrt(P.distance(R) ** 2 * Q.distance(R) ** 2 \
                                 - (P-R).dot(Q-R) ** 2)
def myTaiseki(PTO,PTA,PTB,PTC):
    return abs(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 mySphereFormula(myP,myR):
    print(myP)
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
O=Point(0,      0 ,0)                      #;print("# O",O)
A=Point(3,-sqrt(3),0)                      #;print("# A",A)
B=Point(3, sqrt(3),0)                      #;print("# B",B)
r=2*sqrt(3)
# C=Point(solve([mySphereFormula(O,r),               \
#                mySphereFormula(A,r),               \
#                mySphereFormula(B,r)],[x,y,z])
#         )
# print(C)
S=Triangle(O,A,B).circumcircle.args[0]
C=Point(S.x,S.y,sqrt(r**2-Point(S.x,S.y,0) .distance(O)**2))
p=C.x
q=C.z
print("#(1)",p,q)
T=Point(3/2,0,q/2)
D=O+(T-O)*2                                #;print("# D",D)
E=A+(T-A)*2                                #;print("# E",E)s
F=B+(T-B)*2                                #;print("# F",F)
G=C+(T-C)*2                                #;print("# G",G)
H=Plane(A,B,C).intersection(Line(D,G))[0]  #;print("# H",H)
I=Plane(D,E,G).intersection(Line(C,B))[0]  #;print("# I",I)
J=Plane(D,F,G).intersection(Line(C,A))[0]  #;print("# J",J)
print("#(2)",myMensekiVector3D(C,J,H)  +myMensekiVector3D(C,H,I)  ,myTaiseki(C, J,H,I))
print("#(3)",myTaiseki        (C,J,H,G)+myTaiseki        (C,H,I,G))
#(1) 2 2*sqrt(2)
#(2) sqrt(3)/6 0
#(3) 2*sqrt(6)/27
```

# FreeCADのマクロで作図
・計算部分は、ver0.2をコピー貼り付けです。
```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
# ver0.2
from sympy import *
var('x,y,z',real=True)
var('p,q'  ,real=True,nonnegative=True)
def myMensekiVector3D(P,Q,R):
    return Rational(1,2)*sqrt(P.distance(R)**2*Q.distance(R)** 2 \
           -(P-R).dot(Q-R)**2)
def myTaiseki(PTO,PTA,PTB,PTC):
    return abs(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 myP1P2P3(P1,P2):
    return P1+(P2-P1)*2 
O,A,B,C=map(Point,[(0,0,0),(3,-sqrt(3),0),(3,sqrt(3),0),(p,0,q)])                
L      =2*sqrt(3)
sol    =solve( [Eq(C.distance(O),L),Eq(C.distance(A),L)],[p,q])    
rep    ={p:sol[0][0],q:sol[0][1]}      
T      =Point(3/2,0,q/2)
C,T    =C.subs(rep),T.subs(rep)
D,E,F,G=myP1P2P3(O,T),myP1P2P3(A,T),myP1P2P3(B,T),myP1P2P3(C,T)                
H=Plane(A,B,C).intersection(Line(D,G))[0]  
I=Plane(D,E,G).intersection(Line(C,B))[0]  
J=Plane(D,F,G).intersection(Line(C,A))[0] 
#print("#(1)",p,q)
#print("#(2)",myMensekiVector3D(C,J,H)  +myMensekiVector3D(C,H,I)  ,myTaiseki(C, J,H,I))
#print("#(3)",myTaiseki        (C,J,H,G)+myTaiseki        (C,H,I,G))
#(1) p q
#(2) sqrt(3)/6 0
#(3) 2*sqrt(6)/27
##########################################################################################################
### 作図用
##########################################################################################################
### 3D作図
def myXYZ2Txt(A):
    return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxtXYZ(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.15 mm'
    return
def myTxtXYZ_S(*xy_tx):
    for i in range(1,int(len(xy_tx)/2)+1):
        myTxtXYZ(xy_tx[2*i-2],xy_tx[2*i-1]+myXYZ2Txt(xy_tx[2*i-2]) )
    return
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)
    return
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return
#
myLine_S  (O,T)
myLine_S  (A,B,C,A)
myLine_S  (O,A)
myLine_S  (O,B)
myLine_S  (O,C)
myLine_S  (C,J,H,I,C)
myLine_S  (G,C)
myLine_S  (G,J)
myLine_S  (G,H)
myLine_S  (G,I)
myTxtXYZ_S(O,"O",A,"A",B,"B",C,"C")
myTxtXYZ_S(T,"T")
myTxtXYZ_S(D,"D",E,"E",F,"F",G,"G")
myTxtXYZ_S(H,"H",I,"I",J,"J")
#
doc = App.activeDocument()
App.ActiveDocument.addObject("App::Origin", "Origin")
App.ActiveDocument.getObject('Origin').Visibility = True
App.ActiveDocument.recompute()
Gui.activeDocument().activeView().viewAxonometric()
Gui.SendMsgToActiveView("ViewFit")
```
# 作図
・isometric方向?です。
![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/c6febbae-a473-337e-8a4e-60a1abf53145.png)

・拡大図
・Origin非表示はCAD操作です。
・着色はCAD操作です。
![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/d21690a7-3142-a48a-e8f7-a2e8deab9ba6.png)

# いつもの? sympyの実行環境と 参考のおすすめです。
[(テンプレート)](https://qiita.com/mrrclb48z/items/00dd08b0317069be9342)

# いつもと違うおすすめです。
#### 参考過去問

https://qiita.com/mrrclb48z/items/2bece4f967bcd3bf126a