LoginSignup
0
0
@mrrclb48z

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

Last updated at Posted at 2024-02-23

・３次元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]

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だからいいような気もしますが？

# 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

# 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をコピー貼り付けです。

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

・拡大図
・Origin非表示はCAD操作です。
・着色はCAD操作です。
2.png

いつもの? sympyの実行環境と 参考のおすすめです。

(テンプレート)

いつもと違うおすすめです。

参考過去問

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0