アポロニウスの球？「2023岡山大学前期数学I・数学II・数学A・数学B第3問」をsympyとFreeCADでやってみたい。

sympyは、2次式にならなくていいね。と思いながら、コーディングしていると、
だんだん雲行きが怪しくなってきました。アドバイスをお願いします。
①b=0? a^2=...?
②作図もしてみたい？

**修正してみた。2023/05/10**



オリジナル
T氏様

https://mathmathmass.exblog.jp/33292525/

kamelink様

https://kamelink.com/blog/2023/04/10/6-6-23-okayamadai-bunkei-3/

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

https://mathexamtest.web.fc2.com/2023/202310701/2023107010100mj.html#top-0107

# sympyで(オリジナルのやり方)
勉強中。

# sympyで(sympy的?安易なやり方)
怪しいです。
```python
from sympy import *
r    =symbols('r'    ,real=True,positive=True)
a,b,c=symbols('a b c',real=True)
cc   =symbols('cc'   ,real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a, b,c)
#
r=1
ans=solve([
          Eq( A.distance(P),r*O.distance(P) ),
          Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])[0]
OP=Point(ans[0],ans[1],c).distance(O)
c2=solve(Eq(OP, minimum(OP,c,Interval(-oo,oo)) ))[0]
print("#(1)",ans[0],ans[1],c2)
#########################################################
r=sqrt(3)/2
ans=solve([
           Eq( A.distance(P),r*O.distance(P) ),
           Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])
amin=min(
         minimum(ans[0][0],c,Interval(-oo,oo)),
         minimum(ans[1][0],c,Interval(-oo,oo)) 
        )
amax=max(
         maximum(ans[0][0],c,Interval(-oo,oo)),
         maximum(ans[1][0],c,Interval(-oo,oo)) 
        )
print("#(2)",amin,",",amax)          
#########################################################
print("#   b =",ans[0][1],ans[1][1])          
c01=solve(Eq( A.distance(P),r*O.distance(P) ),c)          
print("#   a2=",c01[0].subs({b:0})**2)
print("#   a2=",c01[1].subs({b:0})**2)
mySub={b:ans[0][1],c:c01[0].subs({b:0})}
#########################################################
OPAP=P.dot(P-A).subs(mySub).expand()
print("#(3)",          
         maximum(OPAP,a,Interval(amin,amax)),",",
         minimum(OPAP,a,Interval(amin,amax))
     ) 
#(1) 1 0 0
#(2) 4 - 2*sqrt(2) , 2*sqrt(2) + 4
#   b = 0 0
#   a2= -a**2 + 8*a - 8
#   a2= -a**2 + 8*a - 8
#(3) 14*sqrt(2) + 20 , 20 - 14*sqrt(2)
```

①私の環境は,pycharmです。
②よく聞くのは、Jupyterです。
③web上で、上記のソースを「SymPy Live shell」に、コピー貼り付けでもできました。
タブレットで、コピー貼り付けができない障害が発生しています。
黒背景の右上に、マウスを移動すると、コピーマークが発生します。

https://live.sympy.org/


# FreeCADのマクロで作図(r=1)
```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
r    =symbols('r'    ,real=True,positive=True)
a,b,c=symbols('a b c',real=True)
cc   =symbols('cc'   ,real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a, b,c)
#
r=1
ans=solve([
          Eq( A.distance(P),r*O.distance(P) ),
          Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])[0]
OP=Point(ans[0],ans[1],c).distance(O)
c2=solve(Eq(OP, minimum(OP,c,Interval(-oo,oo)) ))[0]
print("#(1)",ans[0],ans[1],c2)
#########################################################
PP=Point( ans[0],ans[1],c2 )
##########################################################################################################
### 3D作図
def my_Circle(x,y,z,a):
    zaxis = App.Vector(0, 0, 1)
    p5 = App.Vector(float(x),float(y),float(z))
    place5 = App.Placement(p5, App.Rotation(zaxis, 0))
    circle5 = Draft.make_circle(float(a), placement=place5)
def my_Sphere(myx,myy,myz,myr):
    sphere = doc.addObject("Part::Sphere", "mySphere2")
    sphere.Radius = float(myr)
    sphere.Angle1 = -90
    sphere.Angle2 =  90
    sphere.Angle3 = 360
    sphere.Placement = App.Placement(App.Vector(float(myx),float(myy),float(myz)), App.Rotation(0, 90, 0))
    FreeCADGui.getDocument('Unnamed').getObject('mySphere2').Transparency = 90
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.1 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)
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return 0
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_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
#
myTxtXYZ_S(PP,"P",A,"A",B,"B",O,"O")
myLine_H  (PP, A,B,O)
#
doc = App.activeDocument()
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")
```
![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/b9ab48a7-15da-d1ff-6961-f6a9d18ff275.png)

# ????? FreeCADのマクロで作図(r=sqrt(3)/2)
何がなんだか、
わかりませんでした。
```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
r    =symbols('r'    ,real=True,positive=True)
a,b,c=symbols('a b c',real=True)
cc   =symbols('cc'   ,real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a, b,c)
#
r=1
ans=solve([
          Eq( A.distance(P),r*O.distance(P) ),
          Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])[0]
OP=Point(ans[0],ans[1],c).distance(O)
c2=solve(Eq(OP, minimum(OP,c,Interval(-oo,oo)) ))[0]
print("#(1)",ans[0],ans[1],c2)
#########################################################
r=sqrt(3)/2
ans=solve([
           Eq( A.distance(P),r*O.distance(P) ),
           Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])
amin=min(
         minimum(ans[0][0],c,Interval(-oo,oo)),
         minimum(ans[1][0],c,Interval(-oo,oo)) 
        )
amax=max(
         maximum(ans[0][0],c,Interval(-oo,oo)),
         maximum(ans[1][0],c,Interval(-oo,oo)) 
        )
print("#(2)",amin,",",amax)          
#########################################################
print("#   b =",ans[0][1],ans[1][1])          
c01=solve(Eq( A.distance(P),r*O.distance(P) ),c)          
print("#   a2=",c01[0].subs({b:0})**2)
print("#   a2=",c01[1].subs({b:0})**2)
#########################################################
#?????????????????????????????????????????????????
myEq=Eq((a-1)**2+c**2+1,(a**2+c**2)*Rational(3,4))
myEq1=myEq.subs({a:amin})
ans=solve(myEq1,c)
Pmin=Point(amin,0,ans[0])
myEq2=myEq.subs({a:(amin+amax)/2})
ans=solve(myEq2,c)
Pminmax=Point((amin+amax)/2,0,ans[0])
myEq3=myEq.subs({a:amax})
ans=solve(myEq3,c)
Pmax=Point(amax,0,ans[0])
#print(Pmin)
#print(Pminmax)
#print(Pmax)
##########################################################################################################
### 3D作図
def my_Circle(x,y,z,a):
    zaxis = App.Vector(0, 0, 1)
    p5 = App.Vector(float(x),float(y),float(z))
    place5 = App.Placement(p5, App.Rotation(zaxis, 0))
    circle5 = Draft.make_circle(float(a), placement=place5)
def my_Sphere(myx,myy,myz,myr):
    sphere = doc.addObject("Part::Sphere", "mySphere2")
    sphere.Radius = float(myr)
    sphere.Angle1 = -90
    sphere.Angle2 =  90
    sphere.Angle3 = 360
    sphere.Placement = App.Placement(App.Vector(float(myx),float(myy),float(myz)), App.Rotation(0, 90, 0))
    FreeCADGui.getDocument('Unnamed').getObject('mySphere2').Transparency = 90
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.1 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)
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return 0
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_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
#
myTxtXYZ_S(A,"A",B,"B",O,"O")
myLine_S  (A,B)
#
myTxtXYZ_S(Pmin,"Pmin",Pminmax,"Pminmax",Pmax,"Pmax")
myLine_H  (Pmin,    A,B,O)
myLine_H  (Pminmax, A,B,O)
myLine_H  (Pmax,    A,B,O)
#
doc = App.activeDocument()
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")
```
![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/fe7b0cf9-56c8-c8e3-8b8f-79c6c51aedf2.png)

# 参考
勉強中。
3Dプロット参考にしました。ありがとうございました。
2023岡山大学前期【３】解答＜岡山進研学院様

https://www.shinken.ac.jp/file/okadaikaitourei/kaitourei/2023/2023okadaibunkeisugakukaitourei.pdf#page=3

https://www.shinken.ac.jp/information/okadai-answer2021.html

# アポロニウスの球？
？勉強

https://www.google.com/search?q=%E3%82%A2%E3%83%9D%E3%83%AD%E3%83%8B%E3%82%A6%E3%82%B9%E3%81%AE%E7%90%83&rlz=1CAUHAM_enJP1038JP1038&oq=%E3%82%A2%E3%83%9D%E3%83%AD%E3%83%8B%E3%82%A6%E3%82%B9&aqs=chrome.7.69i57j0i512l9.13473j0j7&sourceid=chrome&ie=UTF-8


# 修正してみた。2023/05/10

オリジナル(youtube)

https://youtu.be/_aHNOjsqCV0

```python
from sympy import *
a,b,c,r=symbols('a b c r',real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a ,b,c)
AP=Point(A).distance(P)
BP=Point(B).distance(P)
OP=Point(O).distance(P)
#
r_inp=1
ans=solve([Eq(AP,BP),Eq(BP,r*OP).subs({r:r_inp})],[a,b,c])
P_ans=P.subs({a:ans[0][0],b:ans[0][1]})
OP_ans=Point(O).distance(P_ans)
P_ans=P_ans.subs({c:solve(Eq(OP_ans, minimum(OP_ans,c,Interval(-oo,oo)) ))[0]})
print("#(1)",P_ans)
#
r_inp=Rational(1,2)*sqrt(3)
ans=solve([Eq(AP,BP),Eq(BP,r*OP).subs({r:r_inp})],[a,b,c])
mySub0={a:ans[0][0],b:ans[0][1]}
mySub1={a:ans[1][0],b:ans[1][1]}
P0=P.subs(mySub0)
P1=P.subs(mySub1)
# plot(P0.x,P1.x)
a_min=min( minimum(P0.x,c,Interval(-oo,oo)) , minimum(P1.x,c,Interval(-oo,oo))  )
a_max=max( maximum(P0.x,c,Interval(-oo,oo)) , maximum(P1.x,c,Interval(-oo,oo))  )
print("#(2)",a_min,",",a_max) 
dot_min=min( minimum((P0-O).dot(P0-A),c,Interval(-oo,oo)) , minimum((P1-O).dot(P1-A),c,Interval(-oo,oo))  ).simplify()
dot_max=max( maximum((P0-O).dot(P0-A),c,Interval(-oo,oo)) , maximum((P1-O).dot(P1-A),c,Interval(-oo,oo))  ).simplify()
print("#(3)",dot_min,",",dot_max) 
#(1) Point3D(1, 0, 0)
#(2) 4 - 2*sqrt(2) , 2*sqrt(2) + 4
#(3) 20 - 14*sqrt(2) , 14*sqrt(2) + 20
```

？？？？plotコマンドは、縦横比？異なります。

![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/5511e66f-97a9-ee63-ae56-501639fc427e.png)