# -*- coding: utf-8 -*-
"""
Лаба 3. Конструирование тоннельных обделок сводчатого очертания.
Расчёт по методике (шаги 1–13), SymPy.

Все длины в мм, углы выводятся в градусах.
"""

from sympy import Rational, sin, cos, asin, deg, rad, N, nsimplify


# ---------------------------------------------------------------- справочные данные

# Табл. 1: ширина проезжей части Г (м) по длине тоннеля и категории дороги
G_TABLE = {
    ("до 300", "I-II"): 9.0,
    ("до 300", "III"): 8.5,
    ("до 300", "IV"): 8.0,
    ("более 300", "I-II"): 8.5,
    ("более 300", "III"): 8.0,
    ("более 300", "IV"): 7.0,
}

# Табл. 7: замковое сечение hз (см) — диапазоны; строка = категория дороги
# (в скобках у габарита — значение Г для тоннелей длиной более 300 м той же категории)
HZ_TABLE = {
    "I-II": [(90, 85), (80, 75), (75, 70), (65, 60), (60, 55), (50, 45), 40],  # Г-9 (8,5)
    "III": [(85, 80), (75, 70), (70, 60), (60, 55), (55, 50), 40, 40],         # Г-8,5 (8)
    "IV": [(80, 75), (70, 65), (60, 55), (55, 50), (50, 45), 40, 40],          # Г-8 (7)
}

# Табл. 7: множители к hз для пяты свода / стены / обреза фундамента по группам крепости
def mult_group(f):
    if f <= 3:
        return {"hп": (1.4, 1.3), "hст": (1.6, 1.5), "hо": (2.0, 1.8)}
    if f <= 6:
        return {"hп": (1.3, 1.1), "hст": (1.5, 1.3), "hо": (1.8, 1.6)}
    return {"hп": (1, 1), "hст": (1.1, 1.0), "hо": (1.5, 1.5)}


# Варианты для самостоятельной работы: (длина, категория, свод, порода, f)
VARIANTS = {
    1:  ("до 300", "I-II", "пологий", "скальная", 4),
    2:  ("до 300", "III", "подъемистый", "скальная", 4),
    3:  ("до 300", "I-II", "пологий", "скальная", 8),
    4:  ("более 300", "III", "подъемистый", "слабая", 2),
    5:  ("до 300", "I-II", "пологий", "скальная", 6),
    6:  ("более 300", "III", "подъемистый", "слабая", 3),
    7:  ("до 300", "I-II", "пологий", "слабая", 3),
    8:  ("до 300", "I-II", "пологий", "скальная", 7),
    9:  ("до 300", "I-II", "пологий", "скальная", 5),
    10: ("до 300", "III", "подъемистый", "скальная", 6),
    11: ("до 300", "I-II", "пологий", "слабая", 2),
    12: ("до 300", "I-II", "пологий", "слабая", 1),
    13: ("более 300", "III", "подъемистый", "скальная", 4),
    14: ("более 300", "III", "подъемистый", "слабая", 1),
    15: ("более 300", "I-II", "пологий", "скальная", 9),
    16: ("до 300", "IV", "подъемистый", "скальная", 4),
    17: ("до 300", "IV", "подъемистый", "скальная", 7),
    18: ("более 300", "III", "подъемистый", "скальная", 7),
    19: ("более 300", "I-II", "подъемистый", "слабая", 2),
    20: ("до 300", "I-II", "пологий", "слабая", 1),
}

PASSAGE = 750    # служебный проход П, мм
STRIP = 500      # защитная полоса ЗП, мм
I_CROSS = 0.02   # поперечный уклон проезжей части, 20 промилле


# ---------------------------------------------------------------- расчёт

def compute(length="до 300", category="I-II", vault="пологий", rock="скальная", f=4,
            beta=None, hz=None, hp=None, dz=None, db=None, round_deg=True):
    """Шаги 1–13 методики. Неизвестные параметры подбираются из таблиц по умолчанию
    (диапазоны берутся с верхним значением, как в примере методички)."""
    out = {}

    # --- исходные данные по таблицам
    G = Rational(int(G_TABLE[(length, category)] * 1000))          # ширина проезжей части, мм
    if dz is None:
        dz = 150 if rock == "слабая" else 100                      # зазор: слабая 15 см, скальная 10 см
    dz = Rational(dz)
    if db is None:
        db = 200 if (vault == "пологий" and rock == "скальная") else 0   # бермы — только пологий свод на скальной
    db = Rational(db)

    # --- 1) пролёт в свету
    if category in ("I-II", "I или II"):
        l0 = G + 2 * PASSAGE + 2 * dz + 2 * db
    else:
        l0 = G + PASSAGE + STRIP + 2 * dz

    # --- 2) отношение пролёта к стреле подъёма
    if beta is None:
        beta = Rational(3) if vault == "пологий" else Rational(9, 4)   # пологий 2,5–4; подъемистый 2,1–2,4

    # --- 3) стрела подъёма в свету
    f0 = l0 / beta

    # --- 4) радиус внутреннего очертания свода
    r0 = f0 * (4 + beta**2) / 8

    # --- 5) предварительный угол наклона пяты свода
    sin_phi_pr = l0 / (2 * r0)
    phi_pr = N(deg(asin(sin_phi_pr)))
    phi_pr_use = round(float(phi_pr)) if round_deg else float(phi_pr)

    # --- 6) сечения обделки по табл. 7
    hz_row = HZ_TABLE[category]
    hz_range = hz_row[min(f, 7) - 1]
    if hz is None:
        hz = hz_range[0] if isinstance(hz_range, tuple) else hz_range  # верх диапазона, см
    hz = Rational(hz * 10)                                             # -> мм
    grp = mult_group(f)
    if hp is None:
        hp = round(grp["hп"][0] * float(hz))                           # верх диапазона множителя
    hp = Rational(hp)

    # --- 7) предварительные осевые стрела и пролёт
    f_pr = f0 + Rational(1, 2) * hz - Rational(1, 2) * hp * cos(rad(phi_pr_use))
    l_pr = l0 + hp * sin(rad(phi_pr_use))

    # --- 8) уточнённый угол пяты
    sin_phi = 4 * f_pr * l_pr / (4 * f_pr**2 + l_pr**2)
    phi = N(deg(asin(sin_phi)))
    phi_use = round(float(phi)) if round_deg else float(phi)

    # --- 9) окончательные осевые стрела и пролёт
    f_ = f0 + Rational(1, 2) * hz - Rational(1, 2) * hp * cos(rad(phi_use))
    l_ = l0 + hp * sin(rad(phi_use))

    # --- 10) расчётный (осевой) радиус свода
    r = Rational(1, 2) * l_ / sin(rad(phi_use))

    # --- 11) стрела и пролёт наружного очертания
    f1 = f_ + Rational(1, 2) * hz - Rational(1, 2) * hp * cos(rad(phi_use))
    l1 = l_ + hp * sin(rad(phi_use))

    # --- 12) радиус наружного очертания
    r1 = (4 * f1**2 + l1**2) / (8 * f1)

    # --- 13) сточный треугольник по оси проезжей части
    h_t = Rational(1, 2) * G * Rational(int(I_CROSS * 1000), 1000)

    # сечения по табл. 7 (для оформления)
    hст = round(grp["hст"][0] * float(hz))
    hо = round(grp["hо"][0] * float(hz))
    hл = round(0.8 * float(hz))

    def v(x):
        return round(float(N(x)), 1)

    out["Г"] = float(G); out["П"] = PASSAGE; out["δз"] = int(dz); out["δб"] = int(db)
    out.update({
        "1) l₀": v(l0), "2) β": float(beta), "3) f₀": v(f0), "4) r₀": v(r0),
        "5) sinφпр": round(float(sin_phi_pr), 4), "5) φпр": phi_pr_use,
        "6) hз,мм": float(hz), "6) hп,мм": float(hp), "6) hст,мм": hст, "6) hо,мм": hо, "6) hл,мм": hл,
        "7) fпр": v(f_pr), "7) lпр": v(l_pr),
        "8) sinφ": round(float(sin_phi), 4), "8) φ": phi_use,
        "9) f": v(f_), "9) l": v(l_),
        "10) r": v(r),
        "11) f₁": v(f1), "11) l₁": v(l1),
        "12) r₁": v(r1),
        "13) hт": v(h_t),
    })
    return out


def report(title, res):
    print(f"\n=== {title} ===")
    for k, val in res.items():
        print(f"  {k:14} = {val}")


if __name__ == "__main__":
    import sys

    if len(sys.argv) > 1 and sys.argv[1] == "all":
        cols = ["1) l₀", "3) f₀", "4) r₀", "5) φпр", "6) hз,мм", "6) hп,мм",
                "8) φ", "9) f", "9) l", "10) r", "11) f₁", "11) l₁", "12) r₁", "13) hт"]
        hdr = "N  " + " ".join(f"{c.split(') ')[1]:>9}" for c in cols)
        print(hdr)
        for n, (length, cat, vault, rock, f) in VARIANTS.items():
            res = compute(length, cat, vault, rock, f)
            print(f"{n:<3}" + " ".join(f"{res[c]:>9}" for c in cols))
    elif len(sys.argv) > 1 and sys.argv[1].isdigit():
        n = int(sys.argv[1])
        length, cat, vault, rock, f = VARIANTS[n]
        report(f"Вариант {n}: {length}, {cat}, {vault} свод, {rock} порода, f={f}",
               compute(length, cat, vault, rock, f))
    else:
        # Проверка на примере методички (до 300, I–II, пологий, скальная, f=4;
        # в примере методички hп = 1,4·hз = 910 мм, β = 3)
        report("Пример методички", compute(
            length="до 300", category="I-II", vault="пологий", rock="скальная", f=4,
            beta=3, hz=65, hp=910, dz=100, db=200))