"""Deterministic timebase sensitivity and illustrative spiral geometry."""
import json,math,csv
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
FS=48000
DURATION_S=2.0
ROTATIONS_PER_SECOND=2.0
CUTOFF_HZ=20000.0

def tone(frequency: float) -> list[float]:
    n = int(FS * DURATION_S)
    return [0.8 * math.sin(2.0 * math.pi * frequency * i / FS) for i in range(n)]

def read_with_linear_track_drift(source: list[float], end_phase_error_deg: float) -> list[float]:
    n = len(source)
    end_delay = end_phase_error_deg / 360.0 / ROTATIONS_PER_SECOND
    delay_samples = end_delay * FS
    result = []
    for i in range(n):
        position = i - delay_samples * i / max(1, n - 1)
        position = min(max(position, 0.0), n - 1.0)
        left = int(position)
        right = min(left + 1, n - 1)
        fraction = position - left
        result.append(source[left] * (1.0 - fraction) + source[right] * fraction)
    return result

def one_pole_lowpass(values: list[float], cutoff_hz: float) -> list[float]:
    a = 1.0 - math.exp(-2.0 * math.pi * cutoff_hz / FS)
    out = []
    previous = 0.0
    for value in values:
        previous += a * (value - previous)
        out.append(previous)
    return out

def metrics(reference: list[float], output: list[float], frequency: float) -> dict:
    start = int(0.1 * FS)
    x = reference[start:]
    y = output[start:]
    n = len(y)
    omega = 2.0 * math.pi * frequency / FS
    ss = cc = sc = sy = cy = 0.0
    for i, value in enumerate(y, start=start):
        s = math.sin(omega * i)
        c = math.cos(omega * i)
        ss += s * s
        cc += c * c
        sc += s * c
        sy += value * s
        cy += value * c
    det = ss * cc - sc * sc
    a = (sy * cc - cy * sc) / det
    b = (cy * ss - sy * sc) / det
    fitted = [a * math.sin(omega * i) + b * math.cos(omega * i) for i in range(start, len(output))]
    signal_power = sum((v * v for v in fitted)) / n
    residual_power = sum(((v - h) ** 2 for v, h in zip(y, fitted))) / n
    mx = sum(x) / n
    my = sum(y) / n
    covariance = sum(((u - mx) * (v - my) for u, v in zip(x, y)))
    vx = sum(((u - mx) ** 2 for u in x))
    vy = sum(((v - my) ** 2 for v in y))
    corr = covariance / math.sqrt(vx * vy) if vx and vy else 0.0
    snr = 10.0 * math.log10(signal_power / max(residual_power, 1e-30))
    thdn = 100.0 * math.sqrt(residual_power / max(signal_power, 1e-30))
    return {'correlation': corr, 'fitted_tone_signal_to_residual_dB': snr, 'relative_residual_RMS_percent': thdn}
cases=[]
for phase in [0,.1,1,10]:
 for freq in [100,1000,10000]:
  source=tone(freq);output=one_pole_lowpass(read_with_linear_track_drift(source,phase),CUTOFF_HZ)
  score=metrics(source,output,freq);delay=phase/(360*ROTATIONS_PER_SECOND)
  cases.append(dict(scanner_endpoint_phase_deg=phase,tone_Hz=freq,endpoint_time_lag_s=delay,audio_endpoint_phase_lag_deg=360*freq*delay,**score))
spiral=[]
for pitch in [0,.002]:
 for k in range(801):
  theta=8*math.pi*k/800;r=.01+pitch*theta/(2*math.pi)
  spiral.append(dict(pitch_m_per_turn=pitch,theta_rad=theta,radius_m=r,x_m=r*math.cos(theta),y_m=r*math.sin(theta),time_s=theta/(2*math.pi*ROTATIONS_PER_SECOND)))
for name,rows in [('timebase-cases.csv',cases),('spiral.csv',spiral)]:
 with (P/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
identity=max(abs(r['relative_residual_RMS_percent']-100*10**(-r['fitted_tone_signal_to_residual_dB']/20)) for r in cases);assert identity<1e-10
assert all(abs(math.hypot(r['x_m'],r['y_m'])-r['radius_m'])<1e-15 for r in spiral)
assert all(r['radius_m']==.01 for r in spiral if r['pitch_m_per_turn']==0)
assert abs(spiral[-1]['radius_m']-.018)<1e-15
assert len(cases)==12 and sum(r['scanner_endpoint_phase_deg']==0 for r in cases)==3
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','pdf.fonttype':42})
fig,ax=plt.subplots(1,2,figsize=(11,4.6))
for freq,color in [(100,'#087e91'),(1000,'#a14786'),(10000,'#7451a0')]:
 rr=[r for r in cases if r['tone_Hz']==freq];ax[0].plot(range(4),[r['fitted_tone_signal_to_residual_dB'] for r in rr],'-o',color=color,label=f'{freq:g} Hz tone')
 ax[1].plot([r['endpoint_time_lag_s']*1000 for r in rr],[r['audio_endpoint_phase_lag_deg'] for r in rr],'-o',color=color,label=f'{freq:g} Hz tone')
ax[0].set(xticks=range(4),xticklabels=['0','0.1','1','10'],xlabel='Accumulated scanner phase error (degrees)',ylabel='Fitted-tone signal-to-residual ratio (dB)',title='A  |  Drift spoils a fit at the original frequency');ax[0].legend(fontsize=8)
ax[1].set(xlabel='Endpoint time lag (ms)',ylabel='Unwrapped audio phase lag (degrees)',title='B  |  Scanner phase and audio phase differ')
for a in ax:a.grid(alpha=.2)
fig.suptitle('Deterministic timing drift in twelve synthetic tone cases',fontsize=14);fig.text(.5,.05,'48 kHz sampling, 2 s duration, 2 scanner turns/s. No simulated noise. Zero-drift residuals are dominated by numerical precision.',ha='center',fontsize=8);fig.text(.5,.015,'The ratio includes frequency mismatch, not just noise or harmonic distortion. These cases are not hardware acceptance trials.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.09,1,.93])
for ext in ['png','svg','pdf']:fig.savefig(P/f'timebase-drift.{ext}',dpi=190)
fig,ax=plt.subplots(1,2,figsize=(10,4.6))
for a,pitch,title in zip(ax,[0,.002],['A  |  Zero pitch: a circle','B  |  Changing radius: a spiral']):
 rr=[r for r in spiral if r['pitch_m_per_turn']==pitch];a.plot([r['x_m']*1000 for r in rr],[r['y_m']*1000 for r in rr],color='#087e91');a.set(xlabel='x (mm)',ylabel='y (mm)',title=title);a.set_aspect('equal');a.grid(alpha=.2)
fig.suptitle('Quadrature coordinates plus a radius law',fontsize=14);fig.text(.5,.05,'r = 10 mm + pitch × turns. Four turns, pitch 0 or 2 mm/turn; these illustrative paths do not simulate a physical recording medium.',ha='center',fontsize=8);fig.text(.5,.015,'At fixed angular speed the spot speed changes with radius. Pitch alone does not establish storage capacity or readable resolution.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.09,1,.93])
for ext in ['png','svg','pdf']:fig.savefig(P/f'spiral-geometry.{ext}',dpi=190)
meta={'title':'Deterministic timebase benchmark and separate spiral geometry','status':'Reproduced twelve synthetic tone calculations; no hardware, noise source or recording medium','rows':{'tone_cases':12,'spiral':1602},'parameters':{'sampling_Hz':FS,'duration_s':2,'tone_amplitude':.8,'tones_Hz':[100,1000,10000],'scanner_turns_per_s':2,'endpoint_scanner_phase_deg':[0,.1,1,10],'one_pole_coefficient':1-math.exp(-2*math.pi*CUTOFF_HZ/FS),'coefficient_frequency_parameter_Hz':20000,'discard_initial_s':.1,'spiral_start_radius_m':.01,'spiral_pitch_m_per_turn':[0,.002],'spiral_turns':4},'metric_definitions':{'fitted_tone_signal_to_residual_dB':'10log10(mean(fit²)/mean((output-fit)²)); fit is sine/cosine at original tone frequency','relative_residual_RMS_percent':'100sqrt(mean((output-fit)²)/mean(fit²))','correlation':'Pearson input/output correlation, phase-sensitive and not a listening score'},'checks':{'all_pass':True,'ratio_residual_identity_max_absolute_error':identity,'circle_and_spiral_radius_identity':True,'twelve_shared_cases':True},'limits':['Scanner phase accumulates linearly over the finite track, causing deterministic time stretching; it is not random jitter.','Sample positions use linear interpolation with endpoint clamping; the first 0.1 s is discarded before fitting.','The 20 kHz parameter sets a discrete one-pole coefficient; it is not an exact digital -3 dB cutoff.','No random noise is injected. Fitted residual includes original-frequency mismatch; zero-drift tiny residuals primarily reflect numerical precision.','The two residual metrics are algebraically dependent; neither is independently measured noise or harmonic distortion.','The same twelve cases support recording and scanning lessons, not two independent experiments. No hardware pass/fail rate is claimed.','Spiral geometry is a separate illustration, not the physical trajectory used in the timebase simulation. No spot, optics, emulsion, servo, cross-talk or medium bandwidth is modeled.']}
(P/'study.json').write_text(json.dumps(meta,indent=2)+'\n');print(meta['checks'])