正八面体「2023広島大学前期【３】」をsympyとFreeCADでやってみた。

FreeCADのマクロで作図しました。
ワイヤーフレームです。ソリッドを勉強中。

オリジナル

https://youtu.be/xAiui8qBoJk

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

https://mathexamtest.web.fc2.com/2023/202310721/2023107210100mj.html#top-0103


解答　広島大学【３】＜2023年解答速報

https://kaisoku.kawai-juku.ac.jp/nyushi/honshi/23/hr1-22a.pdf#page=3

https://kamelink.com/exam/index.php?2023%C7%AF%B2%F2%C5%FA%C2%AE%CA%F3#w5951b80

# sympyで(オリジナルのやり方)
片山司朗様の方法を勉強中。

# sympyで(sympy的?安易なやり方)
・点E,Fは、solveで、3球面の交点。答えがでるとは思いませんでした。
・全部の点の座標を計算後、(1)(2)(3)
・(2)が1行でできなかったので、.xと.yと.zに分けました。
　他の大学の問題でも頻出するなら、関数を自作したらいいかも。
・今までならsympy微分diffです。今回微分を使わず、minimumです。不便でした。
　関数(合成？関数)を自作したらいいかも。
・sympyなので、最初にL=1をせずに、最後にsubs({L:1})でもできるかも。
・varは、複素数でした。正しくは,symbolのrealです。

```python
from sympy import *
var('x y z s t')
var('p q r')
t=Symbol('t',real=True,positive=True)
def mySphereFormula(myP,myR):
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
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()
L=1
a=Rational(1,2)*L
M=Point( 0, 0, 0)
A=Point(-a,-a, 0)
B=Point( a,-a, 0)
C=Point( a, a, 0)
D=Point(-a, a, 0)
Fxyz,Exyz=solve([mySphereFormula(A,L),
                 mySphereFormula(B,L),
                 mySphereFormula(C,L)],[x,y,z]
               )
E=Point(Exyz)
F=Point(Fxyz)
print("#(1)",(B-A).dot(D-A))
print("#(1)",(B-A).dot(E-A))
print("#(1)",(D-A).dot(E-A))
#できませんでした。
#pqr=solve([
#           Eq(F-A,p*(B-A)+q*(D-A)+r*(E-A))
#          ],
#          [p,q,r])
pqr=solve([
           Eq((F-A).x,p*(B-A).x+q*(D-A).x+r*(E-A).x),
           Eq((F-A).y,p*(B-A).y+q*(D-A).y+r*(E-A).y),
           Eq((F-A).z,p*(B-A).z+q*(D-A).z+r*(E-A).z)
          ],
          [p,q,r])
print("#(2)",pqr)
G=B+Rational(1,3)*(E-B)
H=C+           t *(F-C)
AGH=myMenseki3D(A,G,H)
print("#(2)",solveset(Eq(AGH,minimum(AGH,t,Interval(0,1))),t,Interval(0,1)).args[0],
                             minimum(AGH,t,Interval(0,1))
     )
#(1) 0
#(1) 1/2
#(1) 1/2
#(2) {p: 1, q: 1, r: -1}
#(2) 1/3 sqrt(42)/18
```
①私の環境は,pycharmです。
②よく聞くのは、Jupyterです。
③web上で、上記のソースを「SymPy Live shell」に、コピー貼り付けでもできました。
黒背景の右上に、マウスを移動すると、コピーマークが発生します。
計算時間は約1分です。資源の無駄使い。ごめんなさい。

https://live.sympy.org/

# FreeCADのマクロで作図

```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('x y z s t')
var('p q r')
t=Symbol('t',real=True,positive=True)
def mySphereFormula(myP,myR):
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
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()
L=1
a=Rational(1,2)*L
M=Point( 0, 0, 0)
A=Point(-a,-a, 0)
B=Point( a,-a, 0)
C=Point( a, a, 0)
D=Point(-a, a, 0)
Fxyz,Exyz=solve([mySphereFormula(A,L),
                 mySphereFormula(B,L),
                 mySphereFormula(C,L)],[x,y,z]
               )
E=Point(Exyz)
F=Point(Fxyz)
print("#(1)",(B-A).dot(D-A))
print("#(1)",(B-A).dot(E-A))
print("#(1)",(D-A).dot(E-A))
pqr=solve([
           Eq((F-A).x,p*(B-A).x+q*(D-A).x+r*(E-A).x),
           Eq((F-A).y,p*(B-A).y+q*(D-A).y+r*(E-A).y),
           Eq((F-A).z,p*(B-A).z+q*(D-A).z+r*(E-A).z)
          ],
          [p,q,r])
print("#(2)",pqr)
G=B+Rational(1,3)*(E-B)
H=C+           t *(F-C)
AGH=myMenseki3D(A,G,H)
print("#(2)",solveset(Eq(AGH,minimum(AGH,t,Interval(0,1))),t,Interval(0,1)).args[0],
                             minimum(AGH,t,Interval(0,1))
     )
##########################################################################################################
### 作図用
H=H.subs({t:solveset(Eq(AGH,minimum(AGH,t,Interval(0,1))),t,Interval(0,1)).args[0]})
##########################################################################################################
### 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
#
myLine_C  (A,B,C,D)
myLine_H  (E, A,B,C,D)
myLine_H  (F, A,B,C,D)
myLine    (A,C)
myLine    (B,D)
myLine_C  (A,G,H)
#
myTxtXYZ_S(M,"M")
myTxtXYZ_S(A,"A",B,"B",C,"C",D,"D")
myTxtXYZ_S(E,"E",F,"F")
myTxtXYZ_S(G,"G",H,"H")
#
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")
```
Originを非表示しています。isometric方向?です。
マウスの操作勉強中です。
△AGHを、マウスで赤色にしました。

![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/2b7f278e-5016-096b-de6e-adb9a606b850.png)

# 参考

https://qiita.com/mrrclb48z/items/10af3ce483f9711f3749

https://qiita.com/mrrclb48z/items/dd58b9e1e292de68bfe4

minimumです。

https://math-fun.net/20220119/21562/