%% coupled_microstrip_xtalk_demo.m
% 純Matlabスクリプト：結合マイクロストリップのNEXT/FEXTを計算・表示
% - 幾何・材質・長さ・周波数レンジはここで記述（引数なし）
% - 閉形式の偶奇モード(Z0e/Z0o, eps_e/eps_o) + 簡易損失 -> 2ポート(偶/奇) -> 合成

clear; clc;

%% ── パラメータ（SI単位） ─────────────────────────────────────────────
w    = 0.20e-3;    % トレース幅 [m]
s    = 0.20e-3;    % トレース間 隙間（エッジtoエッジ）[m]
h    = 0.10e-3;    % 基板厚 [m]
t    = 0.035e-3;   % 銅厚 [m]
er   = 4.5;        % 比誘電率
tand = 0.04;       % 誘電正接
sigma= 5.8e7;      % 銅導電率 [S/m]
rough= 1.0e-6;     % 表面粗さの代表長さ [m]（0で無視）
len  = 10e-3;      % 結合区間の長さ [m]
z0_term = 50;      % 端子インピーダンス（全ポート同一想定）
f    = linspace(10e6, 15e9, 4001);  % 周波数ベクトル [Hz]

%% ── 物理定数など ────────────────────────────────────────────────────
c0  = 299792458;
mu0 = 4*pi*1e-7;
f   = f(:).';                % 行ベクトル
w_h = w/h; s_h = s/h; t_h = t/h;

%% ── 単線のベース（Hammerstad–Jensen）で有効幅と有効比誘電率の種を作成 ──
% 厚銅の有効幅補正（簡易）
if t>0
    du = (t/pi) * (1 + 1/er) / 2;    % ざっくり補正
else
    du = 0;
end
we = w + du;                         % 有効幅（ラフ）
u  = we/h;

% Hammerstad–Jensen の eeff（単線）
a = 1 + (1/49)*log((u^4 + (u/52)^2)/(u^4 + 0.432)) + (1/18.7)*log(1 + (u/18.1)^3);
b = 0.564*((er-0.9)/(er+3))^(0.053);
eeff_single = (er+1)/2 + ((er-1)/2) * (1 + 10/u)^(-a*b);

% 単線Z0の粗い種
Z0s = 60/sqrt(eeff_single) * log(8/u + 0.25*u);

%% ── 結合線の偶奇モード（Wadell系近似） ───────────────────────────────
k = w_h/(w_h + 2*s_h);
k = min(max(k, 1e-3), 0.999);
% 有効比誘電率（偶/奇）近似
ee = er - (er - eeff_single).*exp(-(1.2*s_h).^0.8) .* (1 - 0.5*k);
eo = er - (er - eeff_single).*exp(-(0.6*s_h).^0.9) .* (1 - 1.2*k.^0.7);
ee = min(max(ee, 1.0001), er);
eo = min(max(eo, 1.0001), er);

% モード特性インピーダンス（偶/奇）近似
ce = 1 + 0.5*exp(-1.2*s_h) + 0.08/(1+s_h^1.2);
co = 1 - 0.5*exp(-0.8*s_h) - 0.05/(1+s_h^1.3);
Z0e = max( Z0s / sqrt(ee/eeff_single) / ce, 5 );
Z0o = max( Z0s / sqrt(eo/eeff_single) / co, 5 );

%% ── 損失モデル（周波数依存）────────────────────────────────────────
% 表面抵抗
Rs = sqrt(pi*mu0*f/sigma);           % [ohm/square]
we_eff = max(we, 1e-12);
% 導体損（ラフ、モード別にZ0でスケール）
alpha_ce = (Rs ./ (2*Z0e*we_eff));   % [Np/m]
alpha_co = (Rs ./ (2*Z0o*we_eff));   % [Np/m]

% 粗さの簡易係数（任意の軽い増幅）
if rough>0
    Kr = 1 + (2/pi)*atan( (2*pi*rough)./sqrt(2).*sqrt(f/(mu0*sigma)) );
    alpha_ce = alpha_ce .* Kr;
    alpha_co = alpha_co .* Kr;
end

% 誘電損： alpha_d ≈ (pi*f/c)*sqrt(eps_eff)*tanδ
alpha_de = (pi .* f ./ c0) .* sqrt(ee) .* tand;
alpha_do = (pi .* f ./ c0) .* sqrt(eo) .* tand;

alpha_e = alpha_ce + alpha_de;       % [Np/m]
alpha_o = alpha_co + alpha_do;       % [Np/m]

% 位相定数
beta_e  = 2*pi.*f .* sqrt(ee) / c0;  % [rad/m]
beta_o  = 2*pi.*f .* sqrt(eo) / c0;  % [rad/m]
gamma_e = alpha_e + 1j*beta_e;
gamma_o = alpha_o + 1j*beta_o;

%% ── 偶/奇モード 2ポート（長さlen）→ S（各モード）─────────────────────
% ABCD（モードごと、ベクトル化）
Ae = cosh(gamma_e*len);
Be = Z0e .* sinh(gamma_e*len);
Ce = (1./Z0e) .* sinh(gamma_e*len);
De = Ae;

Ao = cosh(gamma_o*len);
Bo = Z0o .* sinh(gamma_o*len);
Co = (1./Z0o) .* sinh(gamma_o*len);
Do = Ao;

Zs = z0_term; Zl = z0_term;

% ABCD→S（モード）: S11 = (A+B/Zl - C*Zs - D*Zs/Zl)/den, S21 = 2/den, den= A+B/Zl + C*Zs + D*Zs/Zl
dene = (Ae + Be./Zl + Ce.*Zs + De.*(Zs./Zl));
Te11 = (Ae + Be./Zl - Ce.*Zs - De.*(Zs./Zl)) ./ dene;
Te21 = 2 ./ dene;

deno = (Ao + Bo./Zl + Co.*Zs + Do.*(Zs./Zl));
To11 = (Ao + Bo./Zl - Co.*Zs - Do.*(Zs./Zl)) ./ deno;
To21 = 2 ./ deno;

%% ── 単一励振（Port1励振、他マッチ）での合成：NEXT/FEXT ────────────────
% 標準関係（対称二線路）：
%  S11 = 0.5*(Te11+To11), S31 = 0.5*(Te11-To11)  → NEXT = |S31|
%  S21 = 0.5*(Te21+To21), S41 = 0.5*(Te21-To21)  → FEXT = |S41|
S31 = 0.5*(Te11 - To11);   % NEXT
S41 = 0.5*(Te21 - To21);   % FEXT

NEXT_dB = 20*log10(abs(S31)+eps);
FEXT_dB = 20*log10(abs(S41)+eps);

%% ── プロット ──────────────────────────────────────────────────────────
figure; plot(f/1e9, NEXT_dB, 'LineWidth',1.5); grid on;
xlabel('Frequency [GHz]'); ylabel('NEXT [dB]'); title('Near-End Crosstalk');

figure; plot(f/1e9, FEXT_dB, 'LineWidth',1.5); grid on;
xlabel('Frequency [GHz]'); ylabel('FEXT [dB]'); title('Far-End Crosstalk');

%% ── 参考出力（偶奇モードパラメータの一部）──────────────────────────────
fprintf('Z0e = %.2f ohm, Z0o = %.2f ohm\n', Z0e, Z0o);
fprintf('eeff_even = %.3f, eeff_odd = %.3f\n', ee, eo);
