# -*- coding: utf-8 -*-
"""Equilibrium identities, and the two load paths against each other.

None of these need a textbook: they are things that must be true of any
correct answer. They are the tests that catch a sign error in a case nobody
wrote a worked example for.
"""
import pytest

from conftest import LENGTH_FT, aspire_pockets, make_panel

from panellift import pci, statics
from panellift.model import Blockout, Lifter, Panel, Rigging, Stage


def _lifters(xs, width=10.0):
    out = []
    for i, x in enumerate(xs):
        for j, y in enumerate((0.25 * width, 0.75 * width)):
            out.append(Lifter(mark='L{0}{1}'.format(i, j), family='T',
                              x=x, y=y, role='strip', capacity_lb=8000.0))
    return out


# ------------------------------------------------------------ equilibrium

@pytest.mark.parametrize('mode', ['rolling', 'static'])
def test_reactions_sum_to_the_factored_weight(solid_panel, stripping,
                                              rolling, static_lift, mode):
    rigging = rolling if mode == 'rolling' else static_lift
    result = statics.analyse(solid_panel, stripping, rigging,
                             reconcile_weight=False)
    total = sum(line.reaction_lb for line in result.rig.lines)
    assert abs(total - result.profile.total_lb) / result.profile.total_lb < 1e-6


@pytest.mark.parametrize('mode', ['rolling', 'static'])
def test_moment_closes_at_both_ends(solid_panel, stripping, rolling,
                                    static_lift, mode):
    rigging = rolling if mode == 'rolling' else static_lift
    result = statics.analyse(solid_panel, stripping, rigging,
                             reconcile_weight=False)
    _, peak = result.diagram.max_positive
    assert abs(result.diagram.moment_at(0.0)) < peak * 1e-3
    assert abs(result.diagram.moment_at(solid_panel.length_ft)) < peak * 1e-3


# ------------------------------------------------------- the closure invariant
#
# `statics.diagram`'s docstring makes a promise about itself: deriving V and M
# from equilibrium rather than from the element forces means "M(L) coming back
# to zero is then a real check on the reactions rather than a tautology". It
# was failing that check. `solve_static` sampled the load once per finite
# element while the diagram integrated the 2,000-station profile, and
# `_nodes` put 1e-4 ft sliver elements at every blockout edge, so the reactions
# missed the load by up to 4.05% and M at the free right end came back at up to
# 70% of peak M. 280 of 549 stage-cases on the 184-file reference sweep closed
# worse than 0.1%; the `rolling` path, which is pure statics, always closed
# exactly. Measured 2026-09-01, fixed the same day.
#
# The tolerance below is deliberately tight. This is an identity, not a
# convergence criterion: the reactions either sum to the load the diagram
# integrates or they do not. 1e-6 of peak M is roughly a hundred times the
# residual the sweep's worst case now shows (5.8e-8) and about four orders
# tighter than anything the old code could have passed.
CLOSURE_TOL = 1e-6


def _closure(result):
    """(sum(R) - integral(w))/integral(w), and |M(L)|/peak|M|."""
    total = result.profile.total_lb
    reactions = sum(line.reaction_lb for line in result.rig.lines)
    peak = max(abs(m) for m in result.diagram.moment_lbft)
    return ((reactions - total) / total if total else 0.0,
            abs(result.diagram.moment_lbft[-1]) / peak if peak else 0.0)


def _closure_cases():
    """One case per (panel shape, stage, load path) worth distinguishing.

    Named rather than generated so a failure says which shape broke: it is the
    blockout cases that used to fail, and a solid panel closes even under the
    old midpoint sampling because w is constant and the one sample is right.
    """
    cases = []
    for panel_name, blockouts in (('solid', None),
                                  ('pocketed', 'aspire'),
                                  ('one big opening', 'big'),
                                  ('opening off a support', 'at support')):
        for stage in (Stage('stripping', 1.4, 4000.0),
                      Stage('transportation', 1.5, 5000.0)):
            for mode in ('rolling', 'static'):
                cases.append(pytest.param(
                    blockouts, stage, mode, None,
                    id='{0}-{1}-{2}'.format(panel_name, stage.name, mode)))
        cases.append(pytest.param(
            blockouts, Stage('transportation', 1.5, 5000.0), 'static',
            [4.0, 27.0],
            id='{0}-on dunnage'.format(panel_name)))
    return cases


def _closure_panel(kind):
    if kind is None:
        return make_panel()
    if kind == 'aspire':
        return make_panel(blockouts=aspire_pockets())
    if kind == 'big':
        # Most of the width gone over an eighth of the span, so w(x) drops by
        # 60% at two stations that no uniform mesh would land on.
        return make_panel(blockouts=[Blockout(name='O1', x0=5.3333, y0=0.0,
                                              length=4.3542, width=6.0)])
    # An edge landing on top of a support line is the case where the mesh has
    # to choose between refining and putting the reaction in the right place.
    left, _ = pci.optimal_two_line_supports(LENGTH_FT)
    return make_panel(blockouts=[Blockout(name='O1', x0=left, y0=0.0,
                                          length=3.0, width=6.0)])


@pytest.mark.parametrize('blockouts,stage,mode,supports', _closure_cases())
def test_every_stage_case_closes_the_moment_at_the_free_end(blockouts, stage,
                                                            mode, supports):
    panel = _closure_panel(blockouts)
    if supports is not None:
        stage = Stage(stage.name, stage.multiplier, stage.fc_psi,
                      supports=supports)
    rigging = Rigging(mode=mode, sling_angle_deg=90.0)
    result = statics.analyse(panel, stage, rigging, reconcile_weight=False)

    reaction_error, end_moment = _closure(result)
    assert abs(reaction_error) < CLOSURE_TOL, (
        'sum(R) misses the load the diagram integrates by {0:.4%}'
        .format(reaction_error))
    assert end_moment < CLOSURE_TOL, (
        'M at the free end is {0:.4%} of peak M; it should be zero'
        .format(end_moment))


def _wp313_transportation():
    """The worst-closing case on the reference sweep, rebuilt from its job file.

    WP313 on dunnage: a 22.5 ft piece with an opening that takes 61% of the
    width out over 3.35 ft, starting 1.33 ft in -- so the load drops by nearly
    two thirds a foot and a half OUTSIDE the left support and comes back three
    feet inside it. Both edges fall between uniform mesh stations, and the
    smaller opening puts a third and fourth edge in mid-span.

    That shape is what turned a one-sample-per-element load into a 4.05%
    reaction error and left M(22.5 ft) at 19,914 lb-ft against a peak of
    28,580 -- 69.7% of peak, on the stage that carries 106 of the 147 capacity
    exceedances the free-end stations produce.
    """
    panel = Panel(
        mark='WP313', length_ft=22.5, width_ft=11.9167, thickness_in=8.0,
        unit_weight_pcf=150.0, fci_psi=3500.0, fc_psi=6000.0,
        finish='smooth', product_type='flat', weight_lb=25711.0,
        lifters=[Lifter(mark='L1', family='T', x=5.1863, y=9.6667,
                        role='strip', capacity_lb=8000.0),
                 Lifter(mark='L2', family='T', x=17.3, y=2.4167,
                        role='strip', capacity_lb=8000.0)],
        blockouts=[Blockout(name='OPENING 1', x0=1.3333, y0=0.0,
                            length=3.3542, width=7.2708),
                   Blockout(name='OPENING 2', x0=8.75, y0=10.3333,
                            length=1.0, width=0.5)])
    stage = Stage('transportation', 1.5, 6000.0, supports=[4.725, 17.775])
    return panel, stage


def test_the_worst_case_on_the_sweep_closes():
    """Regression on the shape that measured 4.05%, with the weight scaling on.

    `reconcile_weight` is left ON here, unlike the parametrised cases above,
    because that is how the sweep runs it: the profile is rescaled to the
    model's own cast weight and the FE loads have to follow the SAME rescaled
    profile. Two representations of the load is the defect; a scaling applied
    to one of them and not the other would be the same defect again.
    """
    panel, stage = _wp313_transportation()
    result = statics.analyse(panel, stage, Rigging(mode='rolling',
                                                   sling_angle_deg=60.0),
                             reconcile_weight=True)
    reaction_error, end_moment = _closure(result)
    assert abs(reaction_error) < CLOSURE_TOL, (
        '{0:.4%}; this case measured 4.0477% before 2026-09-01'
        .format(reaction_error))
    assert end_moment < CLOSURE_TOL, (
        '{0:.4%} of peak M; this case measured 69.67% before 2026-09-01'
        .format(end_moment))


def test_the_element_loads_integrate_the_profile_the_diagram_integrates():
    """The invariant one layer down, so a failure says WHICH half moved.

    The closure tests above can only report that the reactions and the load
    disagree. This one asks the element loads directly: summed over the mesh
    they must reproduce the profile's own total and first moment, because
    those two numbers are the whole of what `diagram` needs to bring M(L) back
    to zero.
    """
    panel, stage = _wp313_transportation()
    profile = statics.load_profile(panel, stage.multiplier)
    rig = statics._support_rig(stage.supports)
    nodes = statics._nodes(panel, rig)

    force = moment = 0.0
    for i in range(len(nodes) - 1):
        x0, x1 = nodes[i], nodes[i + 1]
        f = statics._element_load(profile, x0, x1)
        force -= f[0] + f[2]
        moment -= f[0] * x0 + f[2] * x1 + f[1] + f[3]

    assert abs(force - profile.total_lb) < 1e-6 * profile.total_lb
    assert abs(moment - profile.Q[-1]) < 1e-6 * profile.Q[-1]


def test_the_mesh_never_leaves_a_sliver_element():
    """What made the reactions drift even where the loads were right.

    A node 1e-4 ft from its neighbour beside elements 0.19 ft long is a
    stiffness ratio of 6e9 across two adjacent rows of the matrix, and a dense
    Gaussian solve gives up about that many digits. The mesh is allowed a node
    exactly on a support and exactly on a section jump; it is not allowed one
    a hair away from either.
    """
    panel, stage = _wp313_transportation()
    rig = statics._support_rig(stage.supports)
    nodes = statics._nodes(panel, rig)

    for x in stage.supports:
        assert min(abs(n - x) for n in nodes) < 1e-12, (
            'no node on the support at x={0}'.format(x))
    floor = (panel.length_ft / 120.0) * statics.MIN_ELEMENT_FRACTION
    shortest = min(nodes[i + 1] - nodes[i] for i in range(len(nodes) - 1))
    assert shortest >= floor * (1.0 - 1e-9), (
        'shortest element {0:.3e} ft against a floor of {1:.3e}'
        .format(shortest, floor))


def test_factored_weight_scales_with_the_multiplier(solid_panel, rolling):
    light = statics.analyse(solid_panel, Stage('a', 1.0, 4000.0), rolling,
                            reconcile_weight=False)
    heavy = statics.analyse(solid_panel, Stage('b', 1.3, 4000.0), rolling,
                            reconcile_weight=False)
    _, m_light = light.diagram.max_positive
    _, m_heavy = heavy.diagram.max_positive
    assert abs(m_heavy / m_light - 1.3) < 1e-6


# ------------------------------------------------- rolling versus static

def test_two_supports_agree_between_the_two_load_paths(solid_panel, stripping,
                                                       rolling, static_lift):
    """With two support lines the beam is determinate, so stiffness cannot matter."""
    a = statics.analyse(solid_panel, stripping, rolling, reconcile_weight=False)
    b = statics.analyse(solid_panel, stripping, static_lift,
                        reconcile_weight=False)
    for la, lb in zip(a.rig.lines, b.rig.lines):
        assert abs(la.reaction_lb - lb.reaction_lb) / la.reaction_lb < 1e-3


def test_three_lines_are_indeterminate_and_the_paths_diverge(solid_panel,
                                                             stripping,
                                                             rolling,
                                                             static_lift):
    """A middle support takes far more load than equalised sheaves would give it."""
    solid_panel.lifters = _lifters([4.0, 16.0, 28.0])

    equalised = statics.analyse(solid_panel, stripping, rolling,
                                reconcile_weight=False)
    continuous = statics.analyse(solid_panel, stripping, static_lift,
                                 reconcile_weight=False)

    middle_equalised = equalised.rig.lines[1].reaction_lb
    middle_continuous = continuous.rig.lines[1].reaction_lb
    assert middle_continuous > middle_equalised * 1.2


# --------------------------------------------------------------- grouping

def test_four_lines_split_into_two_cable_groups(solid_panel, rolling):
    solid_panel.lifters = _lifters([3.0, 8.0, 24.0, 29.0])
    rig = statics.group_lifters(solid_panel, rolling)
    assert len(rig.lines) == 4
    assert rig.points == 8
    assert [line.group for line in rig.lines] == [0, 0, 1, 1]
    assert 8.0 < rig.grouped_at < 24.0


def test_three_lines_pair_the_two_that_are_close(solid_panel, rolling):
    solid_panel.lifters = _lifters([4.0, 24.0, 28.0])
    rig = statics.group_lifters(solid_panel, rolling)
    assert [line.group for line in rig.lines] == [0, 1, 1]


def test_lifters_at_the_same_station_collapse_to_one_line(solid_panel, rolling):
    rig = statics.group_lifters(solid_panel, rolling)
    assert len(rig.lines) == 2
    assert rig.points == 4
    assert all(line.count == 2 for line in rig.lines)


def test_one_support_line_is_refused(solid_panel, rolling):
    solid_panel.lifters = _lifters([16.0])
    with pytest.raises(ValueError, match='one line'):
        statics.group_lifters(solid_panel, rolling)


# ------------------------------------------------------------ sling angle

def test_vertical_sling_does_not_amplify(rolling):
    forces = statics.sling_forces(10000.0, rolling)
    assert abs(forces.amplification - 1.0) < 1e-9
    assert abs(forces.horizontal_lb) < 1e-6


def test_sixty_degree_sling_amplifies_by_two_over_root_three(rolling):
    rolling.sling_angle_deg = 60.0
    forces = statics.sling_forces(10000.0, rolling)
    assert abs(forces.line_force_lb - 10000.0 / 0.8660254) < 0.1
    assert abs(forces.horizontal_lb - 5773.5) < 1.0


def test_shallow_sling_is_much_worse(rolling):
    rolling.sling_angle_deg = 30.0
    forces = statics.sling_forces(10000.0, rolling)
    assert abs(forces.line_force_lb - 20000.0) < 1.0


def test_a_zero_sling_angle_refuses_loudly_rather_than_dividing_by_zero(
        rolling):
    rolling.sling_angle_deg = 0.0
    with pytest.raises(ValueError, match='greater than zero'):
        statics.sling_forces(10000.0, rolling)


def test_a_negative_sling_angle_refuses_loudly(rolling):
    """Not physically meaningful, and math.sin of it is negative -- caught
    by the same guard as zero, not a NaN or a silently flipped sign."""
    rolling.sling_angle_deg = -30.0
    with pytest.raises(ValueError, match='greater than zero'):
        statics.sling_forces(10000.0, rolling)


# --------------------------------------------- weight reconciliation

def test_reconciliation_scales_the_load_to_the_stated_weight(solid_panel,
                                                             stripping,
                                                             rolling):
    geometric = statics.load_profile(solid_panel, 1.0).total_lb
    solid_panel.weight_lb = geometric * 1.05

    result = statics.analyse(solid_panel, stripping, rolling,
                             reconcile_weight=True)
    expected = solid_panel.weight_lb * stripping.multiplier
    assert abs(result.profile.total_lb - expected) / expected < 1e-6
    assert abs(result.weight_error + 0.047619) < 1e-4