"""Deterministic teaching models, not measurements or speech recognition scores.
Run: python reproduce.py. Requires matplotlib; Python standard library otherwise.
"""
import csv,json,math
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
ROOT=Path(__file__).resolve().parent
C=343.0

def noise_level(background,talker):
    return 10*math.log10(10**(background/10)+10**(talker/10))
def residual_ratio(f,x,r=1.0,phase_deg=0.0):
    theta=2*(2*math.pi*f/C)*x+math.radians(phase_deg)
    # Stable equivalent of 1+r*r-2*r*cos(theta).
    return math.sqrt((1-r)**2+4*r*math.sin(theta/2)**2)
def level(ratio):
    return None if ratio==0 else 20*math.log10(ratio)
def dump_csv(name,rows):
    with (ROOT/name).open('w',newline='',encoding='utf-8') as out:
        writer=csv.DictWriter(out,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows)
def save(fig,name):
    for suffix in ('png','svg','pdf'):
        fig.savefig(ROOT/f'{name}.{suffix}',dpi=200,bbox_inches='tight')
    plt.close(fig)

def main():
    plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'axes.grid':True,'grid.alpha':.2,'figure.facecolor':'white'})
    snr=[]
    base_noise=noise_level(50,50)
    for i in range(61):
        gain=i*.25
        snr.append({'gain_db':gain,'voice_reference_db':60,'background_reference_db':50,'talker_reference_db':50,'combined_noise_reference_db':base_noise,'source_only_snr_db':60+gain-base_noise,'postmix_voice_db':60+gain,'postmix_noise_db':base_noise+gain,'postmix_snr_db':(60+gain)-(base_noise+gain)})
    noise=[]
    for i in range(161):
        talker=30+i*.25
        total=noise_level(50,talker)
        noise.append({'background_db':50,'talker_db':talker,'combined_noise_db':total,'voice_db':60,'snr_db':60-total})
    anc=[]
    for f in (125,250,500,1000):
        for i in range(2001):
            x=-.5+i*.0005
            ratio=residual_ratio(f,x)
            db=level(ratio)
            anc.append({'frequency_hz':f,'position_m':x,'amplitude_ratio_secondary_primary':1,'phase_error_deg':0,'pressure_rms_ratio':ratio,'level_change_db':db,'exact_null':ratio==0,'plot_level_change_db':max(-60,db) if db is not None else -60})
    errors=[]
    for r in (.8,.9,1,1.1,1.2):
        for phase in range(-90,91):
            ratio=residual_ratio(250,0,r,phase)
            db=level(ratio)
            errors.append({'frequency_hz':250,'position_m':0,'amplitude_ratio_secondary_primary':r,'phase_error_deg':phase,'pressure_rms_ratio':ratio,'level_change_db':db,'exact_null':ratio==0,'plot_level_change_db':max(-60,db) if db is not None else -60})
    for name,rows in [('snr-gain.csv',snr),('noise-sum.csv',noise),('anc-position.csv',anc),('anc-errors.csv',errors)]:dump_csv(name,rows)
    checks={
        'equal_noise_addition_db':noise_level(50,50)-50,
        'equal_noise_expected_db':10*math.log10(2),
        'postmix_snr_max_change_db':max(abs(row['postmix_snr_db']-snr[0]['postmix_snr_db']) for row in snr),
        'null_pressure_ratio':residual_ratio(250,0),
        'quarter_wavelength_level_db':level(residual_ratio(250,C/250/4)),
        'quarter_wavelength_expected_db':20*math.log10(2),
        'no_secondary_ratio':residual_ratio(250,.123,0),
        'phase_180_level_db':level(residual_ratio(250,0,1,180)),
    }
    # Mean-square sum includes coherent cross term; energy need not add locally.
    # At equal opposite waves, spatial average over integer half-wavelengths
    # of squared pressure ratio = 1+r^2, despite local nulls and maxima.
    n=10000;r=.8;f=250;wavelength=C/f
    spatial_mean=sum(residual_ratio(f,j*wavelength/n,r)**2 for j in range(n))/n
    checks['spatial_mean_pressure_ratio_squared']=spatial_mean
    checks['spatial_mean_expected']=1+r*r
    checks['frequency_position_scaling_error']=abs(residual_ratio(125,.2)-residual_ratio(1000,.025))
    assert abs(checks['equal_noise_addition_db']-checks['equal_noise_expected_db'])<1e-12
    assert checks['postmix_snr_max_change_db']<1e-12
    assert checks['null_pressure_ratio']==0
    assert abs(checks['quarter_wavelength_level_db']-checks['quarter_wavelength_expected_db'])<1e-12
    assert abs(spatial_mean-(1+r*r))<1e-12
    assert checks['no_secondary_ratio']==1
    assert checks['frequency_position_scaling_error']<1e-12
    metadata={
        'title':'Illustrative received SNR and narrowband coherent wave interference',
        'evidence_status':'New deterministic analytical teaching study; no physical measurements, listening tests or clinical outcomes.',
        'snr_parameters':{'voice_level_db':60,'background_level_db':50,'competing_talker_level_db':50,'level_convention':'Illustrative unweighted band-matched received sound pressure levels with common reference. Noise components assumed mutually uncorrelated. No absolute calibration or safety implication.'},
        'anc_parameters':{'sound_speed_m_s':C,'frequencies_hz':[125,250,500,1000],'geometry':'Two oppositely travelling 1D plane waves in a lossless medium; equal pressure amplitudes for spatial plot; secondary antiphase at x=0.','plot_floor_db':-60,'null_representation':'JSON null / empty CSV level_change_db with exact_null=true means negative infinity, not missing measurement.'},
        'equations':{'noise':'10 log10(10^(Lb/10)+10^(Lt/10))','snr':'Lvoice-Lnoise','waves':'p1=A cos(omega*t-k*x); p2=-r*A cos(omega*t+k*x+epsilon)','rms_ratio':'sqrt((1-r)^2+4*r*sin((2*k*x+epsilon)/2)^2)','level_change':'20 log10(rms_ratio)','mean_square_sum':'1+r^2-2*r*cos(2*k*x+epsilon); spatial mean over full wavelength = 1+r^2'},
        'limitations':['SNR is not a speech recognition or intelligibility score.','Common gain invariance assumes no clipping, quantization or added noise.','Two coherent waves cannot be added as independent powers locally.','The cancellation plot is single-frequency and one-dimensional, with prescribed waves, not a room, active controller, broadband field, energy-flux model or clinical device.','Plane waves have no geometric spreading; reflecting walls, adaptation, latency and transducer limits are absent.','The -60 dB display floor is graphical only, not measured attenuation.'],
        'row_counts':{'snr_gain':len(snr),'noise_sum':len(noise),'anc_position':len(anc),'anc_errors':len(errors)},'checks':checks,
        'sources':['https://www.nist.gov/itl/iad/mltg/nist-speech-signal-noise-ratio-measurements','https://resource.isvr.soton.ac.uk/spcg/tutorial/tutorial/Tutorial_files/Web-inter-applications.htm','https://www.osha.gov/otm/section-3-health-hazards/chapter-5']}
    (ROOT/'study.json').write_text(json.dumps(metadata,indent=2)+'\n',encoding='utf-8')
    fig,ax=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    ax[0].plot([q['gain_db'] for q in snr],[q['source_only_snr_db'] for q in snr],label='Increase voice only',color='#285ec7')
    ax[0].plot([q['gain_db'] for q in snr],[q['postmix_snr_db'] for q in snr],label='Amplify whole mixture',color='#ad367f',linestyle='--')
    ax[0].set(xlabel='Applied increase (dB)',ylabel='Received signal-to-noise ratio (dB)',title='A. Where gain is applied matters');ax[0].legend(fontsize=8)
    ax[1].plot([q['talker_db'] for q in noise],[q['combined_noise_db'] for q in noise],color='#285ec7',label='Combined noise')
    ax[1].axhline(50,color='gray',linestyle=':',label='Fixed background: 50 dB')
    ax[1].plot([30,70],[30,70],color='#ad367f',linestyle='--',label='Competing talker alone')
    ax[1].set(xlabel='Competing talker level (dB)',ylabel='Noise level (dB)',title='B. Independent noise powers add');ax[1].legend(fontsize=8)
    fig.suptitle('Received-level model • illustrative inputs, not listening-test data',fontsize=11)
    save(fig,'snr-study')
    fig,ax=plt.subplots(1,2,figsize=(10,4.4),layout='constrained')
    for f,color in zip((125,250,500,1000),('#285ec7','#ad367f','#13876d','#bd6a11')):
        rows=[q for q in anc if q['frequency_hz']==f]
        ax[0].plot([q['position_m'] for q in rows],[q['plot_level_change_db'] for q in rows],label=f'{f} Hz',color=color,lw=1.2)
    ax[0].set(xlabel='Position relative to tuned null (m)',ylabel='Level change from primary alone (dB)',title='A. Cancellation depends on location',ylim=(-62,8));ax[0].legend(fontsize=8,ncol=2)
    for r,color in zip((.8,1,1.2),('#285ec7','#ad367f','#13876d')):
        rows=[q for q in errors if q['amplitude_ratio_secondary_primary']==r]
        ax[1].plot([q['phase_error_deg'] for q in rows],[q['plot_level_change_db'] for q in rows],label=f'Amplitude ratio {r:g}',color=color)
    ax[1].set(xlabel='Phase error at x = 0 (degrees)',ylabel='Level change (dB)',title='B. Matching phase and amplitude',ylim=(-62,8));ax[1].legend(fontsize=8)
    for a in ax:a.axhline(0,color='black',lw=.7,linestyle=':')
    fig.suptitle('Opposite-travelling plane waves • exact nulls drawn at −60 dB floor',fontsize=11)
    save(fig,'anc-study')
    print(json.dumps({'rows':metadata['row_counts'],'checks':checks},indent=2))
if __name__=='__main__':main()