Решение мехматской стереометрии

вычислительная_стереометрия.py

from sympy import *
from sympy.abc import x
from sympy.vector import *
CoordSystem = CoordSys3D("C")
oX, oY, oZ = CoordSystem.i, CoordSystem.j, CoordSystem.k

A = oX * 14  # Удачным образом сторона 6-угольника равна радиусу окружости, в которую он вписан
B = (oX / 2 + oY * sqrt(3) / 2) * 14
C = (-oX / 2 + oY * sqrt(3) / 2) * 14
D, E, F = -A, -B, -C

s = symbols('s', positive=True)  # Высота пирамиды
S = s * oZ
K = (3 * C + 4 * B) / 7
# Нормаль к плоскости 𝜋 ортогональна AB (так как плоскость параллельна),
# а также векторному произведению SD на SE (так как плоскость ей ортогональна)
# А значит её нормаль — векторное произведение этих нормале
SDEnormal = (S - D).cross(S - E)
pi_normal = (B - A).cross(SDEnormal)
def intersect_pi_with_segm(A1, A2, pi_pt, pi_normal):
    # Плоскость проходит через pi_pt.
    # Поэтому выбираем такую точку на ребре, чтобы вектор (K->точка) был перпендикулярен нормали
    lin_sum = A1 * x + A2 * (1 - x)
    t = solve((lin_sum - pi_pt).dot(pi_normal), x)[0]
    return lin_sum.subs(x, t).expand().simplify()

L_CS = intersect_pi_with_segm(C, S, K, pi_normal)
M_DS = intersect_pi_with_segm(D, S, K, pi_normal)
N_ES = intersect_pi_with_segm(E, S, K, pi_normal)
O_FS = intersect_pi_with_segm(F, S, K, pi_normal)
P_AF = intersect_pi_with_segm(A, F, K, pi_normal)

# По условию отрезки K-CS и AF-FS параллельны, то есть векторное произведение равно нулю
must_be_zero1 = (K - L_CS).cross(P_AF - O_FS).expand().simplify()

# Дополнительно вспомним, что пересечение плоскости 𝜋 с BCS и AFS — параллельные прямые.
# Если 2 параллельные прямые лежат в пересекающихся плоскостях, то они параллельны прямой их пересечения
# Которую можно получить как векторное произведение нормалей к плоскостям
AFS_BCS_vec = ((B - S).cross(C - S)).cross((A - S).cross(F - S))
must_be_zero2 = (K - L_CS).cross(AFS_BCS_vec).expand().simplify()

# Достаточно приравнять Z-компоненту к нулю
s_val1 = solve(must_be_zero1.dot(oZ), s)[0]
s_val2 = solve(must_be_zero2.dot(oZ), s)[0]
assert s_val1 == s_val2

# Пока всё работает. Подставляем найденное S в координаты точек пересечения плоскости с рёбрами
S = S.subs(s, s_val1)
L_CS = L_CS.subs(s, s_val1)
M_DS = M_DS.subs(s, s_val1)
N_ES = N_ES.subs(s, s_val1)
O_FS = O_FS.subs(s, s_val1)
P_AF = P_AF.subs(s, s_val1)


print('Высота пирамиды', s_val1)
print('Плоскость 𝜋 делит CS в отношении', (L_CS - C).magnitude() / (L_CS - S).magnitude())
print('Плоскость 𝜋 делит DS в отношении', (M_DS - D).magnitude() / (M_DS - S).magnitude())
print('Длина боковых ребёр пирамиды равна', (S-C).magnitude())

# Проверим, что PO и KL ортогональны PK
assert (P_AF-O_FS).dot(P_AF-K) == (K-L_CS).dot(P_AF-K)

Q_KP = K/2 + P_AF/2
R_ED = D/2 + E/2
H_MN = M_DS/2 + N_ES/2;  assert (Q_KP - H_MN).dot(R_ED - S) == 0  # Действительно, QH — высота QRS
print('Длина отрезка SR (в решении 3/2a=21):', (S-R_ED).magnitude())
print('sin HQR (в решении √(1/3)) =', (R_ED - H_MN).magnitude() / (Q_KP - R_ED).magnitude())
print('cos HQR (в решении √(2/3)) =', (Q_KP - H_MN).magnitude() / (Q_KP - R_ED).magnitude())
# Синус угла между гранью и основанием равен векторному произведению нормированных нормалей
ABC_normal = ((B-A).cross(B-C)).normalize()
BCS_normal = ((B-S).cross(B-C)).normalize()
sin_e = ABC_normal.cross(BCS_normal).magnitude(); cos_e = sqrt(1-sin_e**2).expand().simplify()
print('sin угла междую гранью и основанием равен =', sin_e)
print('cos угла междую гранью и основанием равен =', cos_e)

S_p = 0*oZ
L_p = L_CS.dot(oX)*oX + L_CS.dot(oY)*oY
M_p = M_DS.dot(oX)*oX + M_DS.dot(oY)*oY
print("Длина CL' равна (в решении 8/2=4) =", (L_p - C).magnitude())
print("Длина S'M' равна (в решении (2*14-8)/3=20/3) =", (S_p - M_p).magnitude())
print("Длина S'L' равна (в решении (14-8/2)-10) =", (S_p - L_p).magnitude())
print("Площадь треугольника S'M'L' (в решении 50/√3) равна =", (S_p-M_p).cross(S_p-L_p).magnitude() / 2)

# Теперь нужно найти площадь треугольника CS-DS-S (LMS). Это — половина модуля векторного произведения
area = (S - L_CS).cross(S - M_DS)
area = area.magnitude() / 2
print('Искомая площадь (LMS в решении) равна =', area)