"""Small-angle rays in an imposed linear refractive-index gradient."""
import csv,json,math
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
N0=1.00027

def ray(x,g,y0=0,slope0=0):
 if not all(math.isfinite(v) for v in [x,g,y0,slope0]) or not 0<=x<=2 or abs(g)>.002 or abs(y0)>.01 or abs(slope0)>.005:raise ValueError('Inputs outside the bounded small-angle example')
 y=y0+slope0*x+g*x*x/(2*N0);slope=slope0+g*x/N0
 return dict(x_m=x,gradient_per_m=g,y_m=y,slope_dy_dx=slope,index_at_ray=N0+g*y)
rows=[ray(i/100,g) for g in [-.002,0,.002] for i in range(201)]
temp=[dict(temperature_K=t,index_fixed_pressure=1+.00027*293/t) for t in range(273,334)]
# Explicit Euler as a deliberately first-order independent integration method.
conv=[]
for steps in [10,20,40,80,160,320]:
 dx=2/steps;y=0;s=0;g=.002
 for _ in range(steps):y,s=y+s*dx,s+g/N0*dx
 exact=ray(2,g)
 conv.append(dict(steps=steps,step_m=dx,euler_endpoint_y_m=y,analytic_endpoint_y_m=exact['y_m'],absolute_error_m=abs(y-exact['y_m']),endpoint_slope_error=abs(s-exact['slope_dy_dx'])))
assert all(abs(conv[i]['absolute_error_m']/conv[i+1]['absolute_error_m']-2)<1e-9 for i in range(5))
assert max(r['endpoint_slope_error'] for r in conv)<1e-15
# Independent finite-difference acceleration from sampled analytic curve.
fd_error=0
positive=[r for r in rows if r['gradient_per_m']==.002]
for a,b,c in zip(positive,positive[1:],positive[2:]):
 fd=(c['y_m']-2*b['y_m']+a['y_m'])/.01**2;fd_error=max(fd_error,abs(fd-.002/N0))
assert fd_error<1e-12
assert all(abs(rows[i]['y_m']+rows[402+i]['y_m'])<1e-15 for i in range(201))
assert all(r['y_m']==0 for r in rows[201:402])
assert all(r['index_at_ray']>0 and abs(r['slope_dy_dx'])<.0041 for r in rows)
assert 1+.00027*293/293==N0
assert all(a['index_fixed_pressure']>b['index_fixed_pressure'] for a,b in zip(temp,temp[1:]))
for name,data in [('rays.csv',rows),('temperature-index.csv',temp),('convergence.csv',conv)]:
 with (P/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=data[0]);w.writeheader();w.writerows(data)
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=(10.5,4.8))
for g,color in [(-.002,'#a04b85'),(0,'#777'),(.002,'#087e91')]:
 rr=[r for r in rows if r['gradient_per_m']==g];ax[0].plot([r['x_m'] for r in rr],[r['y_m']*1000 for r in rr],label=f'g = {g:+.3f} per metre',color=color,lw=2)
ax[0].set(xlabel='Forward distance x (m)',ylabel='Transverse position y (mm)',title='A  |  A ray bends toward higher refractive index');ax[0].legend(fontsize=8);ax[0].text(.04,.93,'Vertical scale exaggerated',transform=ax[0].transAxes,fontsize=8)
ax[1].loglog([r['step_m'] for r in conv],[r['absolute_error_m']*1000 for r in conv],'-o',color='#75469c');ax[1].set(xlabel='Euler step length (m)',ylabel='Endpoint error relative to analytic ray (mm)',title='B  |  Halving the step halves this solver error')
for a in ax:a.grid(alpha=.2)
fig.suptitle('An imposed index gradient: a bounded small-angle model',fontsize=15)
fig.text(.5,.055,'n(y) = 1.00027 + g y; approximate ray equation y\'\' = g / 1.00027. Curves solve this approximate equation analytically.',ha='center',fontsize=9)
fig.text(.5,.015,'No measured medium, temperature field, lens design or atmospheric prediction. Numerical convergence does not validate the physical approximation.',ha='center',fontsize=8)
fig.tight_layout(rect=[0,.10,1,.92])
for ext in ['png','svg','pdf']:fig.savefig(P/f'gradient-rays.{ext}',dpi=190)
meta={'title':'Paraxial rays in a prescribed linear transverse index','status':'Hypothetical analytical illustration and independent numerical integration','rows':{'rays':603,'temperature_index':61,'convergence':6},'parameters':{'n0':N0,'gradient_per_m':[-.002,0,.002],'distance_m':[0,2],'initial_y_m':0,'initial_slope':0,'temperature_reference_K':293,'index_reference':N0,'temperature_range_K':[273,333]},'equations':{'index':'n(y)=n0+g*y','approximate_ODE':'d²y/dx²=g/n0','analytic_reference':'y=y0+slope0*x+g*x²/(2*n0); dy/dx=slope0+g*x/n0','temperature_index':'n(T)-1=(nref-1)*Tref/T at fixed pressure/composition'},'checks':{'all_pass':True,'finite_difference_acceleration_max_abs_error_per_m':fd_error,'Euler_first_order_error_halving':True,'gradient_sign_symmetry':True,'zero_gradient_straight_ray':True,'index_positive':True,'temperature_reference_identity':True,'temperature_index_monotone':True},'endpoint_positive_gradient':ray(2,.002),'limitations':['The analytic reference exactly solves the stated approximate ODE, not the full eikonal ray equation.','Assumes small slope and small index variation so n in the denominator can be replaced by n0; curvature terms and higher-order corrections are omitted.','The prescribed linear gradient is not a computed temperature, density or flow field. No turbulence, diffraction, boundaries or material dispersion is included.','The optional temperature relation assumes weak refractivity proportional to density, ideal-gas density at fixed pressure/composition, and fixed optical wavelength. Reference index is illustrative, not an air refractivity standard.','Temperature relation does not specify a spatial temperature gradient and is not used to generate the rays.','Solver convergence verifies numerical implementation only; it cannot establish device feasibility or atmospheric accuracy.'],'applications':['Established context: understanding gradient-index optics and qualitative ray bending in nonuniform media.','Potential research use: a first sensitivity calculation before measured index maps, wave optics and device validation; no improvement or performance is demonstrated.']}
(P/'study.json').write_text(json.dumps(meta,indent=2)+'\n');print(json.dumps(meta['checks']))