#!/usr/bin/env python3
"""Offline Cassini UVIS v3.0 archived binned-count replay, Python3.10+/NumPy.

Spatial wavelength is radial spacing of a ring structure, not UV photon wavelength.
No private repository, Lab code, network request or raw telemetry is required.
"""
import argparse
import hashlib
import json
import math
from pathlib import Path
import numpy as np

HERE = Path(__file__).resolve().parent
ROOT = HERE.parents[2] if len(HERE.parents) >= 3 else HERE
REPOSITORY_INPUT = ROOT/'src/features/labs/planetary-rings/data/cassini-uvis-alpara-2008-261-cutout.tab'
DEFAULT_INPUT = REPOSITORY_INPUT if REPOSITORY_INPUT.exists() else HERE/REPOSITORY_INPUT.name
INPUT_HASH = '6baadf1d91701ad8e868745f7aad9e49d853cec169ce45fea51b9a4f30fdd94c'
OPENING_DEGREES = 54.4


def derive(rows):
    radius, mean, published, star, background, count, flag = [rows[:, column] for column in [0, 3, 4, 8, 9, 10, 11]]
    valid = (flag == 0) & (count > 0) & (mean > background) & (star > 0)
    transmission = np.full_like(mean, np.nan)
    np.divide(mean-background, star, out=transmission, where=valid)
    sine_opening = math.sin(math.radians(OPENING_DEGREES))
    normal_depth = np.full_like(mean, np.nan)
    normal_depth[valid] = -sine_opening*np.log(transmission[valid])
    # Mean counts per raw integration, with n samples averaged in this radial bin.
    mean_variance = np.zeros_like(mean)
    np.divide(mean, count, out=mean_variance, where=valid)
    sigma_depth = np.full_like(mean, np.nan)
    np.divide(sine_opening*np.sqrt(mean_variance), mean-background, out=sigma_depth, where=valid)
    return {'radius_km': radius, 'transmission': transmission, 'normal_depth': normal_depth, 'published_normal_depth': published, 'sigma_normal_depth_counting_only': sigma_depth, 'valid': valid}


def local_search(radius, depth, sigma, lower, upper, seed):
    selected = (radius >= lower) & (radius < upper)
    x = radius[selected]-lower
    values = depth[selected]
    errors = sigma[selected]
    minimum_wavelength = max(3, 2 * float(np.median(np.diff(radius[selected]))))
    frequencies = np.linspace(1/60, 1/minimum_wavelength, 1200)
    baseline = np.column_stack([np.ones(len(x)), x/100, (x/100)**2])/errors[:,None]
    orthogonal = np.linalg.qr(baseline)[0]
    def remove_baseline(vectors):
        return vectors-np.einsum('ij,kj->ki',orthogonal,np.einsum('ki,ij->kj',vectors,orthogonal,optimize=False),optimize=False)
    cosine = remove_baseline(np.cos(2*np.pi*frequencies[:,None]*x)/errors)
    sine = remove_baseline(np.sin(2*np.pi*frequencies[:,None]*x)/errors)
    cosine /= np.linalg.norm(cosine,axis=1)[:,None]
    sine -= np.sum(sine*cosine,axis=1)[:,None]*cosine
    sine /= np.linalg.norm(sine,axis=1)[:,None]
    z = values/errors
    def improvements(observations):
        return np.einsum('ij,jk->ik',cosine,observations,optimize=False)**2+np.einsum('ij,jk->ik',sine,observations,optimize=False)**2
    score = improvements(z[:,None])[:,0]
    best = int(np.argmax(score))
    frequency = frequencies[best]
    matrix = np.column_stack([np.ones(len(x)), x/100, (x/100)**2, np.cos(2*np.pi*frequency*x), np.sin(2*np.pi*frequency*x)])/errors[:,None]
    fit = np.linalg.lstsq(matrix,z,rcond=None)[0]
    residual = z-np.einsum('ij,j->i',matrix,fit,optimize=False)
    random = np.random.default_rng(seed)
    perturbed = z[:,None]+random.normal(size=(len(x),256))
    wavelengths = 1/frequencies[np.argmax(improvements(perturbed),axis=0)]
    return {'radius_window_km': [lower,upper], 'points': len(x), 'minimum_searched_wavelength_km': minimum_wavelength, 'dominant_spatial_wavelength_km': float(1/frequency), 'frequency_grid_cycles_per_km': frequencies.tolist(), 'fit_improvement': score.tolist(), 'counting_noise_perturbation_10_90percent_km': np.quantile(wavelengths,[.1,.9]).tolist(), 'perturbation_trials': 256, 'seed': seed, 'reduced_chi_squared': float(np.sum(residual**2)/(len(x)-5)), 'scope': 'Dominant local single-harmonic component plus quadratic baseline. Strong nonlinear crests, overlapping waves, calibration and finite-source/diffraction effects are not described by this fit. Perturbation interval is counting-noise sensitivity, not a physical-wave confidence interval or surface-density inference.'}


def peak_spacing(radius, depth):
    smooth = np.convolve(depth,np.ones(3)/3,mode='same')
    candidates = [index for index in range(2,len(radius)-2) if 125295<=radius[index]<=125430 and smooth[index]>smooth[index-1] and smooth[index]>=smooth[index+1] and smooth[index]>.85]
    retained = []
    for index in sorted(candidates,key=lambda item:smooth[item],reverse=True):
        if all(abs(radius[index]-radius[other])>=7 for other in retained): retained.append(index)
    positions = sorted(float(radius[index]) for index in retained)
    return {'method': 'three-bin boxcar, local maximum above0.85 normal depth, greedily retain highest peaks separated by>=7km, declared125295..125430km window', 'crest_radius_km': positions, 'consecutive_spacing_km': np.diff(positions).tolist(), 'scope': 'Feature/threshold dependent descriptive spacing; one-km samples and three-km smoothing limit crest localization. A missed crest doubles a spacing. No claim that every crest belongs to one linear mode.'}


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--input',type=Path,default=DEFAULT_INPUT)
    parser.add_argument('--out',type=Path,default=ROOT/'docs/astronomy/ac10/cassini-replay.json' if REPOSITORY_INPUT.exists() else HERE/'cassini-replay-rerun.json')
    args = parser.parse_args()
    raw = args.input.read_bytes()
    if hashlib.sha256(raw).hexdigest()!=INPUT_HASH: raise ValueError('Frozen exact source-row hash changed')
    rows = np.loadtxt(args.input,delimiter=',')
    if rows.shape!=(1001,12) or not np.all(rows[:,11]==0): raise ValueError('Unexpected dimensions or quality flags')
    measured = derive(rows)
    maximum_difference = float(np.nanmax(np.abs(measured['normal_depth']-rows[:,4])))
    windows = [(125310,125390),(125410,125490),(125510,125590)]
    searches = [local_search(measured['radius_km'],measured['normal_depth'],measured['sigma_normal_depth_counting_only'],lower,upper,1010+index) for index,(lower,upper) in enumerate(windows)]
    # Coarsen actual counts and models with their numbers of raw integrations as weights.
    pairs = rows[:-1].reshape(-1,2,12)
    total_counts = pairs[:,:,10].sum(axis=1)
    coarsened = pairs[:,0,:].copy()
    for column in [0,3,8,9]: coarsened[:,column]=np.sum(pairs[:,:,column]*pairs[:,:,10],axis=1)/total_counts
    coarsened[:,10]=total_counts
    paired = derive(coarsened)
    pair_searches = [local_search(paired['radius_km'],paired['normal_depth'],paired['sigma_normal_depth_counting_only'],lower,upper,1110+index) for index,(lower,upper) in enumerate(windows)]
    result = {'kind':'measured processed Cassini UVIS binned photon profile, not a browser simulation','product_id':'UVIS_HSP_2008_261_ALPARA_E_TAU01KM.TAB','dataset_version':'3.0,2020','input_sha256':INPUT_HASH,'numpy_version':np.__version__,'shape':[1001,12],'radial_sample_km':1,'uv_band_nm':[110,190],'opening_angle_degrees':OPENING_DEGREES,'normal_depth_from_counts_formula':'-abs(sin B)*log((mean_signal-background_model)/star_model)','maximum_abs_archive_tau_difference':maximum_difference,'difference_explanation':'Printed54.4-degree opening angle, rounded per-bin counts/models and rounded archive tau; the unrounded original geometry is not in this binned cutout.','measured':{key:value.tolist() for key,value in measured.items()},'quality_flags':rows[:,11].astype(int).tolist(),'raw_samples_per_bin':rows[:,10].astype(int).tolist(),'detectable_tau_limit':rows[:,5].tolist(),'local_searches':searches,'two_km_rebin_searches':pair_searches,'peak_spacing':peak_spacing(measured['radius_km'],measured['normal_depth']),'failure_cases':['Negative net signal is censored, never logged or clipped into an optical depth.','Saturated or nonzero-flag bins require an explicit treatment; no fabricated replacement.','The one-harmonic local fit has large reduced chi-squared: narrow counting-only intervals omit model error.','A 1km radial grid cannot resolve structures near or below2km, and paired bins further suppress short structure.','Fixed110–190nm UV band is not a multi-color measurement; reported wavelength change is spatial.','This profile alone does not determine surface mass density, viscosity or ring age.']}
    args.out.parent.mkdir(parents=True,exist_ok=True)
    args.out.write_text(json.dumps(result,indent=2,allow_nan=False)+'\n')
    print(json.dumps({'max_tau_difference':maximum_difference,'local_searches':[{key:item[key] for key in ['radius_window_km','dominant_spatial_wavelength_km','counting_noise_perturbation_10_90percent_km','reduced_chi_squared']} for item in searches],'two_km_wavelengths':[item['dominant_spatial_wavelength_km']for item in pair_searches],'peak_spacing':result['peak_spacing']},indent=2))


if __name__=='__main__':main()