ms-09: آیا منحرف‌کنندهٔ پرتو پهنای باند کاربردی دارد؟

from pathlib import Path
import sys, json, hashlib, platform
from dataclasses import asdict, replace
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
ROOT = Path.cwd()
if not (ROOT/'simulations').is_dir(): ROOT = ROOT.parent
sys.path.insert(0, str(ROOT/'simulations'))
from design_studies import (Design, order_metrics, spectral_sweep, bandwidth,
                           local_prediction, optimize_widths, grcwa_check)
from unit_cell_te import Cell, solve_cell, sample_ey
plt.rcParams.update({'figure.dpi': 125, 'font.size': 10, 'axes.spines.top': False,
                     'axes.spines.right': False})
OUT = ROOT/'data/design-studies'; OUT.mkdir(parents=True, exist_ok=True)
def metadata(design):
    import scipy
    return dict(parameters=asdict(design), units='um, degrees; incident-normalized power',
      conventions='exp(-i omega t), TE Ey, r at z=0, t at z=height',
      sources={name:hashlib.sha256((ROOT/'simulations'/name).read_text(encoding='utf-8').encode()).hexdigest()
               for name in ['design_studies.py','unit_cell_te.py','metasurface_synthesis.py']},
      environment={'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__})
def geometry(ax, design, title):
    L=len(design.widths_um)*design.pitch_um; h=design.height_um
    ax.add_patch(Rectangle((0,h),L,.32,fc='#c3ddeb'))
    for j,w in enumerate(design.widths_um):
        ax.add_patch(Rectangle(((j+.5)*design.pitch_um-w/2,0),w,h,fc='#167b9b',ec='white'))
        ax.text((j+.5)*design.pitch_um,-.07,str(j+1),ha='center',fontsize=9)
    ax.axvline(0,c='#9866aa',ls='--');ax.axvline(L,c='#9866aa',ls='--')
    ax.set(xlim=(-.1,L+.1),ylim=(h+.32,-.4),xlabel='x (um)',ylabel='z (um)',title=title)
    ax.text(.05,.04,'TE: E along y; k along +z',transform=ax.transAxes,fontsize=8)
    ax.text(.05,.13,'Repeated in x; invariant in y',transform=ax.transAxes,fontsize=8)

# Physical inputs, performance requirements, and numerical controls are separate.
design = Design()
anchor_um = design.wavelength_um
wavelengths = np.linspace(1.36, 1.74, 161)  # Includes 1.55; avoids exact grazing at 1.4625.
efficiency_threshold = .80
angle_tolerance_deg = .5
M_compare, M_reference = 96, 144
coarse_harmonics = spectral_sweep(design, wavelengths, M_compare)
rows = spectral_sweep(design, wavelengths, M_reference)
max_power_change = max(abs(a['T']-b['T']) for a,b in zip(rows,coarse_harmonics))
assert max_power_change < 1e-3
assert max(r['balance_error'] for r in rows) < 1e-9
bands = []
for stride in [4,2,1]:
    subset=rows[::stride]
    bands.append(dict(samples=len(subset),efficiency=bandwidth(subset,anchor_um,efficiency_threshold),
        usable=bandwidth(subset,anchor_um,efficiency_threshold,angle_tolerance_deg)))
print('Maximum target-power change, M=96 to 144:',max_power_change)
print(json.dumps(bands,indent=2))

lam=np.array([r['wavelength_um'] for r in rows]); power=np.array([r['T'] for r in rows])
angles=np.array([r['angle_deg'] for r in rows]); i0=np.argmin(abs(lam-anchor_um))
fig,ax=plt.subplots(1,3,figsize=(13,3.8),layout='constrained')
geometry(ax[0],design,'Fixed six-ridge geometry')
ax[1].plot(lam,power,label='Target T(+1)',c='#127f93')
ax[1].plot(lam,[r['total_T'] for r in rows],label='Total transmitted power',c='#888',ls='--')
ax[1].axhline(efficiency_threshold,c='#b35540',ls=':',label='Required efficiency')
for band,color,label in [(bands[-1]['efficiency'],'#d4e7f3','Efficiency band'),(bands[-1]['usable'],'#a5d8bd','Useful band')]:
    if band: ax[1].axvspan(band['lower_um'],band['upper_um'],color=color,alpha=.6,label=label)
ax[1].set(xlabel='Vacuum wavelength (um)',ylabel='Incident-power fraction',ylim=(0,1.05),title='Efficiency versus usable bandwidth')
ax[1].legend(fontsize=7)
ax[2].plot(lam,angles,c='#7450a2');ax[2].axhspan(angles[i0]-.5,angles[i0]+.5,color='#a5d8bd',alpha=.5)
ax[2].axvline(1.4625,c='#888',ls=':',label='m = +/-4 substrate cutoff')
ax[2].set(xlabel='Vacuum wavelength (um)',ylabel='Target beam angle (deg)',title='The beam direction also disperses')
ax[2].legend(fontsize=7)
plt.show()

def modal_delay(selected):
    ls=np.array([r['wavelength_um'] for r in selected])[::-1]
    ts=np.array([complex(r['t_real'],r['t_imag']) for r in selected])[::-1]
    omega=2*np.pi*.299792458/ls  # rad/fs; c in um/fs
    phase=np.unwrap(np.angle(ts))
    delay=np.gradient(phase,omega,edge_order=2)
    valid=np.array([r['T'] for r in selected])[::-1]>.05
    # Also mask derivative stencils touching a weak sample.
    valid[1:-1] &= valid[:-2].copy() & valid[2:].copy()
    valid[[0,-1]]=False
    return ls,np.where(valid,delay,np.nan)
fig,ax=plt.subplots(1,2,figsize=(10,3.4),layout='constrained')
for stride in [2,1]:
    ls,delay=modal_delay(rows[::stride]); ax[0].plot(ls,delay,label=f'{len(ls)} wavelengths')
ax[0].set(xlabel='Vacuum wavelength (um)',ylabel='Modal group delay (fs)',title='Derivative: inspect sampling sensitivity');ax[0].legend()
ax[1].semilogy(lam,np.maximum(abs(power-np.array([r['T'] for r in coarse_harmonics])),1e-16),c='#127f93')
ax[1].set(xlabel='Vacuum wavelength (um)',ylabel='Absolute target-power change',title='Independent harmonic refinement')
plt.show()
report=dict(metadata=metadata(design),half_orders=[M_compare,M_reference],rows=rows,
    comparison=coarse_harmonics,bandwidths=bands,efficiency_threshold=efficiency_threshold,
    angle_tolerance_deg=angle_tolerance_deg,max_power_change=max_power_change)
(OUT/'dispersion.json').write_text(json.dumps(report,indent=2)+'\n')