import math
def integrate_linear_pm_energy(pm_labels, pm_mu, pm_hc, api=None):
    """Integrate unsmoothed constant-element fields; restore smoothing to on.

    pm_labels entries: x, y, direction ('theta', 'theta+180', or numeric deg),
    numeric_angle_deg. Nearest angular label identifies these four separated
    PM sectors. The PM material is uniquely identified by pm_mu in these files.
    SI: Hc in A/m. mo_getprobleminfo returns depth in METRES; mo_getelement
    returns area in squared model units, so volume = area * scale**2 * depth_m.
    This differs from mi_getprobleminfo, whose depth is in model units.
    """
    if api is None:
        api = globals()
    problem = api['mo_getprobleminfo']()
    if problem[0] != 0 or problem[1] != 0:
        raise ValueError('This helper requires a planar magnetostatic model')
    depth_m, scale = float(problem[2]), float(problem[3])
    mu0 = 4 * math.pi * 1e-7
    result = dict(raw_density_integral_J=0.0, physical_energy_J=0.0,
                  coenergy_integral_J=0.0, correction_integral_J=0.0,
                  pm_volume_m3=0.0, pm_element_count=0,
                  Hparallel_min_A_m=1e300, Hparallel_max_A_m=-1e300)
    api['mo_clearblock']()
    api['mo_groupselectblock']()
    raw_energy = complex(api['mo_blockintegral'](2)).real
    raw_coenergy = complex(api['mo_blockintegral'](17)).real
    api['mo_clearblock']()
    api['mo_smooth']('off')
    try:
        for k in range(1, int(api['mo_numelements']()) + 1):
            element = api['mo_getelement'](k)
            x, y, area = map(float, element[3:6])
            values = api['mo_getpointvalues'](x, y)
            bx, by = (complex(values[j]).real for j in (1, 2))
            mux, muy = (complex(values[j]).real for j in (9, 10))
            if mux <= 0 or abs(mux - muy) > 1e-8:
                raise ValueError('Only linear isotropic positive mu is supported')
            volume = area * depth_m * scale ** 2
            wc = (bx * bx + by * by) / (2 * mu0 * mux)
            w, correction = wc, 0.0
            if abs(mux - pm_mu) < 1e-8:
                theta = math.atan2(y, x) * 180 / math.pi
                label = min(pm_labels, key=lambda lab: abs(
                    (theta - math.atan2(lab['y'], lab['x']) * 180 / math.pi + 180) % 360 - 180))
                direction = label['direction']
                if direction == 'theta':
                    alpha = theta
                elif direction == 'theta+180':
                    alpha = theta + 180
                else:
                    alpha = float(direction)
                mx, my = math.cos(alpha * math.pi / 180), math.sin(alpha * math.pi / 180)
                raw_angle = label['numeric_angle_deg'] * math.pi / 180
                rx, ry = math.cos(raw_angle), math.sin(raw_angle)
                br = mu0 * mux * pm_hc
                w = ((bx - br * mx) ** 2 + (by - br * my) ** 2) / (2 * mu0 * mux)
                correction = pm_hc * (bx * (rx - mx) + by * (ry - my))
                hparallel = (bx * mx + by * my - br) / (mu0 * mux)
                result['Hparallel_min_A_m'] = min(result['Hparallel_min_A_m'], hparallel)
                result['Hparallel_max_A_m'] = max(result['Hparallel_max_A_m'], hparallel)
                result['pm_volume_m3'] += volume
                result['pm_element_count'] += 1
            result['raw_density_integral_J'] += complex(values[4]).real * volume
            result['physical_energy_J'] += w * volume
            result['coenergy_integral_J'] += wc * volume
            result['correction_integral_J'] += correction * volume
    finally:
        api['mo_smooth']('on')
    result['raw_energy_J'] = raw_energy
    result['raw_coenergy_J'] = raw_coenergy
    result['corrected_from_raw_J'] = raw_energy + result['correction_integral_J']
    agreements = dict(
        raw_density_vs_native_J=abs(result['raw_density_integral_J'] - raw_energy),
        coenergy_vs_native_J=abs(result['coenergy_integral_J'] - raw_coenergy),
        direct_vs_correction_J=abs(result['physical_energy_J'] - result['corrected_from_raw_J']))
    tolerance_J = 1e-9 * max(1.0, abs(raw_energy), abs(raw_coenergy))
    if max(agreements.values()) > tolerance_J:
        raise ValueError('Element integrals disagree: check units, materials and PM label identification')
    result['integral_agreement_abs_J'] = agreements
    result['integral_tolerance_J'] = tolerance_J
    result['integral_agreement_passed'] = True
    result['volume_formula'] = 'element_area * unit_scale**2 * postprocessor_depth_metres'
    result['method'] = 'unsmoothed planar element centroid sum; smoothing restored on'
    return result