Untitled
unknown
plain_text
a month ago
14 kB
6
Indexable
"""
=============================================================================
ELC4039 — Optical Fibers | Assignment 5
Dispersion Curves of a Step-Index Circular Optical Fiber
Dr. H. Kirolous, Spring 2026
=============================================================================
METHODOLOGY — "Experimentation" as required by the assignment
─────────────────────────────────────────────────────────────
The assignment states parameters n1, n2, and a must be selected "carefully
(via experimentation)" so that different modes can be clearly distinguished.
Three parameters control clarity:
1. n1, n2 → index contrast Δ = (n1−n2)/n1
• Too small Δ: modes overlap, cutoffs cluster together
• Too large Δ: fewer modes visible in V = 0→7 range
2. a → physical core radius (µm)
• Does NOT change the b(V) curves (V is normalised)
• Only shifts the secondary optical-frequency axis
• Determines where telecom windows (1310/1550 nm) fall
EXPERIMENTS
───────────
Experiment 1: vary index contrast (5 combinations of n1, n2)
Experiment 2: vary core radius (4 values of a, n1/n2 fixed)
FINAL CHOICE: n1 = 1.50, n2 = 1.36, a = 3.4 µm
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from scipy.special import jv, kv, jvp, kvp
from scipy.optimize import brentq
import os
os.makedirs("outputs/experiments", exist_ok=True)
c = 3e8 # speed of light [m/s]
# ─── Core solver ─────────────────────────────────────────────────────────────
def char_eq(b, V, m, n1, n2):
"""
Eigenvalue equation (Lecture 5, Slide 9).
Returns LHS − RHS; root ↔ a guided mode.
"""
if b <= 1e-9 or b >= 1.0 - 1e-9:
return np.nan
w = V * np.sqrt(b)
u = V * np.sqrt(1.0 - b)
Jm = jv(m, u); Km = kv(m, w)
if abs(Jm) < 1e-30 or abs(Km) < 1e-30:
return np.nan
FJ = jvp(m, u) / (u * Jm)
FK = kvp(m, w) / (w * Km)
beta_n = np.sqrt(n2**2 + b * (n1**2 - n2**2))
lhs = (FJ + FK) * (n1**2 * FJ + n2**2 * FK)
rhs = (beta_n * m)**2 * (1.0/u**2 + 1.0/w**2)**2
return lhs - rhs
def compute_branches(V_array, m, n1, n2, n_b=1500):
"""Find and stitch dispersion branches for azimuthal order m."""
b_grid = np.linspace(1e-5, 1.0 - 1e-5, n_b)
branches = []
for V in V_array:
if V < 0.05:
continue
vals = np.array([char_eq(b, V, m, n1, n2) for b in b_grid])
roots_V = []
for i in range(len(vals)-1):
v0, v1 = vals[i], vals[i+1]
if np.isnan(v0) or np.isnan(v1):
continue
if v0 * v1 < 0:
try:
r = brentq(char_eq, b_grid[i], b_grid[i+1],
args=(V, m, n1, n2), xtol=1e-9, maxiter=200)
roots_V.append(r)
except Exception:
pass
roots_V.sort()
while len(branches) < len(roots_V):
branches.append([])
used = [False]*len(roots_V)
for br in branches:
if not br:
continue
_, b_prev = br[-1]
best_d, best_j = 1e9, -1
for j, rb in enumerate(roots_V):
if not used[j] and abs(rb-b_prev) < best_d:
best_d = abs(rb-b_prev); best_j = j
if best_j >= 0 and best_d < 0.12:
br.append((V, roots_V[best_j])); used[best_j] = True
for j, rb in enumerate(roots_V):
if not used[j]:
branches.append([(V, rb)])
result = []
for br in branches:
if len(br) < 6:
continue
Vs, bs = zip(*br)
result.append((np.array(Vs), np.array(bs)))
result.sort(key=lambda x: x[0][0])
return result
def mode_label(m, idx):
"""Mode naming per Lecture 5 Slides 11 & 21."""
if m == 0:
n = idx+1
return f"$TE_{{0{n}}}$" if idx%2==0 else f"$TM_{{0{n}}}$"
elif m == 1:
return f"$HE_{{1{idx+1}}}$"
else:
if idx == 0:
return f"$HE_{{{m}1}}$"
return f"$EH_{{{m-1}{idx}}}$"
COLORS = {0:"steelblue", 1:"darkorange", 2:"seagreen", 3:"crimson", 4:"purple"}
LEGEND_HANDLES = [
Line2D([0],[0], color="steelblue", lw=2, label="$TE_{0n}$/$TM_{0n}$ (m=0)"),
Line2D([0],[0], color="darkorange", lw=2, label="$HE_{1n}$ (m=1)"),
Line2D([0],[0], color="seagreen", lw=2, label="$HE_{2n}$/$EH_{1n}$ (m=2)"),
Line2D([0],[0], color="crimson", lw=2, label="$HE_{3n}$/$EH_{2n}$ (m=3)"),
Line2D([0],[0], color="purple", lw=2, label="$HE_{4n}$/$EH_{3n}$ (m=4)"),
]
def draw_panel(ax, n1, n2, V_array, a_phys=None, n_b=1500):
"""Draw one dispersion panel on ax, return NA."""
NA = np.sqrt(n1**2 - n2**2)
for m in range(5):
branches = compute_branches(V_array, m, n1, n2, n_b=n_b)
for idx, (Vs, bs) in enumerate(branches):
if Vs[0] > 6.8:
continue
ax.plot(Vs, bs, color=COLORS[m], lw=1.5)
x_a = min(Vs[-1], 6.88)
i_a = np.argmin(np.abs(Vs - x_a))
ax.annotate(mode_label(m, idx),
xy=(Vs[i_a], bs[i_a]),
xytext=(3, 0), textcoords="offset points",
fontsize=7, color=COLORS[m], va="center")
ax.axhline(1.0, color="gray", lw=0.7, ls="--")
ax.axhline(0.0, color="gray", lw=0.7, ls="--")
ax.text(6.93, 1.015, "$n_1$", fontsize=7.5, color="gray", ha="right")
ax.text(6.93, 0.012, "$n_2$", fontsize=7.5, color="gray", ha="right")
ax.set_xlim(0, 7); ax.set_ylim(-0.04, 1.08)
ax.set_xlabel("$V$", fontsize=9); ax.set_ylabel("$b$", fontsize=9)
ax.tick_params(labelsize=8); ax.grid(True, linestyle=":", alpha=0.4)
if a_phys is not None:
ax2 = ax.twiny()
ax2.set_xlim(ax.get_xlim())
Vt = np.array([1,2,3,4,5,6,7])
ft = Vt * c / (2*np.pi*a_phys*NA) * 1e-12
ax2.set_xticks(Vt); ax2.set_xticklabels([f"{f:.0f}" for f in ft], fontsize=6)
ax2.set_xlabel("$f$ [THz]", fontsize=7)
return NA
# V arrays: coarse for experiments (fast), fine for final plot
V_coarse = np.linspace(0.05, 7.0, 280)
V_fine = np.linspace(0.05, 7.0, 600)
# =============================================================================
# EXPERIMENT 1 — Vary index contrast (a = 4 µm fixed)
# =============================================================================
print("Experiment 1: index contrast ...")
contrast_cases = [
(1.50, 1.44, "Exp A: n₂=1.44\n"),
(1.50, 1.40, "Exp B: n₂=1.40\n"),
(1.50, 1.38, "Exp C: n₂=1.38\n"),
(1.50, 1.36, "Exp D: n₂=1.36\n"), # CHOSEN
(1.50, 1.34, "Exp E: n₂=1.34\n"),
]
fig, axes = plt.subplots(2, 3, figsize=(17, 10))
axes = axes.flatten()
for i, (n1_, n2_, comment) in enumerate(contrast_cases):
draw_panel(axes[i], n1_, n2_, V_coarse, a_phys=4e-6, n_b=1200)
fc = "honeydew" if "CHOSEN" in comment else "white"
axes[i].set_facecolor(fc)
axes[i].set_title(comment, fontsize=8.5, loc="left",
pad=22 if i < 3 else 4,
color="darkgreen" if "CHOSEN" in comment else "black")
axes[5].axis("off")
axes[5].legend(handles=LEGEND_HANDLES, loc="upper center",
fontsize=10, title="Mode families", title_fontsize=11)
axes[5].text(0.5, 0.08,
"Selection criterion:\n"
"• All five mode families must be\n"
" individually distinguishable\n"
"• Curves span V = 0 → 7 smoothly\n",
transform=axes[5].transAxes, fontsize=10,
va="bottom", ha="center",
bbox=dict(boxstyle="round", fc="lightyellow", ec="goldenrod", lw=1.5))
fig.suptitle(
"Assignment 5 — Experiment 1: Varying Index Contrast Δ = (n₁−n₂)/n₁"
" [n₁ = 1.50 fixed, a = 4 µm fixed]",
fontsize=13, fontweight="bold")
plt.tight_layout(rect=[0, 0, 1, 0.96])
fig.savefig("outputs/experiments/exp1_index_contrast.png",
dpi=120, bbox_inches="tight")
plt.close()
print(" → exp1_index_contrast.png")
# =============================================================================
# EXPERIMENT 2 — Vary core radius a (n1 = 1.50, n2 = 1.44 fixed)
# =============================================================================
print("Running Experiment 2: core radius ...")
# n2=1.44 is used here (not the final n2=1.36) so that the 1310/1550 nm
# telecom window markers fall inside the V=0→7 window.
# The b(V) curve shapes are independent of a — only the f-axis shifts.
n1_exp2, n2_exp2 = 1.50, 1.44
NA_exp2 = np.sqrt(n1_exp2**2 - n2_exp2**2)
radius_cases = [
(3e-6, "a = 3 µm\n"),
(3.4e-6, "a = 3.4 µm\n"), # CHOSEN
(4e-6, "a = 4 µm\n"),
(5e-6, "a = 5 µm\n"),
]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()
for i, (a_, comment) in enumerate(radius_cases):
chosen = ("CHOSEN" in comment)
# FIX 1: removed 'chosen=chosen' — draw_panel does not accept that argument
draw_panel(axes[i], n1_exp2, n2_exp2, V_coarse, a_phys=a_, n_b=1200)
axes[i].set_title(
comment, fontsize=9, loc="left",
pad=22 if i < 2 else 4,
color="darkgreen" if chosen else "black"
)
# Mark 1310 nm and 1550 nm telecom windows as vertical dashed lines
for lam, col in [(1310e-9, "navy"), (1550e-9, "darkred")]:
V_lam = 2 * np.pi * a_ * NA_exp2 / lam
if 0 < V_lam < 7:
axes[i].axvline(V_lam, color=col, lw=1.5, ls="--", alpha=0.85)
axes[i].text(
V_lam + 0.05, 0.85,
f"{int(lam * 1e9)} nm\nV = {V_lam:.2f}",
fontsize=7.5, color=col, va="top",
bbox=dict(fc="white", ec="none", alpha=0.7, pad=1)
)
fig.suptitle(
"Assignment 5 — Experiment 2: Varying Core Radius a\n"
"(n1=1.50, n2=1.44 used here so telecom windows fall inside V=0-7)\n"
" b(V) curves are IDENTICAL for all a — only the secondary f-axis shifts\n"
" Dashed lines: 1310 nm (blue) and 1550 nm (red) telecom windows",
fontsize=11, fontweight="bold"
)
plt.tight_layout(rect=[0, 0, 1, 0.90])
fig.savefig("outputs/experiments/exp2_core_radius.png",
dpi=120, bbox_inches="tight")
plt.close()
print(" → exp2_core_radius.png")
# =============================================================================
# FINAL OPTIMUM PLOT — n1=1.50, n2=1.36, a=3.4 µm
# =============================================================================
print("Final optimum plot ...")
n1, n2, a_phys = 1.50, 1.36, 3.4e-6
NA = np.sqrt(n1**2 - n2**2)
fig, ax = plt.subplots(figsize=(14, 8))
for m in range(5):
branches = compute_branches(V_fine, m, n1, n2, n_b=3000)
for idx, (Vs, bs) in enumerate(branches):
if Vs[0] > 6.8:
continue
ax.plot(Vs, bs, color=COLORS[m], lw=1.9)
x_a = min(Vs[-1], 6.90)
i_a = np.argmin(np.abs(Vs - x_a))
ax.annotate(mode_label(m, idx),
xy=(Vs[i_a], bs[i_a]),
xytext=(3, 1), textcoords="offset points",
fontsize=9, color=COLORS[m], va="center")
# FIX 2: Guard telecom markers — with n2=1.36, a=3.4µm both windows
# land at V > 7 (1550nm→V≈8.7, 1310nm→V≈10.3) so they won't appear,
# but the if-guard prevents any out-of-range drawing errors.
for lam, col, nm in [(1310e-9, "navy", "1310 nm"), (1550e-9, "darkred", "1550 nm")]:
V_lam = 2 * np.pi * a_phys * NA / lam
if 0 < V_lam < 7:
ax.axvline(V_lam, color=col, lw=1.0, ls=":", alpha=0.6)
ax.text(V_lam + 0.04, 1.07, f"{nm}\n(V={V_lam:.2f})",
fontsize=7.5, color=col, va="top")
ax.axhline(1.0, color="gray", lw=0.9, ls="--")
ax.axhline(0.0, color="gray", lw=0.9, ls="--")
ax.text(6.97, 1.015, "$n_1$", fontsize=10, color="gray", ha="right")
ax.text(6.97, 0.013, "$n_2$", fontsize=10, color="gray", ha="right")
ax2 = ax.twiny()
ax2.set_xlim(ax.get_xlim())
Vt = np.array([1, 2, 3, 4, 5, 6, 7])
ft = Vt * c / (2 * np.pi * a_phys * NA) * 1e-12
ax2.set_xticks(Vt)
ax2.set_xticklabels([f"{f:.1f}" for f in ft], fontsize=8)
# FIX 3: Use :.1f so 3.4 µm displays correctly (:.0f would round to "3 µm")
ax2.set_xlabel(f"Optical frequency $f$ [THz] ($a = {a_phys*1e6:.1f}$ µm)",
fontsize=11)
ax.set_xlim(0, 7); ax.set_ylim(-0.04, 1.14)
ax.set_xlabel(
"Normalized frequency $V = (2\\pi a/\\lambda_o)\\,\\mathrm{NA}$",
fontsize=12)
ax.set_ylabel(
r"Normalized propagation constant $b = \dfrac{\beta/k_o - n_2}{n_1 - n_2}$",
fontsize=12)
ax.set_title(
"Dispersion Curves of a Step-Index Circular Optical Fiber\n"
f"$n_1={n1}$, $n_2={n2}$, $\\Delta={(n1-n2)/n1:.2%}$, "
f"NA$={NA:.4f}$, $a={a_phys*1e6:.1f}$ µm "
"— OPTIMUM parameters (selected after experimentation)\n"
"ELC4039 Assignment 5 | Dr. H. Kirolous, Spring 2026",
fontsize=11)
ax.grid(True, linestyle=":", alpha=0.5)
ax.legend(handles=LEGEND_HANDLES, loc="upper left", fontsize=9.5, framealpha=0.9)
plt.tight_layout()
fig.savefig("outputs/fiber_dispersion_FINAL.png",
dpi=150, bbox_inches="tight")
plt.close()
print(" → fiber_dispersion_FINAL.png")Editor is loading...
Leave a Comment