"""Illustrative one-channel direct-detection optical link; no hardware measurement.
Run python reproduce.py; requires matplotlib. Outputs are written beside script.
"""
import math,json,csv,random
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
ROOT=Path(__file__).resolve().parent
PARAMS=dict(seed=690301,bits_per_case=200000,wavelength_nm=1550,bitrate_bps=1e9,noise_bandwidth_hz=.7e9,transmitted_one_level_dbm=0,fixed_loss_db=3,attenuation_db_km=.2,extinction_ratio=10,responsivity_a_w=.8,dark_current_a=1e-8,load_ohm=50,electron_charge_c=1.602176634e-19,boltzmann_j_k=1.380649e-23,bit_one_probability=.5)
TEMPS=(273.15,298.15,333.15)
def state(distance,temperature):
    loss=3+.2*distance;p1=1e-3*10**(-loss/10);p0=p1/10;i1=.8*p1;i0=.8*p0
    s1=math.sqrt(2*1.602176634e-19*(i1+1e-8)*.7e9+4*1.380649e-23*temperature*.7e9/50)
    s0=math.sqrt(2*1.602176634e-19*(i0+1e-8)*.7e9+4*1.380649e-23*temperature*.7e9/50)
    threshold=(i1+i0)/2; margin=(i1-i0)/2
    ber=.25*(math.erfc(margin/s1/math.sqrt(2))+math.erfc(margin/s0/math.sqrt(2)))
    return dict(distance_km=distance,temperature_k=temperature,received_one_level_dbm=-loss,one_power_w=p1,zero_power_w=p0,one_current_a=i1,zero_current_a=i0,one_sigma_a=s1,zero_sigma_a=s0,threshold_a=threshold,analytic_ber=ber)
def csvwrite(name,rows):
    with (ROOT/name).open('w',newline='',encoding='utf-8') as out:
        w=csv.DictWriter(out,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
def main():
    rng=random.Random(PARAMS['seed']);cases=[];n=PARAMS['bits_per_case']
    for t in TEMPS:
        for d in (0,25,50,75,100,125):
            s=state(d,t);errors=0
            for _ in range(n):
                bit=rng.getrandbits(1)
                current=(s['one_current_a'] if bit else s['zero_current_a'])+rng.gauss(0,s['one_sigma_a'] if bit else s['zero_sigma_a'])
                errors+=int((current>=s['threshold_a'])!=bool(bit))
            p=s['analytic_ber'];tol=5*math.sqrt(p*(1-p)/n)+2/n
            cases.append(dict(**s,bits=n,errors=errors,observed_ber=errors/n,zero_error_upper95=1-.05**(1/n) if errors==0 else '',agreement_tolerance=tol,agreement_within_tolerance=abs(errors/n-p)<=tol))
    curves=[state(d/2,t) for t in TEMPS for d in range(301)]
    csvwrite('monte-carlo.csv',cases);csvwrite('analytic-curves.csv',curves)
    csvwrite('parameters.csv',[dict(parameter=k,value=v) for k,v in PARAMS.items()])
    assert all(c['agreement_within_tolerance'] for c in cases)
    assert all(state(d,t)['analytic_ber']<=state(d+1,t)['analytic_ber'] for t in TEMPS for d in range(150))
    metadata=dict(title='Illustrative optical on-off-keying link',evidence='Gaussian analytical model and seeded Monte Carlo; no hardware measurements',parameters=PARAMS,temperatures_k=TEMPS,curve_rows=len(curves),monte_carlo_rows=len(cases),equations=dict(loss_db='3 + 0.20*d_km',p1_w='1e-3 * 10^(-loss_db/10)',p0_w='p1_w/10',photocurrent='0.8 * optical_power_w',variance='2*q*(I+Idark)*B + 4*k*T*B/R',threshold='(I1+I0)/2',ber='0.5*[Q((I1-threshold)/sigma1)+Q((threshold-I0)/sigma0)]',Q='0.5*erfc(x/sqrt(2))',zero_count_upper95='1-0.05^(1/N)'),limitations=['Transmitted and received power values refer to logical-one levels, not average optical power.','Midpoint threshold and equally likely independent bits; unequal variances mean threshold is not asserted optimal.','The constant dark-current mean is removed; dark current contributes shot noise.','1550 nm identifies an illustrative band; wavelength is not a varied model input.','Temperature affects only Johnson noise.','No dispersion, clock recovery, wavelength multiplexing, laser intensity noise, nonlinearities, error correction or protocol.','All hardware values and attenuation are scenario inputs, not measured specifications or a cable guarantee.','Zero sampled errors is not zero error probability. One-sided bounds shown only for zero-count cases.','Very small analytical tails can numerically underflow to zero; the plot omits such points rather than claiming perfect transmission.'],checks=dict(all_18_monte_carlo_within_declared_tolerance=True,distance_monotonic=True),parity_samples=[state(75,298.15),state(100,298.15),state(125,333.15)],sources=['https://www.itu.int/rec/T-REC-G.652-202408-I/en'])
    (ROOT/'study.json').write_text(json.dumps(metadata,indent=2)+'\n',encoding='utf-8')
    plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'axes.grid':True,'grid.alpha':.2})
    fig,axes=plt.subplots(1,2,figsize=(10,4.6),layout='constrained')
    axes[0].plot([q['distance_km'] for q in curves if q['temperature_k']==298.15],[q['received_one_level_dbm'] for q in curves if q['temperature_k']==298.15],color='#285ec7')
    axes[0].set(title='A. Optical one-level decreases with distance',xlabel='Fibre length (km)',ylabel='Received logical-one power (dBm)')
    ax=axes[1]
    for t,color,offset in zip(TEMPS,('#13876d','#285ec7','#ad367f'),(-.8,0,.8)):
        data=[q for q in curves if q['temperature_k']==t and q['analytic_ber']>=1e-12]
        ax.plot([q['distance_km'] for q in data],[q['analytic_ber'] for q in data],color=color,label=f'{t-273.15:g} °C analytical')
        for q in [c for c in cases if c['temperature_k']==t]:
            x=q['distance_km']+offset
            if q['errors']:
                ax.plot(x,q['observed_ber'],'o',color=color,ms=4)
            else:
                bound=q['zero_error_upper95'];ax.plot(x,bound,'v',color=color,ms=4);ax.annotate('',xy=(x,bound/3),xytext=(x,bound),arrowprops={'arrowstyle':'->','color':color,'lw':.8})
    ax.set(yscale='log',ylim=(1e-12,1),xlim=(-3,153),title='B. Errors rise as current levels overlap',xlabel='Fibre length (km)',ylabel='Bit error probability');ax.legend(fontsize=7,loc='lower right')
    ax.text(.03,.98,'Dots: sampled error rate\n▼ arrows: zero errors, 95% upper bound\n200,000 bits/case; markers offset ≤0.8 km',transform=ax.transAxes,va='top',fontsize=7)
    fig.suptitle('Single-channel optical model • scenario inputs, not a measured link',fontsize=11)
    for ext in ('png','svg','pdf'):fig.savefig(ROOT/f'optical-study.{ext}',dpi=200,bbox_inches='tight')
    plt.close(fig)
    print(json.dumps(metadata['parity_samples'],indent=2))
if __name__=='__main__':main()