"""Independently integrate the existing AC18 model using a different polar chart."""
from pathlib import Path
import json
import numpy as np
from numpy.polynomial.legendre import leggauss

ROOT = Path('.')
record = json.loads((ROOT / 'surface-reference.json').read_text())
source = record['quiet']; a = .45; b = .2
quiet_analytic = np.pi * (1 - a / 3 - b / 6)
assert abs(source['quietFlux'] / quiet_analytic - 1) < 1e-12
polynomial = np.polynomial.Polynomial([1 - a - b, a + 2 * b, -b])
mu_polynomial = np.polynomial.Polynomial([0, 1])
x_second = ((mu_polynomial - mu_polynomial ** 3) * polynomial).integ()(1) * np.pi / quiet_analytic
quiet_variance = source['projectedRotationMps'] ** 2 * x_second
assert abs(source['lineVelocityVariance'] / quiet_variance - 1) < 1e-12
nodes, weights = leggauss(192); mu = (nodes + 1) / 2
phi = (np.arange(960) + .5) * 2 * np.pi / 960
radius = np.sqrt(1 - mu ** 2)
x = radius[:, None] * np.cos(phi)[None, :]
y = radius[:, None] * np.sin(phi)[None, :]
z = -np.broadcast_to(mu[:, None], x.shape)
quiet_weight = (weights / 2 * mu * polynomial(mu))[:, None] * 2 * np.pi / len(phi)
reports = []
for case in record['cases']:
    settings = case['settings']; spot = case['spot']
    sigma = np.deg2rad(settings['spotWidthDegrees'])
    dot = x * spot['x'] + y * spot['y'] + z * spot['z']
    intensity = 1 - settings['spotContrast'] * np.exp(-(1 - dot) / sigma ** 2)
    weighted = quiet_weight * intensity
    flux = np.sum(weighted)
    speed = 2 * np.pi * 695700000 / (settings['rotationDays'] * 86400) * np.sin(np.deg2rad(settings['spinInclinationDegrees']))
    velocity = speed * x
    centroid = np.sum(weighted * velocity) / flux
    variance = np.sum(weighted * velocity ** 2) / flux - centroid ** 2
    relative = flux / quiet_analytic
    flux_error = abs(relative - case['flux'])
    centroid_error = abs(centroid - case['centroidMps'])
    variance_error = abs(variance - case['lineVelocityVariance']) / max(1., abs(variance))
    assert flux_error < 2e-8 and centroid_error < 1e-4 and variance_error < 1e-7, (case['name'], flux_error, centroid_error, variance_error)
    assert abs(case['lineCentroidMps'] - case['centroidMps']) < 1e-8
    assert abs(sum(value * value for value in spot.values()) - 1) < 1e-14
    reports.append({'name': case['name'], 'absolute_flux_difference': flux_error, 'absolute_centroid_difference_mps': centroid_error, 'relative_second_moment_difference': variance_error})
small = []
for case in record['cases'][-3:]:
    settings = case['settings']; sigma = np.deg2rad(settings['spotWidthDegrees']); visible_mu = -case['spot']['z']
    deficit_limit = settings['spotContrast'] * 2 * np.pi * sigma ** 2 * visible_mu * polynomial(visible_mu) / quiet_analytic
    relative_error = abs((1 - case['flux']) / deficit_limit - 1)
    small.append({'width_degrees': settings['spotWidthDegrees'], 'small_spot_deficit_limit': deficit_limit, 'actual_deficit': 1 - case['flux'], 'relative_small_spot_correction': relative_error})
assert small[0]['relative_small_spot_correction'] < .005
assert small[0]['relative_small_spot_correction'] < small[1]['relative_small_spot_correction'] < small[2]['relative_small_spot_correction']
pole = [case['flux'] for case in record['cases'] if case['name'].startswith('i0-lat60')]
assert np.ptp(pole) < 1e-12
report = {'status': 'PASS independent reuse-model flux, first/second moments, line centroid and single-spot limits', 'cases': len(reports), 'quiet_flux_analytic': quiet_analytic, 'quiet_flux_native': source['quietFlux'], 'quiet_x_second_moment': x_second, 'maximum_absolute_flux_error': max(row['absolute_flux_difference'] for row in reports), 'maximum_centroid_error_mps': max(row['absolute_centroid_difference_mps'] for row in reports), 'maximum_relative_variance_error': max(row['relative_second_moment_difference'] for row in reports), 'pole_on_phase_flux_range': float(np.ptp(pole)), 'small_spot_limit': small, 'scope': 'Current rigid, static, single-region AC18 owner checked read-only; no differential/lifetime implementation or age validation is implied.'}
(ROOT / 'results/shared-surface-verification.json').write_text(json.dumps(report, indent=2) + '\n')
print(json.dumps(report, indent=2))