"""Probe the imported, partly gyro-tuned temperature scale without hiding leakage."""
from pathlib import Path
import json
import numpy as np
from age_calibration_prototype import fit_reference, reference_period
from inspect_cluster_samples import read_benchmarks

stars = read_benchmarks(Path('inputs/curtis2020-source/Benchmark_Clusters.txt'))
training = [row for row in stars if row['cluster'] == 'Praesepe' and 4200 <= row['teff_k'] <= 6000]
fit = fit_reference(training)
# Curtis2020 Table4; rounded published coefficients. Apply only its delta,
# preserving the supplied central temperature instead of claiming exact recovery.
color_coefficients = [-416.585, 39780., -84190.5, 85203.9, -48225.9, 15598.5, -2694.76, 192.865]
rows = []
for low, high in [(4200, 4600), (5500, 6000)]:
    selected = [row for row in stars if row['cluster'] == 'NGC 6811' and low <= row['teff_k'] < high]
    original_temperature = np.array([row['teff_k'] for row in selected])
    color = np.array([row['color0'] for row in selected]); period = np.array([row['period_days'] for row in selected])
    original_polynomial = np.polynomial.polynomial.polyval(color, color_coefficients)
    for extinction in [.15, .22, .25]:
        changed_color = color - .415 * (extinction - .15)
        temperature = original_temperature + np.polynomial.polynomial.polyval(changed_color, color_coefficients) - original_polynomial
        valid = (temperature >= 4200) & (temperature <= 6000) & (changed_color >= .55) & (changed_color <= 3.25)
        predicted = .67 * (period[valid] / reference_period(fit, temperature[valid])) ** (1 / .62)
        rows.append({'original_temperature_bin_k': [low, high], 'original_stars': len(selected), 'changed_av_mag': extinction,
                     'retained_in_calibration_domain': int(np.sum(valid)), 'outside_domain_not_clamped': int(np.sum(~valid)),
                     'median_temperature_shift_k': float(np.median(temperature - original_temperature)),
                     'median_predicted_age_gyr': float(np.median(predicted))})
report = {'status': 'COMPLETED external-reddening sensitivity; exposes inherited gyro-tuned color scale, not an independent age validation',
          'conditioning': 'Praesepe fit remains fixed. NGC6811 source Av=.15 was refined using gyrochronology in Curtis2019; imported temperatures therefore carry that assumption. Alternatives.22/.25 are source-reviewed independent reddening scenarios, not a posterior. Same original temperature-bin stars are followed; changed temperatures outside support are explicitly excluded, never clamped.',
          'temperature_rule': 'Use change in rounded Curtis2020 Table4 polynomial at color0-.415*(Av-.15), added to the original supplied temperature',
          'cases': rows}
Path('results/reddening-sensitivity-results.json').write_text(json.dumps(report, indent=2) + '\n')
print(json.dumps(report, indent=2))