"""Illustrative single spherical liquid interface; no measured material data."""
import csv,json,math
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
ROOT=Path(__file__).resolve().parent
GAMMA=.03 # N/m, hypothetical interfacial tension
N1,N2=1.33,1.50 # hypothetical transparent media at one unspecified wavelength
DRHO=200. # kg/m3, hypothetical absolute density contrast
G=9.80665 # m/s2, standard gravity
APERTURE=.00005 # m, radius used for curvature sweep

def interface(curvature):
    return 2*GAMMA*curvature,(N2-N1)*curvature
def bond(a):return DRHO*G*a*a/GAMMA
def write(name,rows):
    with (ROOT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=rows[0].keys());w.writeheader();w.writerows(rows)

rows=[]
for i in range(-100,101):
    c=i*10.;p,power=interface(c)
    rows.append(dict(signed_curvature_per_m=c,signed_radius_m=1/c if c else '',pressure_p2_minus_p1_Pa=p,single_interface_power_per_m=power,aperture_radius_m=APERTURE,aperture_over_abs_radius=abs(APERTURE*c),bond_number=bond(APERTURE)))
apertures=[dict(aperture_radius_m=i*.00005,bond_number=bond(i*.00005)) for i in range(1,101)]
write('interface-curve.csv',rows);write('aperture-curve.csv',apertures)

# Independent checks: hemisphere pressure-force balance, exact Snell rays in
# the paraxial limit, odd symmetry, index match, and quadratic gravity scaling.
checks=[]
for c in (-1000.,-100.,100.,1000.):
    p,power=interface(c);r=abs(1/c)
    assert math.isclose(abs(p)*math.pi*r*r,2*math.pi*r*GAMMA,rel_tol=1e-12)
    h=1e-7/abs(c)
    normal=-math.asin(h*c)
    refracted=normal+math.asin(N1/N2*math.sin(-normal))
    snell_power=-N2*refracted/h
    assert math.isclose(power,snell_power,rel_tol=1e-12)
    assert interface(-c)==(-p,-power)
    checks.append(dict(curvature_per_m=c,power_per_m=power,independent_snell_limit_per_m=snell_power,hemisphere_force_balance_pass=True))
assert interface(0)==(0,0)
assert math.isclose(bond(.002)/bond(.001),4)
assert (N1-N1)*1000==0
metadata=dict(status='new illustrative analytical study; no experimental measurements',parameters=dict(interfacial_tension_N_per_m=GAMMA,refractive_index_incident=N1,refractive_index_transmitted=N2,absolute_density_contrast_kg_per_m3=DRHO,gravity_m_per_s2=G,curvature_sweep_aperture_radius_m=APERTURE),sign_convention='Light travels left to right from medium 1 to medium 2. R is positive when the centre of curvature is to the right of the vertex. Pressure difference is p2-p1.',equations=dict(pressure='2*gamma/R',interface_power='(n2-n1)/R',bond='abs(rho2-rho1)*g*a^2/gamma'),limits=['One spherical interface, not a complete lens or its focal length.','Indices, density contrast and tension are hypothetical; no real liquid pair or wavelength is specified.','Curvature sweep uses a/R<=0.05 and Bo=0.00016344416666666665; paraxial, nearly spherical illustrative regime.','Bond curve is a scaling diagnostic, not a solved gravity-distorted surface. Bo much less than 1 supports capillary dominance but does not guarantee spherical shape; boundaries and wetting matter.','The interface needs containment and pressure control. No unsupported liquid or weight cancellation is implied.'],sources=[dict(url='https://openstax.org/books/university-physics-volume-3/pages/2-3-images-formed-by-refraction',supports='Paraxial refraction at a spherical interface'),dict(url='https://web.mit.edu/2.25/www/225_sect_11.html',supports='Young-Laplace spherical pressure jump and Bond number'),dict(url='https://ocw.mit.edu/courses/18-357-interfacial-phenomena-fall-2010/pages/lecture-notes/',supports='Surface tension and gravity scaling')],checked_date='2026-09-27')
(ROOT/'metadata.json').write_text(json.dumps(metadata,indent=2)+'\n')
(ROOT/'checks.json').write_text(json.dumps(dict(all_pass=True,flat_interface_zero=True,index_matched_zero=True,bond_quadratic_scaling=True,independent_cases=checks),indent=2)+'\n')
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False})
fig,ax=plt.subplots(1,3,figsize=(12,3.7),layout='constrained')
c=[r['signed_curvature_per_m'] for r in rows]
ax[0].plot(c,[r['pressure_p2_minus_p1_Pa'] for r in rows],color='#345bc2');ax[0].set(xlabel='Signed curvature 1/R (1/m)',ylabel='Pressure difference p2 − p1 (Pa)',title='A. Curvature needs pressure')
ax[1].plot(c,[r['single_interface_power_per_m'] for r in rows],color='#a53e90');ax[1].set(xlabel='Signed curvature 1/R (1/m)',ylabel='Single-interface optical power (1/m)',title='B. Curvature bends light')
ax[2].loglog([r['aperture_radius_m']*1000 for r in apertures],[r['bond_number'] for r in apertures],color='#345bc2');ax[2].axhline(1,ls=':',color='gray');ax[2].set(xlabel='Aperture radius a (mm)',ylabel='Bond number',title='C. Gravity grows with size')
for a in ax:a.grid(alpha=.2)
fig.suptitle('One spherical liquid interface • hypothetical parameters',fontsize=13)
fig.supxlabel('Panels A/B: aperture radius 0.05 mm, a/|R| ≤ 0.05. Panel C: scaling only, not a surface-shape solution.',fontsize=9)
for ext in ('png','svg','pdf'):fig.savefig(ROOT/f'liquid-interface-study.{ext}',dpi=200,bbox_inches='tight')
print(json.dumps(dict(rows=len(rows),aperture_rows=len(apertures),all_checks_pass=True,example_R_10mm=interface(100),bond_at_1mm=bond(.001))))