"""Offline verification of actual FEMM outputs; never simulates an FEMM solve.

Reads native JSON, final CSV row, saved model circuit prescriptions and mesh
triangles from the real .ans files, then writes the delivery report/receipt.
"""
from pathlib import Path
import csv
import json
import math
import re
import hashlib
import statistics

folder = Path(__file__).resolve().parent


def read(name):
    return json.loads((folder / name).read_text(encoding="utf-8"))


def zfield(case, field):
    return complex(*case[field]["peak_re_im"])


def mesh_stats(filename):
    lines = (folder / filename).read_text().splitlines()
    start = lines.index("[Solution]") + 1
    nodes = int(lines[start])
    xy = [tuple(map(float, s.split()[:2])) for s in lines[start + 1:start + 1 + nodes]]
    at = start + nodes + 1
    elements = int(lines[at])
    areas = []
    by_label = {}
    for line in lines[at + 1:at + 1 + elements]:
        f = line.split()
        a, b, c = [xy[int(k)] for k in f[:3]]
        area = abs((b[0]-a[0])*(c[1]-a[1]) - (b[1]-a[1])*(c[0]-a[0])) / 2
        areas.append(area)
        by_label.setdefault(int(f[3]), []).append(area)
    return {"nodes": nodes, "elements": elements,
            "total_area_in2": sum(areas), "mean_area_in2": statistics.mean(areas),
            "median_area_in2": statistics.median(areas), "max_area_in2": max(areas),
            "regions": {k: {"elements": len(v), "area_in2": sum(v),
                             "mean_area_in2": statistics.mean(v), "max_area_in2": max(v)}
                        for k, v in by_label.items()}}


def validate(case, prefix):
    i1, i2 = zfield(case, "primary_current_A"), zfield(case, "secondary_winding_current_A")
    v1, v2 = zfield(case, "primary_voltage_V"), zfield(case, "load_voltage_V")
    source = complex(*case["parameters"]["source_voltage_peak_re_im_V"])
    loadz = complex(*case["parameters"]["load_impedance_re_im_ohm"])
    residual = (source-v1, -v2-loadz*i2)
    norm = math.sqrt(sum(abs(z)**2 for z in residual))
    input_power = .5 * (source * i1.conjugate()).real
    load_power = .5 * loadz.real * abs(i2)**2
    winding_power = .5 * (v1*i1.conjugate() + v2*i2.conjugate()).real
    loss = case["volume_losses"]["all"]["total_W"]
    rows = list(csv.DictReader((folder / (prefix + "_iterations.csv")).open()))
    final = rows[-1]
    assert final["stage"] == "final"
    for symbol, value in [("i1", i1), ("i2", i2), ("v1", v1), ("v2", v2)]:
        csv_value = complex(float(final[symbol+"_peak_re_"+("A" if symbol[0] == "i" else "V")]),
                            float(final[symbol+"_peak_im_"+("A" if symbol[0] == "i" else "V")]))
        assert abs(csv_value-value) < 1e-12
    model = (folder / (prefix + "_solved.fem")).read_text()
    for circuit, current in [("Первичная обмотка", i1), ("Вторичная обмотка", i2)]:
        pattern = (r'<CircuitName>\s*=\s*"'+circuit+r'"\s*'
                   r'<TotalAmps_re>\s*=\s*([^\s]+)\s*'
                   r'<TotalAmps_im>\s*=\s*([^\s]+)')
        match = re.search(pattern, model)
        assert match is not None
        assert abs(complex(float(match[1]), float(match[2]))-current) < 1e-12
    assert abs(norm-case["final_kvl_norm_peak_V"]) < 1e-12
    assert abs(input_power-case["input_power_W"]) < 1e-10
    assert abs(load_power-case["load_power_W"]) < 1e-10
    assert abs(winding_power-loss-case["winding_real_power_minus_volume_losses_W"]) < 1e-10
    stats = mesh_stats(prefix + "_solved.ans")
    assert stats["nodes"] == case["nodes"] and stats["elements"] == case["elements"]
    return {"case": case["case"], "all_native_checks_passed": case["all_checks_passed"],
            "csv_and_saved_circuit_currents_match_final_solution": True,
            "kvl_peak_norm_recomputed_V": norm,
            "input_power_recomputed_W": input_power, "load_power_recomputed_W": load_power,
            "winding_power_recomputed_W": winding_power, "volume_loss_W": loss,
            "power_balance_error_W": input_power-load_power-loss,
            "algebraic_residual_power_bound_W": .5*norm*math.sqrt(abs(i1)**2+abs(i2)**2),
            "algebraic_input_load_minus_windings_W": input_power-load_power-winding_power,
            "mesh": stats}


def normalize_historical_metadata(prefix):
    path = folder / (prefix + "_result.json")
    case = read(path.name)
    ratio = case["probe_no_load_voltage_ratio"]
    if "peak_re_im" in ratio:
        raw = folder / (prefix + "_result_run03_raw.json")
        if not raw.exists():
            raw.write_bytes(path.read_bytes())
        case["probe_no_load_voltage_ratio"] = {
            "re_im": ratio["peak_re_im"], "magnitude": ratio["peak_magnitude"],
            "phase_deg": ratio["phase_deg"]}
        case["source_raw_file"] = raw.name
        case["metadata_correction_note"] = "Only the dimensionless voltage-ratio labels were corrected: RMS and peak normalization does not apply to Vsecondary/Vprimary. Solved currents, fields, voltages and powers were not changed; raw run03 JSON is preserved."
        if prefix == "medium":
            case["case"] = "medium"
            case["source_model"] = "mytransformer_medium.fem"
            case["solved_model"] = "medium_solved.fem"
            case["historical_case_alias_note"] = "This was run03's refined case; renamed medium after the strong refinement."
        path.write_text(json.dumps(case, indent=2) + "\n", encoding="utf-8")
    return case


base = normalize_historical_metadata("base")
medium = normalize_historical_metadata("medium")
refined = read("refined_result.json")
verified = [validate(case, prefix) for case, prefix in zip([base, medium, refined], ["base", "medium", "refined"])]
changes = {}
for field in ["primary_current_A", "secondary_winding_current_A", "load_voltage_V"]:
    x, y = zfield(base, field), zfield(refined, field)
    changes[field] = abs(y-x) / abs(y)
for field in ["input_power_W", "load_power_W"]:
    changes[field] = abs(refined[field]-base[field]) / abs(refined[field])
changes["volume_losses_W"] = abs(verified[2]["volume_loss_W"]-verified[0]["volume_loss_W"]) / verified[2]["volume_loss_W"]
element_ratio = refined["elements"] / base["elements"]
mean_area_ratio = verified[2]["mesh"]["mean_area_in2"] / verified[0]["mesh"]["mean_area_in2"]
checks = {
    "all_native_checks_passed": all(c["all_checks_passed"] for c in [base, medium, refined]),
    "actual_element_count_grew_at_least_3fold": element_ratio >= 3,
    "actual_mean_linear_size_shrank_below_0p6": math.sqrt(mean_area_ratio) < .6,
    "currents_and_voltage_mesh_change_below_1_percent": max(changes[k] for k in ["primary_current_A", "secondary_winding_current_A", "load_voltage_V"]) < .01,
    "input_and_load_power_mesh_change_below_2_percent": max(changes[k] for k in ["input_power_W", "load_power_W"]) < .02,
    "loss_power_mesh_change_below_5_percent": changes["volume_losses_W"] < .05,
}
combined = {
    "engine": "FEMM headless engine v5.0 06Oct2026, official portable runtime",
    "interpreter": "micropython", "cases": [base, refined],
    "historical_medium_case": medium, "independent_output_verification": verified,
    "mesh_comparison": {"relative_changes_using_refined_denominator": changes,
                        "element_count_ratio": element_ratio, "mean_area_ratio": mean_area_ratio,
                        "equivalent_mean_linear_size_ratio": math.sqrt(mean_area_ratio),
                        "checks": checks,
                        "note": "One strong refinement repeat; numerical sensitivity validation, not a formal asymptotic convergence study."},
    "checks": checks, "all_checks_passed": all(checks.values()),
    "execution_logs": [{"path": name + ".log", "exit_code": int((folder / (name+".exitcode.txt")).read_text())}
                       for name in ["native_run_03", "native_run_04", "mesh_unit_probe"]],
    "original_source_url": "https://www.femm.info/examples/mytransformer/mytransformer.zip",
    "reference_page_url": "https://www.femm.info/doku/doku.php?id=MyTransformer",
}
(folder / "RESULT.json").write_text(json.dumps(combined, indent=2) + "\n", encoding="utf-8")
report = ["Transformer: native FEMM 5 MicroPython port and real solve verification",
          "Engine: FEMM headless engine v5.0 06Oct2026; internal interpreter micropython.",
          "Execution: native_run_03.log (base + historical medium); native_run_04.log (strong refined only). Both exit 0.",
          "Algorithm: original fixed approximate complex 2x2 impedance correction, 100-iteration limit; residual target tightened from 0.1 to 0.001 V peak. Extra final solve after last current update.",
          "No external femm/COM client; direct native global mi_modifycircprop. Commands become available after open(model).",
          "RMS convention: source peak=120*sqrt(2); all FEMM AC phasors peak. S_average=0.5 V_peak*conj(I_peak). Integrals3/4/6 are already average watts.",
          "Signs: currents enter both windings; Vload=Vsecondary=-Zload*Isecondary; Iload=-Isecondary.",
          "KVL: residual=(Vsource-Vprimary, -Vsecondary-Zload*Isecondary).",
          "Power: Pin=0.5Re(Vsource conj(Iprimary)); Pload=0.5 Re(Zload)|Isecondary|^2; compare Pin-Pload to FEMM volume losses.", ""]
for case, out in zip([base, medium, refined], verified):
    i1, i2 = zfield(case, "primary_current_A"), zfield(case, "secondary_winding_current_A")
    report += ["%s: nodes=%d; elements=%d; iterations=%d" % (out["case"], case["nodes"], case["elements"], case["iteration_count"]),
               "  Iprimary peak=%+.12g%+.12gj A; Isecondary peak=%+.12g%+.12gj A" % (i1.real,i1.imag,i2.real,i2.imag),
               "  Vload=%.12g V RMS; Pin=%.12g W; Pload=%.12g W; Ploss=%.12g W" % (case["load_voltage_V"]["rms_magnitude"], case["input_power_W"],case["load_power_W"],out["volume_loss_W"]),
               "  Final KVL=%.12g V peak; Pin-Pload-Ploss=%.12g W" % (out["kvl_peak_norm_recomputed_V"],out["power_balance_error_W"]),
               "  Mean triangle area=%.12g in^2; CSV, saved .fem currents and actual .ans mesh checked." % out["mesh"]["mean_area_in2"]]
report += ["", "Strong refined/base elements ratio=%.9g; mean area ratio=%.9g; equivalent mean linear size ratio=%.9g" % (element_ratio,mean_area_ratio,math.sqrt(mean_area_ratio)),
           "Relative mesh changes (refined denominator): " + json.dumps(changes),
           "Published legacy peak reference: Iprimary=2.0446-0.1424j A; Isecondary=-5.8061+0.0173j A, yielding |Vload|=41.05550907 V RMS from rounded currents.",
           "Ideal turns ratio=90/260=0.346153846154; ideal secondary=41.5384615385 V RMS. Loaded voltage is lower due to winding/leakage/core effects.",
           "The first medium mesh is preserved. It changed prescribed meshSize from0.2 to0.1 in but actual elements grew only9.31%; strong refinement uses0.025 in and arc2.5 degrees.",
           "No materials, BH curves, coordinates, boundaries, winding signs/turns, source/load, frequency/depth or solver precision changed. Input original copied byte-for-byte; see mesh_derivation.json.",
           "Limits: nonlinear fundamental time-harmonic 2D planar model at60Hz. All Phi_h=0: no nonzero hysteresis loss angle specified; integral3 includes lamination/wire eddy losses. No end-winding/thermal/THD/hardware validation. Two mesh levels are a sensitivity check, not proof of absolute model accuracy.",
           "The canonical base_result.json and medium_result.json have corrected dimensionless-ratio labels. Exact raw native run03 JSON files are preserved as *_result_run03_raw.json; all solved data remain unchanged. The delivered native script writes correct ratio fields.",
           "Overall checks: " + str(all(checks.values())),
           "Official model/lesson: https://www.femm.info/examples/mytransformer/mytransformer.zip ; https://www.femm.info/doku/doku.php?id=MyTransformer"]
(folder / "RESULT.txt").write_text("\n".join(report) + "\n", encoding="utf-8")
files = [p for p in folder.iterdir() if p.is_file() and p.name not in ["files_manifest.json"]]
(folder / "files_manifest.json").write_text(json.dumps({"files": [{"name": p.name, "bytes": p.stat().st_size,
    "sha256": hashlib.sha256(p.read_bytes()).hexdigest()} for p in sorted(files)]}, indent=2) + "\n", encoding="utf-8")
print("Verified real native outputs; overall checks:", all(checks.values()))
print("Actual mesh elements ratio:", element_ratio)
print("Equivalent mean linear size ratio:", math.sqrt(mean_area_ratio))
assert all(checks.values()), checks