(途中)空間図形「(その4/5)2025愛知公立高校入試図形 最後の問題3(3)」様を参考にsympyの空間の座標変換の計算で。

・(その5/5)作成中。積分は未です。

(その1/5) sympyとFreeCADで
(その2/5) 三角柱で
(その3/5) ３つの球の交点計算で
(その4/5) 本ページ
(その5/5) 作成中

・作図はありません。
・空間の座標変換の計算後のz座標を四角錐の高さにしています。
　ここまできたら、みんな同じに見えてきました。
　変換後の全座標で体積を計算した方がよかったカモ。

・次は(その1/5)を見て下さい。
　オリジナル
　sympyのweb上での実行方法
　

sympyで

ver0.1

from sympy import *
def get_plane_to_xy_transform(A: Point3D, B: Point3D, C: Point3D):
    # ベクトル AB, AC を計算
    u = Matrix(B - A)
    v = Matrix(C - A)

    # 法線ベクトル（平面の法線）
    n = u.cross(v)
    n_unit = n / n.norm()

    # Z軸ベクトル
    z_axis = Matrix([0, 0, 1])

    # 回転軸（n と Z軸の外積）
    axis = n_unit.cross(z_axis)
    axis_norm = axis.norm()

    # 回転角
    cos_theta = n_unit.dot(z_axis)
    theta = acos(cos_theta)

    # Rodriguesの回転行列
    def rotation_matrix(axis, theta):
        axis = axis / axis.norm()
        x, y, z = axis
        c = cos(theta)
        s = sin(theta)
        C = 1 - c
        return Matrix([
            [x*x*C + c,   x*y*C - z*s, x*z*C + y*s],
            [y*x*C + z*s, y*y*C + c,   y*z*C - x*s],
            [z*x*C - y*s, z*y*C + x*s, z*z*C + c  ]
        ])

    # 回転行列（軸がゼロ＝すでにXY平面の場合は単位行列）
    R = eye(3) if axis_norm == 0 else rotation_matrix(axis, theta)

    # 平行移動ベクトル（点Aを原点へ）
    t = -R * Matrix(A)

    # 同次変換行列（4x4）
    T = R.row_join(t)
    T = T.col_join(Matrix([[0, 0, 0, 1]]))

    return T

# === 使用例 ===
if __name__ == "__main__":
    OA,AB=12,6
    d=AB*Rational(1,2)
    h=sqrt(OA**2-((d*sqrt(2))**2))   
    O=Point(0,0,0)
    A,B,C,D=map(Point,[(-d,-d,-h),(d,-d,-h),(d,d,-h),(-d,d,-h)]) 
    E=O+Rational(2,3)*(A-O)
    F=O+Rational(1,2)*(B-O)
    G=O+Rational(1,2)*(C-O)
    H=O+Rational(2,3)*(D-O)
    print("#",  Rational(1,3)*(get_plane_to_xy_transform(E,F,G)*Matrix([0, 0, 0, 1]))[2] \
              *( Triangle(sss=(E.distance(F),F.distance(H),H.distance(E))).area
                +Triangle(sss=(F.distance(G),G.distance(H),H.distance(F))).area
              ))
# 7*sqrt(14)

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

(テンプレート)

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

Qiita 内