a_j(k)\phi^j_\ell =:\; e^{\mbox{\tiny u1}}_\ell(k), \qquad\text{and} \qquad
\frac{e^{-ik\overline{z_\ell}}}{-ik}\approx \sum_{j=1}^{L-1} b_j(k)\phi^j_\ell =:\; e^{\mbox{\tiny u2}}_\ell(k),$$ where $z_\ell$ denotes the center of the $\ell$-th electrode. Choosing a scattering radius $R>0$, the integrals in are computed using a Simpson’s rule approximation $$\label{eq:Fran_scatExp_numerics}
\begin{array}{lcl}
\operatorname{S_{12}^{\mbox{\tiny{\textbf{exp}}}}}(k) & \approx & \begin{cases}
\frac{i}{4\pi}\frac{P}{L}\left[e^{-i\bar{k}\overline{\mathbf{z}}}\right]^T\left[\phi\mathbf{L}_\gamma\mathbf{e}^{\mbox{\tiny u2}}(k) + \mathbf{d}_{12}(k)\right] & 0<|k|\leq R\\
0 & |k|>R
\end{cases}\\
&&\\
\operatorname{S_{21}^{\mbox{\tiny{\textbf{exp}}}}}(k) & \approx & \begin{cases}
-\frac{i}{4\pi}\frac{P}{L}\left[e^{i\bar{k}\mathbf{z}}\right]^T\left[\phi\mathbf{L}_\gamma\mathbf{e}^{\mbox{\tiny u1}}(k) + \mathbf{d}_{21}(k)\right]& 0<|k|\leq R\\
0 & |k|>R,
\end{cases}\\
\end{array}$$ where $P$ is the perimeter of the boundary $\operatorname{\partial\Omega}$, $T$ the matrix transpose, and $\mathbf{z}$ is the vector of the centers of the electrodes $z_\ell$. The $\mathbf{d}_{12}$ and $\mathbf{d}_{21}$ terms represent the $32\operatorname{\times}1$ vectors resulting from the action of the tangential derivative map $\partial_\tau$ acting on $\frac{e^{-ik\overline{\mathbf{z}}}}{-ik}$ and $\frac{e^{ik\mathbf{z}}}{ik}$, respectively. These can be computed analytically as $$\partial_\tau\left(\frac{e^{-ik\overline{\mathbf{z}}}}{-ik}\right)=\nabla\left(\frac{e^{-ik\overline{\mathbf{z}}}}{-ik}\right)\cdot\tau=e^{-ik\bar{z}}(1,-i)\cdot(\tau_1,\tau_2) = e^{-ik\bar{z}}\overline{\tau(z)}=-i\overline{\nu(z)}e^{-ik\bar{z}},$$ and $\partial_\tau\left(\frac{e^ik\mathbf{z}}{ik}\right)=i\nu(z)e^{ikz}$. The unit normal vector $\nu(z)$ can be computed via numerical forward-differences using a parameterization $r(\theta)$ of the boundary. The unit tangent vector to the boundary can then be computed as $$\tau(\theta) = \frac{r(\theta+\varepsilon)-r(\theta)}{|r(\theta+\varepsilon)-r(\theta)|}=\tau_1(\theta)+i\tau_2(\theta),$$ for some small $\varepsilon>0$ and then the unit normal $\nu = \tau_2-i\tau_1$. Note that this step is trivially parallelizable in $k$. The scattering data for $k=0$ can be computed by interpolation.

Next, the $\operatorname{{\bar{\partial}}_k}$ system 