دفترچهٔ فصل ۵: سنتز و بررسی منحرف‌کنندهٔ باریکه

from pathlib import Path
from dataclasses import replace,asdict
import sys,json,csv,hashlib,platform,html
import numpy as np
import scipy,matplotlib
import matplotlib.pyplot as plt
from IPython.display import display,HTML
ROOT=next((p for p in [Path.cwd(),*list(Path.cwd().parents)[:3]] if (p/'simulations/metasurface_synthesis.py').is_file()),None)
if ROOT is None:raise FileNotFoundError('Extract the Chapter 5 source bundle, preserving simulations/, data/, and notebooks/.')
sys.path.insert(0,str(ROOT/'simulations'))
from unit_cell_te import Cell,solve_cell,sample_ey,phase_error_deg
from metasurface_synthesis import solve_supercell,staircase_coefficients,select_widths,angular_spectrum_1d,rotated_jones,line_fwhm,SOLVER_VERSION
DATA=ROOT/'data/chapter-05';DATA.mkdir(parents=True,exist_ok=True)
def table(headers,rows):
    def escape(x):return html.escape(f'{x:.7g}' if isinstance(x,(float,np.floating)) else str(x))
    display(HTML('<table><thead><tr>'+''.join('<th>'+escape(x)+'</th>' for x in headers)+'</tr></thead><tbody>'+''.join('<tr>'+''.join('<td>'+escape(x)+'</td>' for x in row)+'</tr>' for row in rows)+'</tbody></table>'))
plt.rcParams.update({'figure.dpi':125,'font.size':10,'axes.spines.top':False,'axes.spines.right':False})
print('Python',platform.python_version(),'NumPy',np.__version__,'SciPy',scipy.__version__,'Matplotlib',matplotlib.__version__)

# Editable physical and numerical parameters, in micrometers and degrees.
N_CELLS=6;PITCH_UM=.65;WAVELENGTH_UM=1.55;N_OUTPUT=1.5
REFERENCE=Cell(period_um=PITCH_UM,wavelength_um=WAVELENGTH_UM,n_substrate=N_OUTPUT)
REFERENCE.validate()
HEIGHT_CANDIDATES=[.6,.9,1.2,1.5,1.8]
SCREEN_WIDTHS=np.linspace(.08,.57,50)
EXTENDED_HEIGHT_UM=1.8;WIDTHS=np.linspace(.08,.57,161)
M_LIBRARY,M_LIBRARY_CHECK=24,32
M_FINAL,M_CHECK=144,192
MINIMUM_T=.8;PHASE_ACCEPT_DEG=5.
PERIOD_UM=N_CELLS*PITCH_UM
centers=(np.arange(N_CELLS)+.5)*PITCH_UM
target_phases=2*np.pi*centers/PERIOD_UM
target_sine=(REFERENCE.n_incident*np.sin(np.deg2rad(REFERENCE.angle_deg))+REFERENCE.wavelength_um/PERIOD_UM)/REFERENCE.n_substrate
if abs(target_sine)>=1:raise ValueError('The desired +1 transmitted order must be open; revise period, wavelength, or incidence.')
target_angle=float(np.rad2deg(np.arcsin(target_sine)))
print('Supercell period',PERIOD_UM,'um; desired transmitted order +1 at',target_angle,'degrees.')

def require_matching_library(parameters, reference, phase_references, label):
    """Reject reuse outside the saved TE library's experiment and face references."""
    expected=asdict(reference)
    keys=('period_um','wavelength_um','n_ridge','n_gap','n_incident','n_substrate','angle_deg')
    mismatches=[]
    for key in keys:
        actual=parameters.get(key) if isinstance(parameters,dict) else None
        if actual is None or not np.isfinite(actual) or not np.isclose(actual,expected[key],rtol=1e-12,atol=1e-12):
            mismatches.append(f'{key}: saved {actual!r}, requested {expected[key]!r}')
    if phase_references.get('r')!='z=0' or phase_references.get('t')!='z=height':
        mismatches.append('reflection/transmission must be referenced to z=0/z=height')
    if mismatches:
        raise ValueError(f"{label} is incompatible with this experiment: "+'; '.join(mismatches)+
                         '. Generate a matching library before changing its pitch, wavelength, materials, or incidence.')

def solve_sequence(widths, height_um, half_order):
    cell=replace(REFERENCE,height_um=height_um)
    return solve_supercell(widths,pitch_um=cell.period_um,height_um=cell.height_um,
                           wavelength_um=cell.wavelength_um,n_ridge=cell.n_ridge,n_gap=cell.n_gap,
                           n_incident=cell.n_incident,n_substrate=cell.n_substrate,
                           angle_deg=cell.angle_deg,half_order=half_order)

original_data=json.loads((ROOT/'data/chapter-03/response-library.json').read_text(encoding='utf-8'))
original_metadata=original_data['metadata']
require_matching_library(original_metadata.get('parameters',{}),REFERENCE,
                         original_metadata.get('phase_references',{}),'Saved Chapter 3 library')
if (original_metadata.get('time_convention')!='exp(-i omega t)' or
        original_metadata.get('spatial_convention')!='exp(i k.r)'):
    raise ValueError('Saved Chapter 3 library uses incompatible phasor conventions.')
ORIGINAL_HEIGHT_UM=float(original_metadata['parameters']['height_um'])
original=original_data['rows']

screen=[]
for height in HEIGHT_CANDIDATES:
    values=[solve_cell(replace(REFERENCE,width_um=float(w),height_um=height),16) for w in SCREEN_WIDTHS]
    t=np.array([v['t'][16] for v in values]);power=np.array([v['T_total'] for v in values])
    screen.append(dict(height_um=height,phase_span_deg=float(np.rad2deg(np.ptp(np.unwrap(np.angle(t))))),minimum_T=float(power.min()),fraction_above_threshold=float(np.mean(power>=MINIMUM_T))))
table(['Height (um)','Sampled phase span (deg)','Minimum T','Fraction with T >= 0.8'],[[r[k] for k in ['height_um','phase_span_deg','minimum_T','fraction_above_threshold']] for r in screen])
extended=[]
for width in WIDTHS:
    c=replace(REFERENCE,height_um=EXTENDED_HEIGHT_UM,width_um=float(width))
    a,b=solve_cell(c,M_LIBRARY),solve_cell(c,M_LIBRARY_CHECK)
    t,r=a['t'][M_LIBRARY],a['r'][M_LIBRARY]
    dt=float(abs(t-b['t'][M_LIBRARY_CHECK]));dp=float(max(abs(a['T_total']-b['T_total']),abs(a['R_total']-b['R_total'])))
    dphi=float(abs(phase_error_deg(t,b['t'][M_LIBRARY_CHECK])))
    extended.append(dict(width_um=float(width),t_real=float(t.real),t_imag=float(t.imag),r_real=float(r.real),r_imag=float(r.imag),T=a['T_total'],R=a['R_total'],phase_deg=float(np.rad2deg(np.angle(t))),delta_t=dt,delta_power=dp,delta_phase_deg=dphi,balance_error=a['balance_error'],accepted=bool(dt<5e-4 and dp<5e-4 and dphi<.05 and a['balance_error']<1e-8)))
assert all(r['accepted'] for r in extended)
print('Extended library:',len(extended),'accepted rows; worst complex refinement',max(r['delta_t'] for r in extended))
fig,axes=plt.subplots(1,2,figsize=(10,3.8),layout='constrained')
for rows,label in [(original,f'h={ORIGINAL_HEIGHT_UM:.2f} um'),(extended,f'h={EXTENDED_HEIGHT_UM:.2f} um')]:
    width=[r['width_um'] for r in rows];t=np.array([complex(r['t_real'],r['t_imag']) for r in rows])
    axes[0].plot(width,np.rad2deg(np.unwrap(np.angle(t))),label=label)
    axes[1].plot(width,[r['T'] for r in rows],label=label)
axes[0].set(xlabel='Width (um)',ylabel='Unwrapped phase (deg)',title='Same width domain, different common height')
axes[1].axhline(MINIMUM_T,color='gray',ls='--');axes[1].set(xlabel='Width (um)',ylabel='T',ylim=(0,1.05),title='Transmission constrains usable phases')
for ax in axes:ax.legend();ax.grid(alpha=.2)
plt.show()

def choose(rows):
    trials=[]
    for offset in np.linspace(-np.pi,np.pi,361,endpoint=False):
        selected,error=select_widths(rows,target_phases,MINIMUM_T,float(offset))
        trials.append((float(np.mean(error**2)),float(offset),selected,error))
    score,offset,selected,error=min(trials,key=lambda x:x[0])
    return dict(offset_rad=offset,selected=selected,errors_deg=np.rad2deg(error).tolist(),rms_phase_error_deg=float(np.rad2deg(np.sqrt(score))))
limited,expanded=choose(original),choose(extended)
table(['Library height (um)','Global offset (deg)','RMS phase error (deg)','Maximum error (deg)','All targets within 5 deg'],[[height,np.rad2deg(v['offset_rad']),v['rms_phase_error_deg'],max(abs(np.array(v['errors_deg']))),bool(max(abs(np.array(v['errors_deg'])))<=PHASE_ACCEPT_DEG)] for height,v in [(ORIGINAL_HEIGHT_UM,limited),(EXTENDED_HEIGHT_UM,expanded)]])
table(['Position','Width (um)','Local T','Phase error (deg)'],[[j,r['width_um'],r['T'],expanded['errors_deg'][j]] for j,r in enumerate(expanded['selected'])])
assert max(abs(np.array(expanded['errors_deg'])))<PHASE_ACCEPT_DEG

case_specs=[('limited',ORIGINAL_HEIGHT_UM,limited['selected']),('expanded',EXTENDED_HEIGHT_UM,expanded['selected']),
            ('uniform',EXTENDED_HEIGHT_UM,[min(extended,key=lambda r:abs(r['width_um']-REFERENCE.period_um/2))]*N_CELLS),
            ('reversed',EXTENDED_HEIGHT_UM,expanded['selected'][::-1])]
cases=[];raw={}
for name,height,selected in case_specs:
    widths=[r['width_um'] for r in selected]
    result=solve_sequence(widths,height,M_FINAL)
    check=solve_sequence(widths,height,M_CHECK)
    ct=staircase_coefficients([complex(r['t_real'],r['t_imag']) for r in selected],result['orders'])
    cr=staircase_coefficients([complex(r['r_real'],r['r_imag']) for r in selected],result['orders'])
    lpaT=result['qt'].real/result['qi'][M_FINAL].real*abs(ct)**2
    lpaR=result['qi'].real/result['qi'][M_FINAL].real*abs(cr)**2
    open_either=(result['qi'].real>0)|(result['qt'].real>0)
    rows=[]
    for index in np.where(open_either)[0]:
        m=int(result['orders'][index]);j=m+M_CHECK
        rows.append(dict(m=m,reflection_open=bool(result['qi'][index].real>0),transmission_open=bool(result['qt'][index].real>0),R=float(result['R'][index]),T=float(result['T'][index]),local_R=float(lpaR[index]),local_T=float(lpaT[index]),t_real=float(result['t'][index].real),t_imag=float(result['t'][index].imag),delta_T=float(abs(result['T'][index]-check['T'][j])),delta_R=float(abs(result['R'][index]-check['R'][j])),delta_complex_t=float(abs(result['t'][index]-check['t'][j]))))
    assert max(r['delta_T'] for r in rows)<1e-4 and max(r['delta_R'] for r in rows)<1e-4
    assert max(r['delta_complex_t'] for r in rows if r['transmission_open'])<1e-4
    assert result['balance_error']<1e-8
    target=next(r for r in rows if r['m']==1)
    record=dict(name=name,height_um=height,widths_um=widths,local_t=[[r['t_real'],r['t_imag']] for r in selected],orders=rows,R_total=result['R_total'],T_total=result['T_total'],target_T=target['T'],conditional_target_fraction=target['T']/result['T_total'],local_target_T=target['local_T'],local_balance=float(lpaT.sum()+lpaR.sum()),balance_error=result['balance_error'])
    cases.append(record);raw[name]=result
table(['Design','T(+1), incident-normalized','T(+1)/Ttotal','Local-mask T(+1)','Rtotal','Ttotal','Local R+T'],[[v[k] for k in ['name','target_T','conditional_target_fraction','local_target_T','R_total','T_total','local_balance']] for v in cases])
table(['Order','R','T','Local R estimate','Local T estimate','Delta T: 144 to 192'],[[r[k] for k in ['m','R','T','local_R','local_T','delta_T']] for r in cases[1]['orders']])
assert cases[2]['target_T']<1e-12
assert abs(next(r['T'] for r in cases[3]['orders'] if r['m']==-1)-cases[1]['target_T'])<1e-10
fig,axes=plt.subplots(1,2,figsize=(11,4),layout='constrained')
orders=[r['m'] for r in cases[1]['orders']]
axes[0].bar(np.array(orders)-.18,[r['local_T'] for r in cases[1]['orders']],width=.36,label='Local staircase estimate')
axes[0].bar(np.array(orders)+.18,[r['T'] for r in cases[1]['orders']],width=.36,label='Full-wave supercell')
axes[0].set(xlabel='Transmitted order',ylabel='Incident-normalized estimate / power',title='Expanded library: geometry interactions matter');axes[0].legend()
for case in cases:axes[1].plot(orders,[r['T'] for r in case['orders']],'o-',label=case['name'])
axes[1].set(xlabel='Transmitted order',ylabel='Full-wave transmitted power',title='Limited, expanded, identical, and reversed sequences');axes[1].legend()
plt.show()

convergence=[]
for M in sorted(set([12,24,36,48,72,96,M_FINAL])):
    result=solve_sequence(cases[1]['widths_um'],cases[1]['height_um'],M)
    convergence.append(dict(M=M,harmonics=2*M+1,target_T=float(result['T'][M+1]),delta_target_vs_check=0.,balance_error=result['balance_error']))
# Use a fresh explicit high-order target for signed-difference-independent errors.
high=solve_sequence(cases[1]['widths_um'],cases[1]['height_um'],M_CHECK)
for row in convergence:row['delta_target_vs_check']=abs(row['target_T']-float(high['T'][M_CHECK+1]))
table(['M','Harmonics','T(+1)',f'Difference vs M={M_CHECK}','Power residual'],[[r[k] for k in ['M','harmonics','target_T','delta_target_vs_check','balance_error']] for r in convergence])

result=raw['expanded']
x=np.linspace(0,2*PERIOD_UM,361);z=np.linspace(-.35,result['cell'].height_um+5,301)
field=sample_ey(result['cell'],result,x,z)
fig,ax=plt.subplots(figsize=(9,5.2),layout='constrained')
image=ax.pcolormesh(x,z,field.real,shading='auto',cmap='RdBu_r');fig.colorbar(image,ax=ax,label='Re(Ey/Einc)')
for repeat in [0,1]:
    for center,width in zip(centers+repeat*PERIOD_UM,cases[1]['widths_um']):ax.add_patch(plt.Rectangle((center-width/2,0),width,result['cell'].height_um,fill=False,color='black',lw=.7))
ax.axhline(result['cell'].height_um,color='black',lw=.7);ax.invert_yaxis();ax.set(xlabel='x (um)',ylabel='z (um)',title=f'Full-wave repeated supercell; +1 order at {target_angle:.2f}° in the substrate')
plt.show()

metadata=dict(topic_id='ms-05-beam-steering',solver=SOLVER_VERSION,units='micrometers/degrees',time_convention='exp(-i omega t)',spatial_convention='exp(i k.r)',phase_references={'r':'z=0','t':'z=common height'},pitch_um=PITCH_UM,period_um=PERIOD_UM,wavelength_um=WAVELENGTH_UM,n_incident=REFERENCE.n_incident,n_substrate=REFERENCE.n_substrate,n_ridge=REFERENCE.n_ridge,n_gap=REFERENCE.n_gap,angle_deg=REFERENCE.angle_deg,library_parameters=asdict(replace(REFERENCE,height_um=EXTENDED_HEIGHT_UM)),library_phase_references={'r':'z=0','t':'z=height'},target_order=1,target_angle_deg=target_angle,half_order=M_FINAL,check_half_order=M_CHECK,library_half_order=M_LIBRARY,library_check_half_order=M_LIBRARY_CHECK,minimum_T=MINIMUM_T,selection='361 fixed global offsets; minimize circular phase RMS using unit-cell data only',source_sha256={name:hashlib.sha256((ROOT/'simulations'/name).read_bytes()).hexdigest() for name in ['unit_cell_te.py','metasurface_synthesis.py']},environment=dict(python=platform.python_version(),numpy=np.__version__,scipy=scipy.__version__),limits=['1D TE, real constant dielectric indices','Infinite periodic supercell','Common ridge height within each design','No FEM/FDTD or experiment','Focusing and polarization are separate scalar/Jones models'])
output=dict(metadata=metadata,height_screen=screen,extended_library=extended,selection={name:{k:v for k,v in choice.items() if k!='selected'} for name,choice in [('limited',limited),('expanded',expanded)]},cases=cases,convergence=convergence)
(DATA/'beam-steering.json').write_text(json.dumps(output,indent=2)+'\n',encoding='utf-8')
with (DATA/'extended-library.csv').open('w',newline='',encoding='utf-8') as handle:
    writer=csv.DictWriter(handle,fieldnames=list(extended[0]));writer.writeheader();writer.writerows(extended)
print('Saved beam-steering.json and extended-library.csv')
print('Expanded target-order power:',cases[1]['target_T'],'conditional transmitted fraction:',cases[1]['conditional_target_fraction'])