from pathlib import Path
import math,csv,json
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
OUT=Path(__file__).resolve().parent
C=299792458.;E=1.602176634e-19;EPS=8.8541878188e-12;ME=9.1093837139e-31

def loss(dkm,fmhz):return 20*math.log10(4*math.pi*dkm*1000*fmhz*1e6/C)
def fp(ne):return math.sqrt(ne*E*E/(EPS*ME))/(2*math.pi)
def secant(fc,deg):
    if not 0<=deg<=75:raise ValueError('Teaching domain: 0..75 degree launch zenith angles in the sub-ionospheric n~1 region only.')
    return fc/math.cos(math.radians(deg))
def write(name,rows):
    with (OUT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
path=[dict(distance_km=float(d),frequency_MHz=f,free_space_loss_dB=loss(d,f),vacuum_one_way_delay_ms=d*1e6/C) for d in np.geomspace(1,100000,101) for f in [10,100,1000]]
plasma=[dict(electron_density_m3=float(ne),critical_frequency_MHz=fp(ne)/1e6) for ne in np.geomspace(1e10,1e12,101)]
angles=[dict(vertical_critical_frequency_MHz=fc,launch_zenith_angle_deg=a,secant_multiplier=1/math.cos(math.radians(a)),ideal_oblique_limit_MHz=secant(fc,a)) for a in range(76) for fc in [3,6,9]]
write('free-space.csv',path);write('plasma-frequency.csv',plasma);write('secant-law.csv',angles)
# Independent derivation: unit isotropic transmitter flux P/(4*pi*d²),
# effective receiving aperture lambda²/(4*pi). Reciprocal received fraction is loss.
errors=[]
for r in path:
    d=r['distance_km']*1000;lam=C/(r['frequency_MHz']*1e6)
    received=1/(4*math.pi*d*d)*lam*lam/(4*math.pi)
    errors.append(abs(-10*math.log10(received)-r['free_space_loss_dB']))
assert max(errors)<1e-12
assert abs(loss(2,100)-loss(1,100)-20*math.log10(2))<1e-12
assert abs(loss(1,200)-loss(1,100)-20*math.log10(2))<1e-12
assert abs(fp(4e11)/fp(1e11)-2)<1e-14
assert secant(6,0)==6 and abs(secant(6,60)-12)<1e-12
assert abs(299792.458*1e6/C-1000)<1e-10
# Rounded engineering plasma rule 8.98*sqrt(ne[m^-3]) Hz.
relative=abs(8.98*math.sqrt(1e12)/fp(1e12)-1);assert relative<.0002
for a in [-1,76,90]:
    try:secant(6,a)
    except ValueError:pass
    else:raise AssertionError('Angle guard missing')
meta={'title':'Radio paths: distance, delay and an ideal plasma limit','kind':'Analytical teaching sweeps, not measured reception or operational forecasts','constants':{'c_m_s':C,'elementary_charge_C':E,'vacuum_permittivity_F_m':EPS,'electron_mass_kg':ME},'formulas':{'free_space_loss_dB':'20 log10(4*pi*d*f/c), SI units','delay_s':'d/c in vacuum','critical_frequency_Hz':'sqrt(ne*e²/(epsilon0*me))/(2*pi)','secant':'f_oblique_limit=fc/cos(theta), theta is the launch zenith angle in the sub-ionospheric region (refractive index approximately 1), measured from the upward layer normal; not the local angle at the turning point'},'row_counts':{'free_space':len(path),'plasma_frequency':len(plasma),'secant':len(angles)},'domains':{'distance_km':[1,100000],'frequency_MHz':[10,100,1000],'electron_density_m3':[1e10,1e12],'launch_zenith_angle_deg':[0,75],'vertical_fc_MHz':[3,6,9]},'validation':{'isotropic_flux_aperture_max_difference_dB':max(errors),'doubling_distance_and_frequency_add_6_0206dB':True,'quadruple_density_doubles_plasma_frequency':True,'angle0_multiplier1_angle60_multiplier2':True,'one_light_second_distance_gives1000ms':True,'rounded_8_98_plasma_rule_relative_difference':relative,'out_of_domain_angles_rejected':True},'limitations':['Free-space model requires an unobstructed far-field path, matched polarisation and antenna gains handled separately. Distances are hypothetical path lengths, not ground coverage.','Frequency comparison holds dimensionless antenna gains fixed; fixed physical antenna apertures can change the frequency dependence.','Delay is vacuum one-way travel time only, excluding equipment, processing and dispersive plasma group delay.','Plasma model is cold, collisionless, unmagnetized with stationary ions and prescribed electron density.','Critical frequency uses peak layer density for the ideal vertical limit.','Secant relation assumes plane stratified geometry. Theta is the sub-ionospheric launch zenith angle at refractive index approximately 1, not local turning-point incidence; 75 degrees is a teaching cap, not a certified validity boundary. Near-horizon extrapolation diverges and cannot represent curved Earth and finite ionosphere.','No ray paths, geographic ranges, absorption, magnetic splitting, weather variation or guaranteed usable frequencies are computed.'],'sources':[{'title':'ITU-R P.525-5 (2024), free-space attenuation','url':'https://www.itu.int/dms_pubrec/itu-r/rec/p/R-REC-P.525-5-202411-I!!PDF-E.pdf'},{'title':'Richard Fitzpatrick, University of Texas, Plasma Frequency eq1.5','url':'https://farside.ph.utexas.edu/teaching/plasma/Plasma/node5.html'},{'title':'NBS Monograph80, Ionospheric Radio Propagation, section4.1 plane Earth and plane ionosphere, secant law eq4.2b','url':'https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nbsmonograph80.pdf'},{'title':'NIST CODATA2022 constants','url':'https://physics.nist.gov/cuu/pdf/all.pdf'}]}
(OUT/'study.json').write_text(json.dumps(meta,indent=2)+'\n',encoding='utf-8')
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False})
fig,axs=plt.subplots(2,2,figsize=(11,8),constrained_layout=True);a,b,c,d=axs.ravel();colors=['#393f94','#087f8c','#b42e8d']
for f,col in zip([10,100,1000],colors):
    r=[x for x in path if x['frequency_MHz']==f];a.semilogx([x['distance_km'] for x in r],[x['free_space_loss_dB'] for x in r],label=f'{f} MHz',color=col)
r=path[::3];b.loglog([x['distance_km'] for x in r],[x['vacuum_one_way_delay_ms'] for x in r],color=colors[0])
c.loglog([x['electron_density_m3'] for x in plasma],[x['critical_frequency_MHz'] for x in plasma],color=colors[1])
for fc,col in zip([3,6,9],colors):
    r=[x for x in angles if x['vertical_critical_frequency_MHz']==fc];d.plot([x['launch_zenith_angle_deg'] for x in r],[x['ideal_oblique_limit_MHz'] for x in r],label=f'fc = {fc} MHz',color=col)
a.set(title='A  |  Geometric spreading',xlabel='Hypothetical path length (km)',ylabel='Free-space loss (dB)');a.legend(fontsize=8)
b.set(title='B  |  A finite travel time',xlabel='Vacuum path length (km)',ylabel='One-way delay (ms)')
c.set(title='C  |  More electrons, higher critical frequency',xlabel='Peak electron density (m⁻³)',ylabel='Ideal critical frequency (MHz)')
d.set(title='D  |  The bounded planar secant model',xlabel='Launch zenith angle below layer (degrees)',ylabel='Ideal oblique frequency limit (MHz)',xlim=(0,75));d.legend(fontsize=8)
for ax in axs.ravel():ax.grid(alpha=.2,which='both')
fig.suptitle('Radio propagation: useful models, different assumptions',fontsize=16,weight='bold')
fig.supxlabel('Analytical examples only. Panels A–B are vacuum references; C–D use an ideal plasma layer, not a forecast of reception.',fontsize=9)
for ext in ['png','svg','pdf']:fig.savefig(OUT/f'radio-paths.{ext}',dpi=180)
print(json.dumps({'rows':meta['row_counts'],'validation':meta['validation'],'example_loss100MHz1km':loss(1,100),'example_fp1e12_MHz':fp(1e12)/1e6}))