四面体「2023群馬大学前期【３】」をsympyとFreeCADでやってみたい。できませんでした。

できませんでした。
FreeCADのマクロで作図しました。

オリジナル 公式ホームページより
医(医学科)学部 理工,情報学部【３】の類題
問題文と解答例

https://www.gunma-u.ac.jp/wp-content/uploads/2023/03/4c2753bf0b954793918a8209864418e9.pdf#page=3

一般選抜（前期日程）、私費外国人留学生選抜＜群馬大学

https://www.gunma-u.ac.jp/admission/adm001/adm001_002/2023%e5%b9%b4%e5%ba%a6/2023%e5%b9%b4%e5%ba%a6%e3%80%80%e4%b8%80%e8%88%ac%e9%81%b8%e6%8a%9c%ef%bc%88%e5%89%8d%e6%9c%9f%e6%97%a5%e7%a8%8b%ef%bc%89%e3%80%81%e7%a7%81%e8%b2%bb%e5%a4%96%e5%9b%bd%e4%ba%ba%e7%95%99%e5%ad%a6

上と同じです。大学入試数学問題集成＞【3】テキスト

https://mathexamtest.web.fc2.com/2023/202310201/2023102010100mj.html#top-0111

# sympyで(公式ページの解答例を参考に)
勉強のしなおしです。

# sympyで(sympy的?安易なやり方)
できませんでした。座標計算です。
・Pointクラスの値の変更がわからなかったので、3回繰り返しました。
　(1),(2),(3)に分けました。
・「同一平面上は...」、点から平面への垂線(projection投影)です。
・絶対値のはずしかた？意味を教えて下さい。

```python
from sympy import *
var('x y z')
ax,ay   =symbols("ax ay"   ,real=True)
cx,cy   =symbols("cx cy"   ,real=True)
dx,dy   =symbols("dx dy "  ,real=True)
dz      =symbols(      "dz",real=True,positive=True)
u,v,w   =symbols("u v w"   ,real=True)
s,t     =symbols("s t"     ,real=True,positive=True)
L       =symbols("L"       ,real=True,positive=True)
def myMenseki3D(O,A,B):
    P=A-O
    Q=B-O
    return Rational(1, 2) * sqrt(P.distance(Point(0, 0, 0)) ** 2 * Q.distance(Point(0, 0, 0)) ** 2 \
                                 - P.dot(Q) ** 2).simplify()
def mySphereFormula(myP,myR):
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
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 myVectorABC(OQ, OA,OB,OC):
    return solve([
                  Eq(OQ.x,u*OA.x+v*OB.x+w*OC.x),
                  Eq(OQ.y,u*OA.y+v*OB.y+w*OC.y),
                  Eq(OQ.z,u*OA.z+v*OB.z+w*OC.z)
                 ],
                 [u,v,w])
###########################################
O=Point(0 ,0 ,0 )
A=Point(ax,ay,0 )
C=Point(cx,cy,0 )
B=A+C
D=Point(dx,dy,dz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
Y=Plane(A,C,D).intersection(Line(O,X))[0]
print("#(1)",myVectorABC(Y, A,C,D))
########################################################
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
print("#(2)",solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0] )
########################################################
L=1
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
DPQ=myMenseki3D(D,P,Q).subs({t:ft})
print("#(3)XXXX",DPQ )
# 間違っています。絶対値をはずしました。
DPQabs=(s - 1)*((3*s - 2)/(4*s - 3) - 1)/2
ans=solve(diff(DPQabs,s))
print("#(3)",ans[0],DPQabs.subs({s:ans[0]}) )
plot(DPQ,(s,solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0]  ))
#(1) {u: 1/4, v: 1/4, w: 1/2}
#(2) 2/3 1
#(3)XXXX Abs((s - 1)*((3*s - 2)/(4*s - 3) - 1))/2
#(3) 1/2 1/8
```
①私の環境は,pycharmです。
②よく聞くのは、Jupyterです。
③web上で、上記のソースを「SymPy Live shell」に、コピー貼り付けでもできました。
黒背景の右上に、マウスを移動すると、コピーマークが発生します。

https://live.sympy.org/

1/3<=s<=1の△DPQの面積のグラフです。
plotコマンド。縦横比が異なります。

![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/75393349-5a8f-c012-0a85-708cfd3837b1.png)


# FreeCADのマクロで作図

```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('x y z')
ax,ay   =symbols("ax ay"   ,real=True)
cx,cy   =symbols("cx cy"   ,real=True)
dx,dy   =symbols("dx dy "  ,real=True)
dz      =symbols(      "dz",real=True,positive=True)
u,v,w   =symbols("u v w"   ,real=True)
s,t     =symbols("s t"     ,real=True,positive=True)
L       =symbols("L"       ,real=True,positive=True)
def myMenseki3D(O,A,B):
    P=A-O
    Q=B-O
    return Rational(1, 2) * sqrt(P.distance(Point(0, 0, 0)) ** 2 * Q.distance(Point(0, 0, 0)) ** 2 \
                                 - P.dot(Q) ** 2).simplify()
def mySphereFormula(myP,myR):
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
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 myVectorABC(OQ, OA,OB,OC):
    return solve([
                  Eq(OQ.x,u*OA.x+v*OB.x+w*OC.x),
                  Eq(OQ.y,u*OA.y+v*OB.y+w*OC.y),
                  Eq(OQ.z,u*OA.z+v*OB.z+w*OC.z)
                 ],
                 [u,v,w])
###########################################
O=Point(0 ,0 ,0 )
A=Point(ax,ay,0 )
C=Point(cx,cy,0 )
B=A+C
D=Point(dx,dy,dz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
Y=Plane(A,C,D).intersection(Line(O,X))[0]
print("#(1)",myVectorABC(Y, A,C,D))
########################################################
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
print("#(2)",solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0] )
########################################################
L=1
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
DPQ=myMenseki3D(D,P,Q).subs({t:ft})
print("#(3)XXXX",DPQ )
# wrong
DPQabs=(s - 1)*((3*s - 2)/(4*s - 3) - 1)/2
ans=solve(diff(DPQabs,s))
print("#(3)",ans[0],DPQabs.subs({s:ans[0]}) )
#plot(DPQ,(s,solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0]  ))
##########################################################################################################
### 作図用
# wrong
P=P.subs({s:1/2})
Q=Q.subs({t:1/2})
##########################################################################################################
### 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.05 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
def myLine_C(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    myLine(args[i],args[0])
    return
def myLine_H(*args):
    for i in range(1,len(args)):
        myLine(args[0],args[i])
    return
#
print("X",X)
print("P",P)
print("Q",Q)
#
myLine_C  (O,A,B,C)
myLine_H  (D, O,A,B,C)
myLine_C  (O,P,X,Q)
myLine_C  (D,P,Q)
myTxtXYZ_S(O,"O",A,"A",B,"B",C,"C",D,"D")
myTxtXYZ_S(X,"X",P,"P",Q,"Q")
#
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方向です。
Originは、非表示です。マウス操作で、線分と文字を赤色にしました。

![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/3840a71c-62d2-78af-580f-275dbfc44768.png)

上図の拡大図です。

![1png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/70753471-b2a5-54a8-25d1-5b177ba23816.png)

# 参考

https://qiita.com/mrrclb48z/items/881a72b88333798afbc1

https://stackoverflow.com/questions/70800934/sympy-subs-cannot-substituting-a-point-to-array








https://chokkan.github.io/python/13class.html