#!/usr/bin/env python3
from typing import Optional
"""
VDI 2230:2003 — Motore di calcolo completo R0-R13
Uso: python vdi2230_calc.py <input.json>
Output: JSON con tutti i risultati e il report testuale
"""

import json
import math
import sys
import os

# ─── Caricamento dati ────────────────────────────────────────────────────────

SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
DATA_DIR   = os.path.join(SCRIPT_DIR, "..", "data")

def load_json(filename):
    path = os.path.join(DATA_DIR, filename)
    with open(path, "r", encoding="utf-8") as f:
        return json.load(f)

THREAD = load_json("thread_data.json")
FM_TAB = load_json("fm_table.json")

# ─── Costanti ────────────────────────────────────────────────────────────────
ES_BOLT = 206000.0   # N/mm² modulo elastico acciaio vite
PI      = math.pi
ALPHA   = PI / 6     # 30° semi-angolo filettatura metrica

# ─── Funzioni ausiliarie ─────────────────────────────────────────────────────

def get_thread_data(size: str, thread_type: str = "coarse") -> dict:
    """Recupera dati filettatura. size es: 'M12', 'M12x1.5'"""
    if thread_type == "coarse":
        data = THREAD["thread_metric_coarse"]
        key = size.upper()
    else:
        data = THREAD["thread_metric_fine"]
        key = size.upper()
    if key not in data:
        raise ValueError(f"Filettatura '{size}' non trovata nei dati ({thread_type}). "
                         f"Disponibili: {list(data.keys())}")
    return data[key]

def get_strength(fk: str) -> dict:
    """Proprietà meccaniche per classe di resistenza '8.8', '10.9', '12.9'"""
    sc = THREAD["strength_classes"]
    if fk not in sc:
        raise ValueError(f"Classe '{fk}' non trovata. Disponibili: {list(sc.keys())}")
    return sc[fk]

def get_fm_table(size: str, strength_class: str, mu: float) -> Optional[dict]:
    """
    Legge FM_Tab e MA da tabella A1 per la µ più vicina disponibile.
    Restituisce {'FM': kN, 'MA': Nm} o None se fuori tabella.
    """
    tab = FM_TAB.get("A1_regular_shank", {})
    size_upper = size.upper()
    if size_upper not in tab:
        return None
    mu_keys = sorted([float(k) for k in tab[size_upper][strength_class].keys()])
    # Trova µ più vicina
    closest = min(mu_keys, key=lambda x: abs(x - mu))
    entry = tab[size_upper][strength_class][f"{closest:.2f}"]
    return {"FM_kN": entry["FM"], "MA_Nm": entry["MA"], "mu_used": closest}

# ─── R0: stima diametro e verifica limite G ──────────────────────────────────

def r0_check_geometry(d: float, lK: float, dW: float, dh: float,
                      DA: float, hmin: float,
                      joint_type: str = "DSV") -> dict:
    """
    joint_type: 'DSV' (vite passante) o 'ESV' (vite in foro filettato)
    Restituisce dizionario con G, cT_max, e flag di validità.
    """
    if joint_type == "DSV":
        G = hmin + dW
        G_type = "DSV"
    else:
        G = 1.75 * dW   # valore medio dell'intervallo (1.5...2)·dW
        G_type = "ESV"

    valid = DA <= G
    return {
        "G": G,
        "G_type": G_type,
        "DA": DA,
        "valid": valid,
        "warning": f"DA={DA:.1f} > G={G:.1f}: calcolo eccentrico potrebbe avere errori significativi!" if not valid else ""
    }

# ─── R1: fattore di serraggio ────────────────────────────────────────────────

def r1_tightening_factor(method: str, mu_class: str = "B") -> float:
    """
    method: 'torque_wrench_classB', 'torque_manual_classC',
            'yield_or_angle_controlled', 'torque_wrench_classA',
            'power_tool_controlled'
    Restituisce αA tipico (valore medio dell'intervallo).
    """
    tf = THREAD["tightening_factors_alphaA"]
    if method not in tf:
        raise ValueError(f"Metodo '{method}' non trovato. Disponibili: {list(tf.keys())}")
    entry = tf[method]
    return (entry["alphaA_min"] + entry["alphaA_max"]) / 2.0

# ─── R2: forza di serraggio minima FKerf ────────────────────────────────────

def r2_min_clamp_force(FQ_max: float = 0.0, MY_max: float = 0.0,
                       mu_T: float = 0.12, qF: int = 1, qM: int = 1,
                       ra: float = 0.0,
                       pi_max: float = 0.0, AD: float = 0.0,
                       FA_max: float = 0.0, ssym: float = 0.0,
                       a: float = 0.0, u: float = 0.0,
                       IBT: float = 0.0, MB_max: float = 0.0) -> dict:
    """
    Calcola FKerf per tre requisiti:
    - FKQ: anti-scorrimento
    - FKP: tenuta
    - FKA: anti-apertura (klaffen)
    Tutti i valori in N (non kN) e N·mm.
    """
    # a) Anti-scorrimento
    FKQ = 0.0
    if FQ_max > 0:
        FKQ += FQ_max / (qF * mu_T)
    if MY_max > 0 and ra > 0:
        FKQ += MY_max / (qM * ra * mu_T)

    # b) Tenuta
    FKP = AD * pi_max if (AD > 0 and pi_max > 0) else 0.0

    # c) Anti-apertura
    FKA = 0.0
    if IBT > 0 and u > 0 and FA_max > 0:
        denom = IBT + ssym * u * AD
        FKA  = FA_max * AD * (a * u - ssym * u) / denom
        FKA += MB_max * u * AD / denom if MB_max > 0 else 0.0
        FKA  = max(0.0, FKA)

    FKerf = max(FKQ, FKP + FKA)

    return {
        "FKQ_N": FKQ,
        "FKP_N": FKP,
        "FKA_N": FKA,
        "FKerf_N": FKerf,
        "governing": "FKQ" if FKQ >= FKP + FKA else "FKP+FKA"
    }

# ─── R3: cedevolezze e fattore di carico Φ ──────────────────────────────────

def r3_bolt_resilience(thread_data: dict, lK: float,
                       lGew: float = 0.0, lGew_sections: list = None,
                       joint_type: str = "DSV",
                       E_internal: float = None) -> float:
    """
    Calcola δS [mm/N] per vite con gambo pieno standard.
    lGew_sections: lista di (li, di) per elementi non-standard nel gambo.
    E_internal: modulo elastico del materiale madrevite (ESV). Default = ES_BOLT.
    """
    d  = thread_data["d"]
    d3 = thread_data["d3"]
    P  = thread_data["P"]
    d2 = thread_data["d2"]
    AN = PI / 4 * d**2
    Ad3= PI / 4 * d3**2

    if E_internal is None:
        E_internal = ES_BOLT

    # Testa (esagonale standard: lSK = 0.5d)
    lSK = 0.5 * d
    dSK = lSK / (ES_BOLT * AN)

    # Gambo (elementi standard — tutto il gambo come sezione piena se non specificato)
    d_shaft = 0.0
    if lGew_sections:
        d_shaft = sum(li / (ES_BOLT * (PI / 4 * di**2)) for li, di in lGew_sections)
    else:
        # Se non ci sono sezioni diverse, usiamo tutto lK come gambo pieno d
        # (approssimazione; nella realtà parte è gambo, parte filetto)
        pass

    # Filetto libero non avvitato
    d_Gew = lGew / (ES_BOLT * Ad3) if lGew > 0 else 0.0

    # Filetto avvitato + zona madrevite
    lG = 0.5 * d
    dG = lG / (ES_BOLT * Ad3)
    lM = 0.4 * d if joint_type == "DSV" else 0.33 * d
    dM = lM / (E_internal * AN)
    dGM = dG + dM

    # Approssimazione: per lK senza dettagli sezioni → tutto il gambo libero in sede piena
    if not lGew_sections:
        l_shaft = lK - lGew   # porzione gambo pieno (cilindrico)
        d_shaft = l_shaft / (ES_BOLT * AN)

    dS = dSK + d_shaft + d_Gew + dGM
    return dS

def r3_plate_resilience(dW: float, dh: float, lK: float, DA: float,
                        EP: float, DA_prime: float = None,
                        joint_type: str = "DSV") -> dict:
    """
    Calcola δP [mm/N] con il modello a cono VDI 2230 §5.1.2.
    Restituisce anche il cono angle φ e il tipo di modello usato.
    """
    w = 1 if joint_type == "DSV" else 2
    if DA_prime is None:
        DA_prime = DA  # approssimazione conservativa

    # Angolo cono
    bL  = lK / dW
    psi = DA_prime / dW

    if joint_type == "DSV":
        tan_phi = 0.362 + 0.032 * math.log(bL / 2 + 1e-9) + 0.153 * math.log(psi)
    else:
        tan_phi = 0.348 + 0.013 * math.log(bL + 1e-9) + 0.193 * math.log(psi)

    # Clamp al range fisico [0.3, 1.0]
    tan_phi = max(0.3, min(tan_phi, 1.0))
    phi_deg = math.degrees(math.atan(tan_phi))

    # Diametro limite
    DA_Gr = dW + w * lK * tan_phi

    if DA >= DA_Gr:
        # Solo coni (no manicotto)
        model = "coni_puri"
        num = 2 * math.log(
            ((dW + dh) * (dW + w * lK * tan_phi - dh)) /
            ((dW - dh) * (dW + w * lK * tan_phi + dh) + 1e-15)
        )
        den = w * EP * PI * dh * tan_phi
        dP  = num / den if den != 0 else 1e-12

    elif dW < DA < DA_Gr:
        # Coni + manicotto
        model = "coni_e_manicotto"
        lV = (DA - dW) / (2 * tan_phi)
        lH = lK - 2 * lV / w

        num_cone = math.log(
            ((dW + dh) * (dW + 2 * lV * tan_phi - dh)) /
            ((dW - dh) * (dW + 2 * lV * tan_phi + dh) + 1e-15)
        )
        dP_cone = num_cone / (EP * PI * dh * tan_phi)
        dP_sleeve = 4 * lH / (EP * PI * (DA**2 - dh**2))
        dP = (2 / w) * dP_cone + dP_sleeve

    else:
        # Solo manicotto (dW >= DA)
        model = "solo_manicotto"
        dP = 4 * lK / (EP * PI * (DA**2 - dh**2))

    return {
        "dP_mm_per_N": abs(dP),
        "tan_phi": tan_phi,
        "phi_deg": phi_deg,
        "DA_Gr": DA_Gr,
        "model": model
    }

def r3_load_factor(dS: float, dP: float, n: float = 1.0,
                   ssym: float = 0.0, a: float = 0.0,
                   lK: float = 0.0, EP: float = 210000.0,
                   IBers: float = None) -> dict:
    """
    Calcola Φ (fattore di carico) e forze addizionali.
    Per caso centrico: ssym=0, a=0.
    Per caso eccentrico: fornire ssym, a, lK, EP, IBers.
    """
    if ssym == 0 and a == 0:
        # Caso centrico
        Phi_n = n * dP / (dS + dP)
        return {"Phi": Phi_n, "type": "centrico",
                "dP_star": dP, "dP_doublestar": dP}

    # Caso eccentrico
    if IBers is None or IBers == 0:
        raise ValueError("IBers richiesto per caso eccentrico (ssym ≠ 0 o a ≠ 0)")

    bP_lK = lK / (EP * IBers)   # flessibilità a flessione (= lK/(EP·IBers))

    dP_star  = dP + ssym**2 * bP_lK
    dP_dstar = dP_star + ssym * a * bP_lK   # può diventare negativo → caso di apertura anticipata

    Phi_star = n * dP_dstar / (dS + dP_star)

    return {
        "Phi": Phi_star,
        "type": "eccentrico",
        "ssym": ssym,
        "a": a,
        "dP_star": dP_star,
        "dP_doublestar": dP_dstar,
        "bP_term": bP_lK
    }

# ─── R4: variazioni di precarico ────────────────────────────────────────────

def r4_preload_changes(fZ_total: float, dS: float, dP: float,
                       lK: float = 0.0,
                       alpha_S: float = 11.5e-6, dT_S: float = 0.0,
                       alpha_P: float = 11.5e-6, dT_P: float = 0.0,
                       E_S_RT: float = 206000.0, E_P_RT: float = 210000.0,
                       E_S_T:  float = None,      E_P_T:  float = None) -> dict:
    """
    fZ_total: settaggio totale [µm]
    Restituisce FZ [N] e delta_F_Vth [N] (semplificato).
    """
    # Settaggio
    FZ = (fZ_total * 1e-3) / (dS + dP)   # fZ in µm → mm / cedevolezze in mm/N → N

    # Variazione termica semplificata (formula R4/2)
    delta_F_Vth = 0.0
    if lK > 0 and (dT_S != 0 or dT_P != 0):
        E_S_T = E_S_T or E_S_RT
        E_P_T = E_P_T or E_P_RT
        num = lK * (alpha_S * dT_S - alpha_P * dT_P)
        den = dS * (E_S_RT / E_S_T) + dP * (E_P_RT / E_P_T)
        delta_F_Vth = num / den if den != 0 else 0.0

    return {
        "FZ_N":           FZ,
        "fZ_total_um":    fZ_total,
        "deltaF_Vth_N":   delta_F_Vth,
        "note_FM_min":    "Usare deltaF_Vth=0 se negativo e non si garantisce sequenza temp→carico" if delta_F_Vth < 0 else ""
    }

# ─── R5/R6: precarico minimo e massimo ──────────────────────────────────────

def r5_r6_assembly_preload(FKerf_N: float, Phi: float, FA_max_N: float,
                            FZ_N: float, dF_Vth_N: float,
                            alpha_A: float) -> dict:
    """
    Calcola FM_min e FM_max. Tutti i valori in N.
    """
    dF_for_min = max(0.0, dF_Vth_N) if dF_Vth_N < 0 else dF_Vth_N
    # Nota: se dF_Vth < 0 e sequenza non garantita → usa 0
    dF_for_min = dF_Vth_N   # l'utente deve gestire il segno prima di chiamare

    FM_min = FKerf_N + (1.0 - Phi) * FA_max_N + FZ_N + dF_for_min
    FM_max = alpha_A * FM_min

    return {
        "FM_min_N":  FM_min,
        "FM_min_kN": FM_min / 1000,
        "FM_max_N":  FM_max,
        "FM_max_kN": FM_max / 1000,
        "alpha_A":   alpha_A
    }

# ─── R7: verifica montaggio e scelta diametro ────────────────────────────────

def r7_assembly_stress(thread_data: dict, strength: dict,
                        FM_max_N: float, mu_G: float,
                        n_util: float = 0.9,
                        d0: float = None) -> dict:
    """
    Verifica che FMzul >= FM_max.
    Calcola anche σred,M e FMzul esatto.
    """
    d  = thread_data["d"]
    d2 = thread_data["d2"]
    d3 = thread_data["d3"]
    P  = thread_data["P"]
    AS = thread_data["AS"]
    AN = PI / 4 * d**2

    Rp02 = strength["Rp0_2min"]

    if d0 is None:
        d0 = (d2 + d3) / 2   # dS per sezione di sforzo standard

    A0 = PI / 4 * d0**2

    # Termine torsionale (formula 5.5/7)
    torsion_ratio = (3 / 2) * (d2 / d0) * (P / (PI * d2) + 1.155 * mu_G)
    FM_zul = A0 * n_util * Rp02 / math.sqrt(1 + 3 * torsion_ratio**2)

    # Lettura tabella (se disponibile)
    # FM_Tab viene letta in run_full_calculation con la classe corretta

    ok = FM_zul >= FM_max_N

    return {
        "FM_zul_N":   FM_zul,
        "FM_zul_kN":  FM_zul / 1000,
        "FM_max_N":   FM_max_N,
        "FM_max_kN":  FM_max_N / 1000,
        "ratio":      FM_zul / FM_max_N,
        "ok":         ok,
        "status":     "✅ OK" if ok else "❌ FALLISCE",
        "remedies":   [] if ok else [
            "Aumentare il diametro nominale (es. M12→M14)",
            "Passare a una classe di resistenza superiore (8.8→10.9→12.9)",
            "Ridurre µG con lubrificante adeguato (classe A o B)",
            "Usare serraggio angolare o al limite di snervamento (αA=1,0)",
        ],
        "d0_used":    d0,
        "torsion_ratio": torsion_ratio
    }

# ─── R8: verifica esercizio ──────────────────────────────────────────────────

def r8_working_stress(thread_data: dict, strength: dict,
                       FM_zul_N: float, Phi: float, FA_max_N: float,
                       mu_G: float, dF_Vth_N: float = 0.0,
                       kt: float = 0.5, d0: float = None) -> dict:
    """
    Calcola σred,B e sicurezza SF. dF_Vth_N: se > 0 usare 0 (salvo garanzia).
    """
    d  = thread_data["d"]
    d2 = thread_data["d2"]
    d3 = thread_data["d3"]
    P  = thread_data["P"]
    Rp02 = strength["Rp0_2min"]

    if d0 is None:
        d0 = (d2 + d3) / 2
    A0 = PI / 4 * d0**2
    WP = PI / 16 * d0**3

    dF_use = 0.0 if dF_Vth_N > 0 else dF_Vth_N
    FS_max = FM_zul_N + Phi * FA_max_N - dF_use

    # Tensioni
    sz_max = FS_max / A0

    # Momento nel filetto e tensione torsionale
    MG = FM_zul_N * (d2 / 2) * (P / (PI * d2) + 1.155 * mu_G)
    tau_max = MG / WP

    # Tensione equivalente
    sigma_red_B = math.sqrt(sz_max**2 + 3 * (kt * tau_max)**2)

    SF = Rp02 / sigma_red_B
    ok = sigma_red_B < Rp02

    return {
        "FS_max_N":      FS_max,
        "FS_max_kN":     FS_max / 1000,
        "sz_max_MPa":    sz_max,
        "tau_max_MPa":   tau_max,
        "sigma_red_B":   sigma_red_B,
        "Rp02_MPa":      Rp02,
        "SF":            SF,
        "ok":            ok,
        "status":        "✅ OK" if ok else "❌ FALLISCE",
        "remedies":      [] if ok else [
            "Aumentare il diametro (riduce σz per maggiore A0)",
            "Aumentare la classe di resistenza (aumenta Rp0.2)",
            "Verificare se il serraggio senza torsione è possibile (kt=0 → σred,B = σz)",
        ],
        "kt":            kt
    }

# ─── R9: fatica ─────────────────────────────────────────────────────────────

def r9_fatigue(thread_data: dict, strength: dict,
               FM_zul_N: float, Phi: float,
               FA_max_N: float, FA_min_N: float,
               thread_treatment: str = "SV") -> dict:
    """
    thread_treatment: 'SV' (rullato prima della tempra) o 'SG' (rullato dopo tempra).
    """
    d   = thread_data["d"]
    AS  = thread_data["AS"]
    Rp02= strength["Rp0_2min"]

    FSA_o = Phi * FA_max_N
    FSA_u = Phi * FA_min_N

    sigma_a = (FSA_o - FSA_u) / (2 * AS)

    # Tensione media
    FSm = (FSA_o + FSA_u) / 2 + FM_zul_N
    F02_min = Rp02 * AS

    ratio = FSm / F02_min
    ratio = max(0.3, min(ratio, 0.999))   # clamp al range di validità

    # Limite di fatica
    sigma_ASV = 0.85 * (150 / d + 45)
    if thread_treatment == "SG":
        sigma_AS = (2 - ratio) * sigma_ASV
    else:
        sigma_AS = sigma_ASV

    SD = sigma_AS / sigma_a if sigma_a > 0 else float("inf")
    ok = sigma_a <= sigma_AS

    return {
        "FSA_o_N":       FSA_o,
        "FSA_u_N":       FSA_u,
        "sigma_a_MPa":   sigma_a,
        "sigma_AS_MPa":  sigma_AS,
        "sigma_ASV_MPa": sigma_ASV,
        "SD":            SD,
        "FSm_ratio":     ratio,
        "treatment":     thread_treatment,
        "ok":            ok,
        "status":        "✅ OK" if ok else "❌ FALLISCE",
        "remedies":      [] if ok else [
            "Aumentare il diametro (AS maggiore → σa minore)",
            "Usare filetto rullato post-tempra SG (+20-30% su σAS)",
            "Aumentare il precarico FM (riduce l'ampiezza relativa FSA)",
            "Verificare se il numero di cicli consente progetto a vita finita",
        ],
        "SD_recommended": "SD ≥ 1,2 raccomandato VDI"
    }

# ─── R10: pressione superficiale ────────────────────────────────────────────

def r10_surface_pressure(FM_zul_N: float, dW: float, dh: float,
                          FA_max_N: float, Phi: float,
                          pG_N_mm2: float,
                          dF_Vth_N: float = 0.0,
                          FV_max_N: float = None,
                          alpha_A: float = None,
                          FM_tab_N: float = None) -> dict:
    """
    Verifica pressione superficiale sotto testa/dado.
    """
    # Area appoggio (anello circolare con svasatura al foro)
    DKi = dh   # diametro interno appoggio ≈ diametro foro
    # dwa = diametro esterno appoggio testa vite ≈ dW (senza rondella)
    Ap_min = PI / 4 * (dW**2 - DKi**2)

    if Ap_min <= 0:
        return {"error": f"Area appoggio negativa: dW={dW}, dh={dh}"}

    # Pressione montaggio
    pM_max = FM_zul_N / Ap_min
    ok_M   = pM_max <= pG_N_mm2

    # Pressione esercizio
    dF_use = 0.0 if dF_Vth_N > 0 else dF_Vth_N
    if FV_max_N is None:
        FV_max_N = FM_zul_N   # conservativo
    FSA_max = Phi * FA_max_N
    pB_max  = (FV_max_N + FSA_max - dF_use) / Ap_min
    ok_B    = pB_max <= pG_N_mm2

    # Per serraggio angolare/snervamento: pmax = FM_Tab/Ap * 1.4
    p_yield = None
    if FM_tab_N is not None:
        p_yield = FM_tab_N / Ap_min * 1.4

    return {
        "Ap_min_mm2":  Ap_min,
        "pG_MPa":      pG_N_mm2,
        "pM_max_MPa":  pM_max,
        "ok_montaggio":ok_M,
        "pB_max_MPa":  pB_max,
        "ok_esercizio":ok_B,
        "SP_montaggio":pG_N_mm2 / pM_max if pM_max > 0 else 999,
        "SP_esercizio": pG_N_mm2 / pB_max if pB_max > 0 else 999,
        "p_yield_MPa": p_yield,
        "status": "✅ OK" if (ok_M and ok_B) else "❌ FALLISCE",
        "remedies": [] if (ok_M and ok_B) else [
            "Aggiungere rondella (aumenta Ap: dWa = dW + 1,6·hS)",
            "Usare materiale più resistente sotto testa (pG maggiore)",
            "Ridurre FM_zul (solo se R7 rimane verificato)",
        ]
    }

# ─── R11: lunghezza avvitamento ──────────────────────────────────────────────

def r11_engagement_length(d: float, meff_actual: float,
                           strength_class: str,
                           internal_material: str = "steel") -> dict:
    """
    Verifica meff_actual >= meff_min.
    internal_material: 'steel', 'cast_iron', 'aluminum'
    """
    el = THREAD["min_engagement_length_ratio"]
    key = f"{internal_material}_{strength_class}"

    if key not in el:
        # fallback
        ratio = el.get(f"steel_{strength_class}", 1.0)
    else:
        ratio = el[key]

    meff_min = ratio * d
    ok = meff_actual >= meff_min

    return {
        "meff_min_mm":    meff_min,
        "meff_actual_mm": meff_actual,
        "ratio":          meff_actual / meff_min,
        "ok":             ok,
        "status":         "✅ OK" if ok else f"❌ FALLISCE — meff attuale {meff_actual:.1f} mm < minimo {meff_min:.1f} mm",
        "remedies":       [] if ok else [
            f"Aumentare la lunghezza di avvitamento a ≥ {meff_min:.1f} mm (approfondire il foro o usare dado più alto)",
            "Usare inserto filettato (Helicoil/Böllhoff Ensat/Kerb-Konus): aumenta meff disponibile, riduce la pressione sul filetto e ripristina la resistenza dell'acciaio anche in alluminio o ghisa",
            "Passare a classe di resistenza inferiore (riduce meff_min richiesta)",
            "Passare a filettatura fine (stesso diametro esterno, passo ridotto → meff_min leggermente minore)",
        ]
    }

# ─── R12: anti-scorrimento e taglio ─────────────────────────────────────────

def r12_slip_shear(FM_zul_N: float, alpha_A: float, Phi: float,
                   FA_max_N: float, FZ_N: float, dF_Vth_N: float,
                   FKQ_erf_N: float,
                   FQ_max_N: float = 0.0, At_mm2: float = None,
                   tau_B_MPa: float = None, d3: float = None) -> dict:
    """
    Calcola FKR_min e verifica SG e SA.
    """
    dF_use = 0.0 if dF_Vth_N < 0 else dF_Vth_N   # caso conservativo
    FKR_min = FM_zul_N / alpha_A - (1 - Phi) * FA_max_N - FZ_N - dF_use

    SG = FKR_min / FKQ_erf_N if FKQ_erf_N > 0 else float("inf")
    ok_slip = FKR_min > FKQ_erf_N

    # Taglio
    ok_shear = True
    SA = None
    tau_Q_max = None
    if FQ_max_N > 0 and At_mm2 and tau_B_MPa:
        tau_Q_max = FQ_max_N / At_mm2
        SA = tau_B_MPa * At_mm2 / FQ_max_N
        ok_shear = SA >= 1.1

    return {
        "FKR_min_N":   FKR_min,
        "FKR_min_kN":  FKR_min / 1000,
        "FKQ_erf_N":   FKQ_erf_N,
        "SG":          SG,
        "ok_slip":     ok_slip,
        "status_slip": "✅ OK" if ok_slip else "❌ FALLISCE",
        "remedies_slip": [] if ok_slip else [
            "Aumentare il precarico FM (cambio metodo di serraggio o diametro)",
            "Aumentare µT con trattamento superficiale (sabbiatura, ecc.)",
            "Aggiungere spine di posizionamento per trasmettere FQ in forma",
        ],
        "SG_limit":    "SG ≥ 1,2 (statico) | ≥ 1,8 (alternato)",
        "tau_Q_max":   tau_Q_max,
        "SA":          SA,
        "ok_shear":    ok_shear,
        "status_shear":"✅ OK" if ok_shear else "❌ FALLISCE (taglio)"
    }

# ─── R13: coppia di serraggio ────────────────────────────────────────────────

def r13_tightening_torque(thread_data: dict, FM_zul_N: float,
                           mu_G: float, mu_K: float,
                           dW: float, dh: float) -> dict:
    """
    Calcola MA e le componenti di coppia.
    """
    d  = thread_data["d"]
    d2 = thread_data["d2"]
    P  = thread_data["P"]

    DKi  = dh   # diametro interno appoggio
    DKm  = (dW + DKi) / 2   # diametro medio appoggio

    # Momento filetto + appoggio (FM[N]*dist[mm] = N·mm → /1000 = N·m)
    MA        = FM_zul_N * (0.16 * P + 0.58 * d2 * mu_G + DKm / 2 * mu_K) / 1000

    # Componenti
    M_thread  = FM_zul_N * (0.16 * P + 0.58 * d2 * mu_G) / 1000
    M_bearing = FM_zul_N * DKm / 2 * mu_K / 1000

    return {
        "MA_Nm":         MA,
        "M_thread_Nm":   M_thread,
        "M_bearing_Nm":  M_bearing,
        "DKm_mm":        DKm,
        "FM_zul_kN":     FM_zul_N / 1000
    }

# ─── Calcolo completo orchestrato ───────────────────────────────────────────

def run_full_calculation(inp: dict) -> dict:
    """
    Esegue R0–R13 dato un dizionario di input.
    """
    results = {"input": inp, "steps": {}, "summary": {}}
    errors  = []

    try:
        # === Dati base ===
        size    = inp["size"]           # es "M12"
        fk      = inp["strength_class"] # "8.8","10.9","12.9"
        thread_type = inp.get("thread_type", "coarse")

        td = get_thread_data(size, thread_type)
        sc = get_strength(fk)

        lK       = inp["lK_mm"]
        DA       = inp["DA_mm"]
        dh       = inp["dh_mm"]
        dW       = td["dW_hex"]
        joint    = inp.get("joint_type", "DSV")   # DSV o ESV

        DA_prime = inp.get("DA_prime_mm", DA * 1.5)
        hmin     = inp.get("hmin_mm", lK / 2)     # approssimazione conservativa
        ssym     = inp.get("ssym_mm", 0.0)
        a_dist   = inp.get("a_mm", 0.0)
        n_factor = inp.get("n_factor", 1.0)

        FA_max   = inp["FA_max_N"]
        FA_min   = inp.get("FA_min_N", FA_max)   # statico di default
        mu_G     = inp.get("mu_G", 0.12)
        mu_K     = inp.get("mu_K", mu_G)
        alpha_A  = inp.get("alpha_A", 1.7)        # torque_wrench_classB tipico

        EP       = inp.get("EP_N_mm2", 210000.0)
        pG       = inp.get("pG_N_mm2", 900.0)
        lGew     = inp.get("lGew_mm", 0.0)
        meff     = inp.get("meff_actual_mm", td["d"] * 1.5)

        FQ_max   = inp.get("FQ_max_N", 0.0)
        MY_max   = inp.get("MY_max_Nmm", 0.0)
        mu_T     = inp.get("mu_T", 0.12)
        qF       = inp.get("qF", 1)
        qM       = inp.get("qM", 1)
        ra       = inp.get("ra_mm", 0.0)
        pi_max   = inp.get("pi_max_N_mm2", 0.0)
        AD       = inp.get("AD_mm2", 0.0)
        IBT      = inp.get("IBT_mm4", 0.0)
        u_dist   = inp.get("u_mm", 0.0)
        MB_max   = inp.get("MB_max_Nmm", 0.0)

        # Variazioni termiche
        fZ_um    = inp.get("fZ_total_um", 10.0)   # settaggio default conservativo
        alpha_S  = inp.get("alpha_S", 11.5e-6)
        alpha_P  = inp.get("alpha_P", 11.5e-6)
        dT_S     = inp.get("dT_S_K", 0.0)
        dT_P     = inp.get("dT_P_K", 0.0)

        int_mat  = inp.get("internal_material", "steel")
        thread_treat = inp.get("thread_treatment", "SV")
        at_shear = inp.get("At_shear_mm2", td["Ad3"])
        tau_B    = inp.get("tau_B_MPa", 0.6 * sc["Rm_min"])

        # Calcolo IBers se non fornito (approssimazione cilindrica)
        IBers_def = PI / 64 * (DA**4 - dh**4)
        IBers    = inp.get("IBers_mm4", IBers_def)

        # R0
        r = r0_check_geometry(td["d"], lK, dW, dh, DA, hmin, joint)
        results["steps"]["R0"] = r

        # R2
        r = r2_min_clamp_force(FQ_max, MY_max, mu_T, qF, qM, ra,
                                pi_max, AD, FA_max, ssym, a_dist,
                                u_dist, IBT, MB_max)
        results["steps"]["R2"] = r
        FKerf = r["FKerf_N"]

        # R3 - cedevolezze
        dS = r3_bolt_resilience(td, lK, lGew, joint_type=joint)
        rP = r3_plate_resilience(dW, dh, lK, DA, EP, DA_prime, joint)
        dP = rP["dP_mm_per_N"]
        results["steps"]["R3_dS"] = {"dS_mm_per_N": dS}
        results["steps"]["R3_dP"] = rP

        # R3 - fattore di carico
        rPhi = r3_load_factor(dS, dP, n_factor, ssym, a_dist, lK, EP, IBers)
        results["steps"]["R3_Phi"] = rPhi
        Phi = rPhi["Phi"]

        # R4
        rChg = r4_preload_changes(fZ_um, dS, dP, lK, alpha_S, dT_S, alpha_P, dT_P)
        results["steps"]["R4"] = rChg
        FZ   = rChg["FZ_N"]
        dFVth= rChg["deltaF_Vth_N"]

        # Gestione segno dFVth per FM_min
        dFVth_for_min = 0.0 if dFVth < 0 else dFVth

        # R5/R6
        rFM = r5_r6_assembly_preload(FKerf, Phi, FA_max, FZ, dFVth_for_min, alpha_A)
        results["steps"]["R5_R6"] = rFM
        FM_min = rFM["FM_min_N"]
        FM_max = rFM["FM_max_N"]

        # R7
        rR7 = r7_assembly_stress(td, sc, FM_max, mu_G)
        # Leggi FM_Tab
        fm_tab = get_fm_table(size, fk, mu_G)
        if fm_tab:
            FM_tab_N = fm_tab["FM_kN"] * 1000
            MA_tab   = fm_tab["MA_Nm"]
            rR7["FM_Tab_kN"] = fm_tab["FM_kN"]
            rR7["mu_used"]   = fm_tab["mu_used"]
            rR7["ok_tab"]    = FM_tab_N >= FM_max
            rR7["status"]    = ("✅ OK" if FM_tab_N >= FM_max else
                                "❌ FALLISCE — aumentare diametro o classe") + " (tabella)"
        results["steps"]["R7"] = rR7
        FM_zul = rR7["FM_zul_N"]

        # R8
        rR8 = r8_working_stress(td, sc, FM_zul, Phi, FA_max, mu_G, dFVth)
        results["steps"]["R8"] = rR8

        # R9
        rR9 = r9_fatigue(td, sc, FM_zul, Phi, FA_max, FA_min, thread_treat)
        results["steps"]["R9"] = rR9

        # R10
        rR10 = r10_surface_pressure(FM_zul, dW, dh, FA_max, Phi, pG, dFVth,
                                     FM_tab_N=FM_tab_N if fm_tab else None)
        results["steps"]["R10"] = rR10

        # R11
        rR11 = r11_engagement_length(td["d"], meff, fk, int_mat)
        results["steps"]["R11"] = rR11

        # R12
        rR12 = r12_slip_shear(FM_zul, alpha_A, Phi, FA_max, FZ, dFVth,
                               FKerf if FQ_max > 0 else 0.0,
                               FQ_max, at_shear, tau_B, td["d3"])
        results["steps"]["R12"] = rR12

        # R13
        rR13 = r13_tightening_torque(td, FM_zul, mu_G, mu_K, dW, dh)
        results["steps"]["R13"] = rR13

        # Sommario verifiche
        verifications = {
            "R7_montaggio":    rR7["ok"],
            "R8_esercizio":    rR8["ok"],
            "R9_fatica":       rR9["ok"],
            "R10_pressione":   rR10["ok_montaggio"] and rR10["ok_esercizio"],
            "R11_avvitamento": rR11["ok"],
            "R12_scorrimento": rR12["ok_slip"]
        }
        if FQ_max > 0:
            verifications["R12_taglio"] = rR12["ok_shear"]

        all_ok = all(verifications.values())
        results["summary"] = {
            "size":            size,
            "strength_class":  fk,
            "thread_type":     thread_type,
            "joint_type":      joint,
            "FM_min_kN":       FM_min / 1000,
            "FM_max_kN":       FM_max / 1000,
            "FM_zul_kN":       FM_zul / 1000,
            "Phi":             Phi,
            "MA_Nm":           rR13["MA_Nm"],
            "MA_tab_Nm":       MA_tab if fm_tab else None,
            "verifications":   verifications,
            "all_ok":          all_ok,
            "global_status":   "✅ DIMENSIONAMENTO OK" if all_ok else "❌ VERIFICHE FALLITE — vedere dettaglio"
        }

    except Exception as e:
        results["error"] = str(e)
        import traceback
        results["traceback"] = traceback.format_exc()

    return results

# ─── Stampa report leggibile ─────────────────────────────────────────────────

def print_report(res: dict):
    s = res.get("summary", {})
    steps = res.get("steps", {})

    print("=" * 70)
    print("  VDI 2230:2003 — REPORT CALCOLO GIUNZIONE BULLONATA")
    print("=" * 70)

    if "error" in res:
        print(f"\n❌ ERRORE: {res['error']}")
        return

    print(f"\n  Vite: {s.get('size','')} — Classe {s.get('strength_class','')} — {s.get('thread_type','')} — {s.get('joint_type','')}")
    print()
    print("  ── RISULTATI CHIAVE ──────────────────────────────────────")
    print(f"  FM min       = {s.get('FM_min_kN',0):.2f} kN")
    print(f"  FM max       = {s.get('FM_max_kN',0):.2f} kN")
    print(f"  FM zul (calc)= {s.get('FM_zul_kN',0):.2f} kN")
    print(f"  Phi (fattore carico) = {s.get('Phi',0):.4f}")
    print(f"  MA (calc)    = {s.get('MA_Nm',0):.1f} Nm")
    if s.get('MA_tab_Nm'):
        print(f"  MA (tabella) = {s.get('MA_tab_Nm'):.1f} Nm")
    print()
    print("  ── VERIFICHE ─────────────────────────────────────────────")
    for k, v in s.get("verifications", {}).items():
        symbol = "✅" if v else "❌"
        print(f"  {symbol} {k}")
    print()
    print(f"  ══ {s.get('global_status','')}")

    # Dettaglio cedevolezze
    if "R3_dS" in steps and "R3_dP" in steps:
        dS = steps["R3_dS"].get("dS_mm_per_N", 0)
        dP = steps["R3_dP"].get("dP_mm_per_N", 0)
        print(f"\n  Cedevolezze: δS = {dS:.3e} mm/N | δP = {dP:.3e} mm/N")
        print(f"  Modello cono: {steps['R3_dP'].get('model','')} | tan φ = {steps['R3_dP'].get('tan_phi',0):.3f}")

    # Dettaglio R8
    if "R8" in steps:
        r = steps["R8"]
        print(f"\n  R8 Esercizio: σz = {r.get('sz_max_MPa',0):.1f} MPa | τ = {r.get('tau_max_MPa',0):.1f} MPa | σred,B = {r.get('sigma_red_B',0):.1f} MPa | SF = {r.get('SF',0):.2f}")

    # Dettaglio R9
    if "R9" in steps:
        r = steps["R9"]
        print(f"  R9 Fatica: σa = {r.get('sigma_a_MPa',0):.1f} MPa | σAS = {r.get('sigma_AS_MPa',0):.1f} MPa | SD = {r.get('SD',0):.2f}")

    # Dettaglio R10
    if "R10" in steps:
        r = steps["R10"]
        print(f"  R10 Pressione: pM = {r.get('pM_max_MPa',0):.0f} MPa | pB = {r.get('pB_max_MPa',0):.0f} MPa | pG = {r.get('pG_MPa',0):.0f} MPa")

    print("=" * 70)

# ─── Entry point ─────────────────────────────────────────────────────────────

if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("Uso: python vdi2230_calc.py <input.json>  [output.json]")
        sys.exit(1)

    with open(sys.argv[1], "r", encoding="utf-8") as f:
        input_data = json.load(f)

    result = run_full_calculation(input_data)
    print_report(result)

    if len(sys.argv) >= 3:
        with open(sys.argv[2], "w", encoding="utf-8") as f:
            json.dump(result, f, indent=2, ensure_ascii=False)
        print(f"\nRisultati salvati in: {sys.argv[2]}")