sympyで(オリジナルの手順で)

等脚台形を勉強中

sympyで(ブラーマグプタの公式で)

from sympy import *
k=Symbol("k", positive=True)
O=Point(0,0,0)
A=Point(4*k,-4*k,-4*sqrt(2)*k)
B=Point(7  ,5   ,  -sqrt(2)  )
C=B+Rational(1,4)*A
#
OA=O.distance(A).simplify()
OB=O.distance(B).simplify()
BC=B.distance(C).simplify()
cosAOB=(A-O).dot(B-O)/(OA*OB)
print("#A",cosAOB)
OACB=Rational(1,2)*(OA+BC)*OB*sqrt(1-cosAOB**2)
print("#I",OACB)
AC=A.distance(C).simplify()
print("#U",AC)
# ブラーマグプタの公式
s=Rational(1,2)*(OA+AC+BC+OB)
ans=solve(Eq(OACB,sqrt((s-OA)*(s-AC)*(s-BC)*(s-OB))),k)[0]
print("#E",ans)
# 以下(3)
A=A.subs({k:ans})
C=C.subs({k:ans})
D=Line(O,B).intersection(Line(A,C)) [0]
M=Rational(1,2)*(O+A)
print("#O",M.distance(D))
#A sqrt(19)/19
#I 30*sqrt(2)*k
#U 2*sqrt(9*k**2 - 6*k + 19)
#E 2/3
#O 8*sqrt(2)

sympyとFreeCADのマクロで(作図)

円,平面αの作図の勉強中

import FreeCAD
import Part
import DraftTools
import Draft
import Mesh
############################################################################
# 計算
from sympy import *
k=Symbol("k", positive=True)
O=Point(0,0,0)
A=Point(4*k,-4*k,-4*sqrt(2)*k)
B=Point(7  ,5   ,  -sqrt(2)  )
C=B+Rational(1,4)*A
#
OA=O.distance(A).simplify()
OB=O.distance(B).simplify()
BC=B.distance(C).simplify()
cosAOB=(A-O).dot(B-O)/(OA*OB)
print("#A",cosAOB)
OACB=Rational(1,2)*(OA+BC)*OB*sqrt(1-cosAOB**2)
print("#I",OACB)
AC=A.distance(C).simplify()
print("#U",AC)
# ブラーマグプタの公式
s=Rational(1,2)*(OA+AC+BC+OB)
ans=solve(Eq(OACB,sqrt((s-OA)*(s-AC)*(s-BC)*(s-OB))),k)[0]
print("#E",ans)
# 以下(3)
A=A.subs({k:ans})
C=C.subs({k:ans})
D=Line(O,B).intersection(Line(A,C)) [0]
M=Rational(1,2)*(O+A)
print("#O",M.distance(D))
print(A)
######################################################################
# 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.5 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
myLine_C  (O,A,C,B)
myLine    (D,M)
#myTxtXYZ_S(O,"O",A,"A",B,"B")
myTxtXYZ_S(O,"O",A,"A",B,"B",C,"C",D,"D",M,"M")
#
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")

真上から見た。

参考