#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
================================================================================
INSTITUTO NACIONAL DE CIÊNCIA E TECNOLOGIA EM TECNOLOGIAS FOTÔNICAS (IN-Foton)
LABORATÓRIO DE TELECOMUNICAÇÕES (LabTel) / DEPARTAMENTO DE ENGENHARIA ELÉTRICA (DEE)
UNIVERSIDADE FEDERAL DO ESPÍRITO SANTO (UFES)
================================================================================
Projeto Educacional: Portal IN-Foton / LabTel (in-foton.org)
Módulo: [MOD-MTL] Análise de Linhas de Transmissão com Multicondutores (Material Bônus)
Coordenação e Autoria: Prof. Marcelo Eduardo Vieira Segatto

Descrição:
Simulação eletromagnética completa e rigorosa da microfita simétrica de três
condutores acoplados (Figura 5.32 da obra canônica de Clayton R. Paul, 
"Analysis of Multiconductor Transmission Lines", 2ª Ed., Wiley, 2008):
1. Extração analítica/numérica das matrizes de capacitância [C], capacitância no
   ar [C0] e indutância [L] pelo Método das Imagens de Clayton Paul;
2. Variação paramétrica dos acoplamentos mútuos em função de s/h (Figuras 5.34 e 5.35);
3. Solução 2D da Equação de Poisson na seção transversal para mapeamento dos três
   modos normais fundamentais (Modo 1, Modo 2 e Modo 3);
4. Análise da dispersão modal (velocidades de fase distintas v1, v2, v3) decorrente
   da interface heterogênea ar/dielétrico;
5. Propagação temporal de pacote de RF em 2.4 GHz e cálculo do comprimento de
   batimento direcional periódico L_beat ~ 59.7 cm.

Licença: Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional (CC BY-SA 4.0)
         e Licença MIT para execução e reutilização do código computacional.
Requisitos: Python 3.8+, numpy, scipy, matplotlib, pillow
Uso: python simular_mtl_3linhas.py
================================================================================
"""

import os
import sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from PIL import Image

# Configuração do diretório de saída
DIR_OUTPUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "outputs")
os.makedirs(DIR_OUTPUT, exist_ok=True)

# Constantes Físicas Universais
EPS0 = 8.8541878128e-12
MU0  = 4.0 * np.pi * 1e-7
SPEED_OF_LIGHT = 1.0 / np.sqrt(MU0 * EPS0)

# Parâmetros Físicos do Substrato FR-4 e Geometria em 2.4 GHz
EPS_R = 4.40       # Permissividade relativa do substrato FR-4
H_SUB = 1.575e-3   # Espessura do dielétrico: 1.575 mm = 62.0 mils
W_LINE = 3.03e-3   # Largura de cada pista: 3.03 mm (Z0 isolada = 50.0 ohms)
S_SEP  = 1.50e-3   # Espaçamento entre bordas: 1.50 mm (equidistante)
F0     = 2.40e9    # Frequência central: 2.40 GHz (banda ISM)

print("="*80)
print("IN-Foton / LabTel / UFES — SIMULAÇÃO MTL DE CLAYTON PAUL (FIGURA 5.32)")
print("Coordenação: Prof. Marcelo Eduardo Vieira Segatto")
print("="*80)
print(f"Substrato: FR-4 (eps_r = {EPS_R}, h = {H_SUB*1e3:.3f} mm)")
print(f"Condutores: w = {W_LINE*1e3:.2f} mm, s = {S_SEP*1e3:.2f} mm, f0 = {F0*1e-9:.2f} GHz")
print(f"Diretório de saída: {DIR_OUTPUT}")
print("="*80)

# -------------------------------------------------------------------------
# 1. MÉTODO DAS IMAGENS DE CLAYTON PAUL (SEÇÃO 5.3.1.2)
# -------------------------------------------------------------------------
def F_kernel(u, a):
    """Núcleo de integração sem singularidades de Clayton Paul (Eq. 5.90b)."""
    if abs(u) < 1e-15:
        return 0.0
    return u * np.log((u**2 + a**2) / (u**2)) + 2.0 * a * np.arctan(u / a)

def calc_mtl_matrices(w=W_LINE, s=S_SEP, h=H_SUB, eps_r=EPS_R, n_lands=3, n_sub=25, n_terms=60):
    w_sub = w / n_sub
    total_sub = n_lands * n_sub
    x_centers = np.zeros(total_sub)
    idx = 0
    for land in range(n_lands):
        x_start = land * (w + s)
        for sub in range(n_sub):
            x_centers[idx] = x_start + (sub + 0.5) * w_sub
            idx += 1

    def compute_P(k_val, eps_eff_val):
        P_mat = np.zeros((total_sub, total_sub))
        n_arr = np.arange(1, n_terms + 1)
        sign_k = (-1.0)**(n_arr + 1) * (k_val**(n_arr - 1))

        for i in range(total_sub):
            xi = x_centers[i]
            for j in range(total_sub):
                xj = x_centers[j]
                d = abs(xi - xj)
                # Termo do ar e imagem no plano de terra
                u1 = d + w_sub
                u2 = d - w_sub
                u0 = d
                P_self = (F_kernel(u1, 0) + F_kernel(u2, 0) - 2*F_kernel(u0, 0)) - \
                         (F_kernel(u1, 2*h) + F_kernel(u2, 2*h) - 2*F_kernel(u0, 2*h))

                # Série infinita de imagens dielétricas de Clayton Paul
                term_sum = 0.0
                for n_idx, n in enumerate(n_arr):
                    an = 2 * n * h
                    val_n = (F_kernel(u1, an) + F_kernel(u2, an) - 2*F_kernel(u0, an))
                    term_sum += sign_k[n_idx] * val_n

                P_mat[i, j] = (P_self + (1.0 + k_val) * term_sum) / (4.0 * np.pi * EPS0 * eps_eff_val * w_sub**2)
        return P_mat

    k_fr4 = (eps_r - 1.0) / (eps_r + 1.0)
    eps_eff_fr4 = (eps_r + 1.0) / 2.0
    P_inhom = compute_P(k_fr4, eps_eff_fr4)
    P_hom   = compute_P(0.0, 1.0)

    P_inv_inhom = np.linalg.inv(P_inhom)
    P_inv_hom   = np.linalg.inv(P_hom)

    C_mat = np.zeros((n_lands, n_lands))
    C0_mat = np.zeros((n_lands, n_lands))

    for i in range(n_lands):
        idx_i = slice(i * n_sub, (i + 1) * n_sub)
        for j in range(n_lands):
            idx_j = slice(j * n_sub, (j + 1) * n_sub)
            C_mat[i, j]  = np.sum(P_inv_inhom[idx_i, idx_j])
            C0_mat[i, j] = np.sum(P_inv_hom[idx_i, idx_j])

    C_pw = np.zeros((n_lands, n_lands))
    C0_pw = np.zeros((n_lands, n_lands))
    for i in range(n_lands):
        for j in range(n_lands):
            if i == j:
                C_pw[i, i]  = np.sum(C_mat[i, :])
                C0_pw[i, i] = np.sum(C0_mat[i, :])
            else:
                C_pw[i, j]  = -C_mat[i, j]
                C0_pw[i, j] = -C0_mat[i, j]

    L_mat = MU0 * EPS0 * np.linalg.inv(C0_pw)
    return C_pw, C0_pw, L_mat

# Executa extração das matrizes canônicas
C_pw, C0_pw, L_mat = calc_mtl_matrices()

print("\n--- MATRIZ DE CAPACITÂNCIA [C] [pF/m] ---")
print(np.array2string(C_pw * 1e12, precision=3, suppress_small=True))

print("\n--- MATRIZ DE INDUTÂNCIA [L] [uH/m] ---")
print(np.array2string(L_mat * 1e6, precision=4, suppress_small=True))

# Autovalores e velocidades de fase
LC = L_mat @ C_pw
eigenvalues, T_V = np.linalg.eig(LC)
idx_sort = np.argsort(eigenvalues)
eigenvalues = eigenvalues[idx_sort]
T_V = T_V[:, idx_sort]

v_modes = 1.0 / np.sqrt(eigenvalues)
eps_eff_modes = (SPEED_OF_LIGHT / v_modes)**2
beta_modes = 2.0 * np.pi * F0 / v_modes

print("\n--- MODOS PRÓPRIOS NORMAIS FUNDAMENTAIS ---")
for m in range(3):
    print(f"Modo {m+1}: v_fase = {v_modes[m]*1e-8:.4f} x 10^8 m/s | eps_eff = {eps_eff_modes[m]:.3f} | beta = {beta_modes[m]:.2f} rad/m")
    v_norm = T_V[:, m] / np.max(np.abs(T_V[:, m]))
    print(f"        Autovetor de Tensão: [{v_norm[0]:+.3f}, {v_norm[1]:+.3f}, {v_norm[2]:+.3f}]^T")

L_beat = 2.0 * np.pi / abs(beta_modes[2] - beta_modes[0])
print(f"\nComprimento de Batimento Espacial (L_beat): {L_beat*1e3:.1f} mm ({L_beat*1e2:.2f} cm)")

# Salva gráfico de curvas paramétricas estilo Clayton Paul
s_over_h_list = np.linspace(0.2, 2.5, 12)
c11_list, c12_list, l11_list, l12_list = [], [], [], []
for soh in s_over_h_list:
    c_p, _, l_p = calc_mtl_matrices(s=soh * H_SUB, n_sub=15, n_terms=40)
    c11_list.append(c_p[0, 0] * 1e12)
    c12_list.append(abs(c_p[0, 1]) * 1e12)
    l11_list.append(l_p[0, 0] * 1e6)
    l12_list.append(l_p[0, 1] * 1e6)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.8), dpi=150)
fig.patch.set_facecolor('#0F172A')

ax1.plot(s_over_h_list, c11_list, 'o-', color='#38BDF8', lw=2, label=r'$C_{11} = C_{33}$')
ax1.plot(s_over_h_list, c12_list, 's--', color='#F43F5E', lw=2, label=r'$|C_{12}| = |C_{23}|$')
ax1.set_title('Capacitâncias por Unidade de Comp.', color='#F1F5F9', fontsize=11, fontweight='bold')
ax1.set_xlabel('Espaçamento Normalizado $s/h$', color='#94A3B8')
ax1.set_ylabel('[pF/m]', color='#94A3B8')
ax1.grid(color='#334155', linestyle=':', alpha=0.6)
ax1.tick_params(colors='#94A3B8')
ax1.legend(facecolor='#1E293B', edgecolor='#334155', labelcolor='#F1F5F9')
ax1.set_facecolor('#1E293B')

ax2.plot(s_over_h_list, l11_list, 'o-', color='#38BDF8', lw=2, label=r'$L_{11} = L_{33}$')
ax2.plot(s_over_h_list, l12_list, 's--', color='#F43F5E', lw=2, label=r'$L_{12} = L_{23}$')
ax2.set_title('Indutâncias por Unidade de Comp.', color='#F1F5F9', fontsize=11, fontweight='bold')
ax2.set_xlabel('Espaçamento Normalizado $s/h$', color='#94A3B8')
ax2.set_ylabel(r'[$\mu$H/m]', color='#94A3B8')
ax2.grid(color='#334155', linestyle=':', alpha=0.6)
ax2.tick_params(colors='#94A3B8')
ax2.legend(facecolor='#1E293B', edgecolor='#334155', labelcolor='#F1F5F9')
ax2.set_facecolor('#1E293B')

plt.suptitle("IN-Foton / LabTel — Curvas de Clayton Paul (Microfita Tripla)", color='#38BDF8', fontsize=12, fontweight='bold')
plt.tight_layout()
fig_out = os.path.join(DIR_OUTPUT, "mtl_clayton_paul_curves.png")
plt.savefig(fig_out, facecolor='#0F172A', dpi=150)
plt.close(fig)

print(f"\n✔ Gráfico salvo com sucesso em: {fig_out}")
print("="*80)