(三角錐の内接球)「2021年度 九州大・理系数学問題1(1)」」をsympyでやってみたい。

オリジナルポスト
(0,0,0),(1,0,0),(0,1,0),(0,0,2)
「2021年度 九州大・理系数学問題1(1)」

https://www.densu.jp/kyusyu/21kyusyuspass.pdf

上のpdfのページより「解答解説のページへ」いけます。


#sympyで

内接球の半径は、体積と表面積から計算しました。3平面の交点計算です。
①コードを短くする方法を教えて下さい。
②右手系、左手系の入力がわかっていません。どちらでも対応する予定です。
➂別解(三角形の断面？2等分線？)にチャレンジしたいので、アドバイスをお願いします。
④目標は、平行三角錐？の体積計算です。またまた、アドバイスをお願いします。

```python
from sympy import *
def myDis(x1,y1,z1,x2,y2,z2):
    return sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1 - z2) ** 2)
def myHeron(x1, y1, z1, x2, y2, z2, x3, y3, z3):
    s = (myDis(x1, y1, z1, x2, y2, z2) + myDis(x2, y2, z2, x3, y3, z3) + myDis(x3, y3, z3, x1, y1, z1)) / 2
    return sqrt(s * (s - myDis(x1, y1, z1, x2, y2, z2))
                  * (s - myDis(x2, y2, z2, x3, y3, z3))
                  * (s - myDis(x3, y3, z3, x1, y1, z1)))
def mySimentaiMenseki(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4):
    return myHeron(x1, y1, z1, x2, y2, z2, x3, y3, z3) \
    +      myHeron(x2, y2, z2, x3, y3, z3, x4, y4, z4) \
    +      myHeron(x3, y3, z3, x4, y4, z4, x1, y1, z1) \
    +      myHeron(x4, y4, z4, x1, y1, z1, x2, y2, z2)
def mySimentaiTaiseki(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4):
    return ((x4 - x1) * abs(y2 - y1) * (z3 - z1) - (z3 - z1) * (y3 - y1)
    +       (y4 - y1) * abs(z2 - z1) * (x3 - x1) - (x2 - x1) * (z3 - z1)
    +       (z4 - z1) * abs(x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)) / 6
def myHankei(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4):
    return mySimentaiTaiseki(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4)*3\
    /      mySimentaiMenseki(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4)
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 my3PaneIntersection(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=Point3D(0,0,0)
PTA=Point3D(1,0,0)
PTB=Point3D(0,1,0)
PTC=Point3D(0,0,2)
myR=myHankei(PTO.x,PTO.y,PTO.z,
             PTA.x,PTA.y,PTA.z,
             PTB.x,PTB.y,PTB.z,
             PTC.x,PTC.y,PTC.z)
print("#",my3PaneIntersection(PTO,PTA,PTB,PTC,myR,myR,myR))
# [Point3D(1/4, 1/4, 1/4)]
```

(2021-11-15)ちょっと直しました。
①ヘロンの公式(github参考)
②三角錘の体積。行列式使用。

https://gist.github.com/wootfish/74114158a367cfd8356b

```python
from sympy import *
# def myTaiseki(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()
def myTaiseki(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 myMenseki(PTO,PTA,PTB,PTC):
    return  herons(PTO.distance(PTA),PTA.distance(PTB),PTB.distance(PTO))\
           +herons(PTO.distance(PTA),PTA.distance(PTC),PTC.distance(PTO))\
           +herons(PTA.distance(PTC),PTC.distance(PTB),PTB.distance(PTA))\
           +herons(PTB.distance(PTO),PTO.distance(PTC),PTC.distance(PTB))
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=Point3D(0,0,0)
PTA=Point3D(1,0,0)
PTB=Point3D(0,1,0)
PTC=Point3D(0,0,2)
myR=myTaiseki(PTO,PTA,PTB,PTC)/myMenseki(PTO,PTA,PTB,PTC)*3
print("#",my3PLaneIntersection(PTO,PTA,PTB,PTC,myR,myR,myR))
# [Point3D(1/4, 1/4, 1/4)]
```


#参考

https://manabitimes.jp/math/665

https://qiita.com/mrrclb48z/items/3f17157fd2a98a85e2a0

https://ja.stackoverflow.com/questions/83009/%e5%9b%9b%e9%9d%a2%e4%bd%93%e4%b8%89%e8%a7%92%e9%8c%90%e3%81%ae%e5%86%85%e6%8e%a5%e7%90%83%e3%81%aesympy%e3%82%92%e6%8e%a2%e3%81%97%e3%81%a6%e3%81%84%e3%81%be%e3%81%99-vba%e3%82%bd%e3%83%ab%e3%83%90%e3%83%bc%e3%81%a7%e3%82%84%e3%81%a3%e3%81%a6%e3%81%bf%e3%81%be%e3%81%97%e3%81%9f/83097#83097

https://qiita.com/mrrclb48z/items/159ea13c0fb80b3a8f0d