eppy.py

"""Calculate eddy currents in a flat plate.

This module provides a basic simulation frameowork to calculate eddy
currents in flat, non-magnetic and isotropic, plates. The functions in
this modulecan be subdivided into three parts:
- Functions to calculate and visualize the magnetic field generated by
  a coil. The field strength is calculated using Biot-Savart.
- Functions to calculate and visualize the eddy currents in a flat
  non-magnetic plate. The eddy current distribution is calculated
  using the approach proposed by Nagel [Nagel2019].
- A text parser that allows running simulations via an input file.

Use of this module via an interactive computing environment (e.g.
Jupyter) is recommended. Alternatively the module can be run from the
command prompt to read an input file with the simulation details:

  > python eppy.py input_file.txt

This module is accompanied by the module coil_geom that can be used to
generate various coil geometries, and by a set of examples.

Nagel, J. Finite-difference simulation of eddy currents in nonmagnetic
sheets via electric vector potential. IEEE Transactions on Magnetics.
DOI: 10.1109/TMAG.2019.2940204

"""

import os
import sys
import cmath
from typing import Any, Union
import numpy as np
import numpy.typing as npt
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.integrate import nquad

from coil_geom import coil_segments


# ----------------------------------------------------------------------
# User defined types
#
ArrayFloat = npt.NDArray[np.float_]
ArrayCompl = npt.NDArray[np.complex_]


# ----------------------------------------------------------------------
# Magnetic field
#

def dB_biot_savart(dl: ArrayFloat, R: ArrayFloat, cur: float = 1.0,
                   mu_0: float = 4*np.pi*1E-7) -> ArrayFloat:
    """Return the magnetic field at R generated by wire element dl.

    Parameters
    ----------
    dl : NDArray(dtype=float, dim=1)
        Vector with the length of a wire element dl, pointing in the
        direction of the wire.
    R : NDArray(dtype=float, dim=1)
        Position vector from the wire element (dl) to the point where
        the field is evaluated.
    cur : float
        The current through the wire in Ampere (defaults to 1).
    mu_0 : float
        Magnatic permeability (the default is the magnetic
        permeability of free space).

    Returns
    -------
    dB : NDArray(dtype=float, dim=1)
        Magnetic field at position R due to a wire element (dl)
        carrying a current (I).

    """
    return (mu_0 * cur * np.cross(dl, R)) / (4*np.pi * np.linalg.norm(R)**3)


def biot_savart(dl: ArrayFloat, R0: ArrayFloat, points: ArrayFloat,
                cur: float=1.0, mu_0: float=4*np.pi*1E-7) -> ArrayFloat:
    """Return magnetic field at a list of points generated by a coil.

    Parameters
    ----------
    dl : NDArray(dtype=float, dim=2)
        Array with length vectors of the coil wire elements.
    R0 : NDArray(dtype=float, dim=2)
        Array with position vectors of the corresponding wire elements.
    points : NDArray(dtype=float, dim=2)
        Array with position vectors of the points where the magnetic
        field is to be evaluated.
    cur : float
        Current in Amperes (defaults to 1 A).
    mu_0 : float
        Magnetic permeability of the surrounding medium (defaults to
        the magnetic permeability of vacuum).

    Returns
    -------
    B : NDArray(dtype=float, dim=2)
        Array of the magnetic field components at points.

    """
    R = points[:, None, :] - R0
    cross = np.cross(dl, R)
    norm = np.linalg.norm(R, axis=2)
    dB = (mu_0 * cur * cross) / (4*np.pi * norm[:, :, None]**3)
    return np.sum(dB, axis=1)


def plot_coil(R: ArrayFloat, dl: ArrayFloat, ax: plt.Axes=None) -> plt.Axes:
    """Plot coil geometry.

    Parameters
    ----------
    R : ndarray(dtype=float, dim=2)
        Array with position vectors for line segments.
    dl : ndarray(dtype=float, dim=2)
        Array with length vectors for line segments.
    ax : matplotlib.axes, optional.
        Axes handle; if not provided, a new figure will be created.

    """
    if ax is None:
        fig = plt.figure()
        ax = fig.gca(projection='3d')
    ax.quiver(R[:, 0], R[:, 1], R[:, 2],
              dl[:, 0], dl[:, 1], dl[:, 2], pivot='middle')
    return ax



def plot_mf(x: ArrayFloat, y: ArrayFloat, B: ArrayFloat, label: str="z",
            levels: int=10, ax: plt.Axes=None) -> tuple[plt.Axes, Any]:
    """Plot magnetic field.

    Parameters
    ----------
    x : NDArray(dtype=float, dim=1)
        Coordinates in x-direction.
    y : NDArray(dtype=float, dim=1)
        Coordinates in y-direction.
    B : NDArray(dtype=float, dim=2) or NDArray(dtype=float, dim=1)
        Can be:
        - an array of vectors with magnetic field strength in x, y, z, or
        - an array with the magnetic field strength in z-direction.
    label : {'x', 'y', 'z', 'norm'}, defaults to 'z'.
        Direction of field to plot; only relevant in case B represents
        a vector field (B.ndim == 2).
    levels : int
        Number of color levels to use in the contour plot.
    ax : matplotlib.axes, optional.
        Axes handle; if not provided, a new figure will be created.

    Returns
    -------
    ax : matplotlib.axes
        Axes handle.

    """
    Nx = len(x)
    Ny = len(y)
    if B.ndim == 2:
        d = {"x": 0, "y": 1, "z": 2}[label]
        B = np.linalg.norm(B, axis=1) if d == "norm" else B[:, d]
    if ax is None:
        _, ax = plt.subplots()
    cf = ax.contourf(x, y, vector2matrix(B, Ny, Nx), levels=levels)
    ax.set_aspect('equal', adjustable='box')
    return ax, cf


# ----------------------------------------------------------------------
# Eddy currents
#

def vector2matrix(V: Union[ArrayFloat, ArrayCompl],
                  Nx: int, Ny: int) -> Union[ArrayFloat, ArrayCompl]:
    """Transform vector to matrix.

    Parameters
    ----------
    V : NDArray(dtype=float, dim=1)
        Vector to be transformed of length Nx*Ny.
    Nx : int
        Number of points in x-direction.
    Ny : int
        Number of points in y-direction.

    Returns
    -------
    M : NDArray(dtype=float, dim=2)
        Matrix.

    """
    return np.reshape(V, (Ny, Nx))


def matrix2vector(M: Union[ArrayFloat, ArrayCompl],
                  Nx: int, Ny: int) -> Union[ArrayFloat, ArrayCompl]:
    """Transform matrix to vector.

    Parameters
    ----------
    M : NDArray(dtype=float, dim=2)
        Matrix of shape (Ny, Nx)
    Nx : int
        Number of points in x-direction.
    Ny : int
        Number of points in y-direction.

    Returns
    -------
    V : NDArray(dtype=float, dim=1)
        Vector.

    """
    return np.reshape(M, Nx*Ny)


def system_matrix(rho: float, dx: float, dy: float,
                  Nx: int, Ny: int) -> ArrayFloat:
    """Return system matrix in case self-inductance is negligible.

    Returns the system matrix [M] required to calculate the electric
    vector potential {T} as a function of the magnetic flux vector
    {b}:

        [M]{T} = {b},

    following the paper of Nagel [cite:Nagel2019]. Please note that
    self-inductance is excluded in this case.

    Parameters
    ----------
    rho : float
        Resistivity of the flat plate.
    dx : float
        X-coordinate spacing.
    dy : float
        Y-coordinate spacing.
    Nx : int
        Number of grid points in X-direction.
    Ny : int
        Number of grid points in Y-direction.

    Returns
    -------
    M : NDArray(dtype=float, dim=2)
        System matrix.

    Nagel2019: Nagel, J. Finite-difference simulation of eddy currents
    in nonmagnetic sheets via electric vector potential. IEEE
    Transactions on Magnetics. DOI: 10.1109/TMAG.2019.2940204

    """
    dxdy = dx/dy
    dydx = dy/dx
    diagonals = [-rho*dxdy*np.ones(Nx*Ny-Nx),
                 -rho*dydx*np.ones(Nx*Ny-1),
                 2*rho*(dxdy+dydx)*np.ones(Nx*Ny),
                 -rho*dydx*np.ones(Nx*Ny-1),
                 -rho*dxdy*np.ones(Nx*Ny-Nx)]
    offsets = [-Nx, -1, 0, 1, Nx]
    M = diags(diagonals, offsets).toarray()
    return M


def rhs(Bz: ArrayFloat, omega: float,
        dx: float, dy: float) -> ArrayCompl:
    """Return magnetic flux vector.

    Returns the magnetic flux vector {b} required to calculate the
    electric vector potential {T} according to:

        [M]{T} = {b},

    following the work of Nagel [cite:Nagel2019]. Please note that {b}
    consists of complex numbers.

    Parameters
    ----------
    Bz : NDArray(dtype=float, dim=1)
        Z-compnent of magnetic field.
    omega : float
        Angular frequency.
    dx : float
        X-coordinate spacing.
    dy : float
        Y-coordinate spacing.

    Returns
    -------
    b : NDArray(dtype=complex, dim=1)
        Magnetic flux vector.

    Nagel2019: Nagel, J. Finite-difference simulation of eddy currents
    in nonmagnetic sheets via electric vector potential. IEEE
    Transactions on Magnetics. DOI: 10.1109/TMAG.2019.2940204

    """
    return Bz*omega*dx*dy*1j


## RENAME THIS FUNCTION TO CELL_CENTERS
def current_coordinates(x: ArrayFloat,
                        y: ArrayFloat) -> tuple[ArrayFloat, ArrayFloat]:
    """Return the coordinates of cell centers.

    Considering the illustration below, the magnetic vector potential
    is defined at the circles (O), while the currents are defined at
    the crosses (x).

         O-------O-------O-------O
         |       |       |       |
         |   x   |   x   |   x   |
         |       |       |       |
         O-------O-------O-------O
         |       |       |       |
         |   x   |   x   |   x   |
    ^ y  |       |       |       |
    |    O-------O-------O-------O
    |
    0----> x

    Parameters
    ----------
    x : NDArray(dtype=float, dim=1)
        Array X-coordinates where the electric vector potentials are
        defined.
    y : NDArray(dtype=float, dim=1)
        Array Y-coordinates where the electric vector potentials are
        defined.

    Returns
    -------
    xc : nd.array(dtype=float, dim=1)
        Array X-coordinates where the currents are defined.
    yc : nd.array(dtype=float, dim=1)
        Array Y-coordinates where the currents are defined.

    """
    xc = (x[:-1] + x[1:])/2
    yc = (y[:-1] + y[1:])/2
    return xc, yc


def derivative_matrices(dx: float, dy: float,
                        Nx: int, Ny: int) -> tuple[ArrayFloat, ArrayFloat]:
    """Return derivative matrices along X and Y.

    The derivative matrices are used to determine the current in X-
    and Y-direction from the electric vector potential {T} as:

        {Jx} = +[Dy]{T}, and: {Jy} = -[Dx]{T},

    following the work of Nagel [cite:Nagel2019].

    In more detail, and considering the illustration below, the
    derivatives at a point along X and Y are defined as:

        d/dx = (P2+P4-P1-P3)/(2 dx),

        d/dy = (P1+P2-P3-P4)/(2 dy).

    The derivative matrices are simply used to express these relations
    for all points in conveniently.

          |       |
        --P1------P2--
          |       |
          |   x   |
          |       |
        --P3------P4--
          |       |

    Nagel2019: Nagel, J. Finite-difference simulation of eddy currents
    in nonmagnetic sheets via electric vector potential. IEEE
    Transactions on Magnetics. DOI: 10.1109/TMAG.2019.2940204

    """
    N_col = Nx*Ny
    N_row = (Nx-1)*(Ny-1)
    Dx = np.zeros((N_row, N_col))
    Dy = np.zeros((N_row, N_col))
    for j in range(0, Ny-1):
        for i in range(0, Nx-1):
            Dx[j*(Nx-1)+i, j*Nx+i] = -1
            Dx[j*(Nx-1)+i, j*Nx+i+Nx] = -1
            Dx[j*(Nx-1)+i, j*Nx+i+1] = 1
            Dx[j*(Nx-1)+i, j*Nx+i+1+Nx] = 1
            Dy[j*(Nx-1)+i, j*Nx+i] = -1
            Dy[j*(Nx-1)+i, j*Nx+i+Nx] = 1
            Dy[j*(Nx-1)+i, j*Nx+i+1] = -1
            Dy[j*(Nx-1)+i, j*Nx+i+1+Nx] = 1
    Dx = Dx/2/np.float64(dx)
    Dy = Dy/2/np.float64(dy)
    return (Dx, Dy)


def mask_bc(Nx: int, Ny: int) -> npt.NDArray[np.bool_]:
    """Return a mask with False values for the domain boundaries."""
    bot = np.arange(0, Nx)
    right = np.arange(Nx-1, Nx*Ny, Nx)
    top = np.arange(Nx*Ny-Nx, Nx*Ny)
    left = np.arange(0, Nx*Ny, Nx)
    ind = np.unique(np.concatenate((bot, right, top, left), axis=None))
    mask = np.ones(Nx*Ny, dtype=bool)
    mask[ind] = False
    return mask


def plot_current_density_cf(x: ArrayFloat, y: ArrayFloat, Jx: ArrayCompl, Jy: ArrayCompl,
                            label: str="mag", ax: plt.Axes=None) -> tuple[plt.Axes, Any]:
    """Plot current density using contourf.

    Parameters
    ----------
    x : ndarray(dtype=float, dim=1)
        Coordinates in x-direction where electric vector potential T
        is defined.
    y : ndarray(dtype=float, dim=1)
        Coordinates in y-direction where electric vector potential T
        is defined.
    Jx : ndarray(dtype=float, dim=1)
        Current density in x-direction.
    Jy : ndarray(dtype=float, dim=1)
        Current density in y-direction.
    label : {'x', 'y', 'mag'}, defaults to 'mag'
        Direction of current to plot.
    ax : matplotlib.axes, optional.
        Axes handle; if not provided, a new figure will be created.

    Returns
    -------
    ax : matplotlib.axes
        Axes handle.
    cs : matplotlib.contour.QuadContourSet
        Contour set handle.

    """
    Nx, Ny = len(x), len(y)
    xc, yc = current_coordinates(x, y)
    if ax is None:
        _, ax = plt.subplots()
    if label == "mag":
        J_mag = np.sqrt(np.real(Jx)**2 + np.real(Jy)**2)
        cf = ax.contourf(xc, yc, vector2matrix(J_mag, Nx-1, Ny-1), levels=10)
    elif label == "x":
        cf = ax.contourf(xc, yc, vector2matrix(np.real(Jx), Nx-1, Ny-1), levels=10)
    elif label == "y":
        cf = ax.contourf(xc, yc, vector2matrix(np.real(Jy), Nx-1, Ny-1), levels=10)
    else:
        cf = None
    ax.set_aspect("equal", adjustable="box")
    return ax, cf


def plot_current_density(x: ArrayFloat, y: ArrayFloat, Jx: ArrayCompl, Jy: ArrayCompl,
                         label: str="mag", ax: plt.Axes=None) -> tuple[plt.Axes, Any]:
    """Plots current density using pcolormesh.

    Parameters
    ----------
    x : ndarray(dtype=float, dim=1)
        Coordinates in x-direction where electric vector potential T
        is defined.
    y : ndarray(dtype=float, dim=1)
        Coordinates in y-direction where electric vector potential T
        is defined.
    Jx : ndarray(dtype=float, dim=1)
        Current density in x-direction.
    Jy : ndarray(dtype=float, dim=1)
        Current density in y-direction.
    label : {'x', 'y', 'mag'}, defaults to 'mag'
        Direction of current to plot.
    ax : matplotlib.axes, optional.
        Axes handle; if not provided, a new figure will be created.

    Returns
    -------
    ax : matplotlib.axes
        Axes handle.
    pc : matplotlib.collections.Quadmesh
        Handle for colormesh.

    """
    Nx, Ny = len(x), len(y)
    X, Y = np.meshgrid(x, y)
    if ax is None:
        _, ax = plt.subplots()
    if label == "mag":
        J_mag = np.sqrt(np.real(Jx)**2 + np.real(Jy)**2)
        pc = ax.pcolormesh(X, Y, vector2matrix(J_mag, Nx-1, Ny-1))
    elif label == "x":
        pc = ax.pcolormesh(X, Y, vector2matrix(np.real(Jx), Nx-1, Ny-1))
    elif label == "y":
        pc = ax.pcolormesh(X, Y, vector2matrix(np.real(Jy), Nx-1, Ny-1))
    else:
        pc = None
    ax.set_aspect("equal", adjustable="box")
    return ax, pc


def plot_current_streamlines(x: ArrayFloat, y: ArrayFloat, Jx: ArrayCompl, Jy: ArrayCompl,
                             ax: plt.Axes=None) -> tuple[plt.Axes, Any]:
    """Plot current streamlines.

    Parameters
    ----------
    x : ndarray(dtype=float, dim=1)
        Coordinates in x-direction where electric vector potential T
        is defined.
    y : ndarray(dtype=float, dim=1)
        Coordinates in y-direction where electric vector potential T
        is defined.
    Jx : ndarray(dtype=float, dim=1)
        Current density in x-direction.
    Jy : ndarray(dtype=float, dim=1)
        Current density in y-direction.
    ax : matplotlib.axes, optional.
        Axes handle; if not provided, a new figure will be created.

    Returns
    -------
    ax : matplotlib.axes, optional.
        Axes handle.
    sp : matplotlib.streamplot.StreamplotSet
        Streamplot handle.

    """
    Nx, Ny = len(x), len(y)
    xc, yc = current_coordinates(x, y)
    if ax is None:
        _, ax = plt.subplots()
    sp = ax.streamplot(xc, yc,
                       vector2matrix(np.real(Jx), Nx-1, Ny-1),
                       vector2matrix(np.real(Jy), Nx-1, Ny-1),
                       density=0.6, linewidth=1, color='white')
    ax.set_aspect('equal', adjustable='box')
    return ax, sp


def contour_matrices(dx: float, dy: float, Nx: int, Ny:int,
                     omega: float) -> tuple[ArrayFloat, ArrayFloat]:
    """Return countour matrices.

    The contour matrices Cx and Cy relate the magnetic vector
    potential A to the self inductance potential p as:

        {p} = [Cx]{a_x} + [Cy]{a_y}

    with a_x and a_y denoting the x- and y-components of the magnetic
    vector potential A, following the work of Nagel [cite:Nagel2019].

    Parameters
    ----------
    dx : float
        Grid spacing in x-direction.
    dy : float
        Grid spacing in y-direction.
    Nx : int
        Number of grid points in x-direction.
    Ny : int
        Number of grid points in y-direction.
    omega : float
        Excitation frequency.

    Returns
    -------
    Cx : nd.array(dtype=float, ndim=2)
        Contour matrix in 'x-direction'.
    Cy : nd.array(dtype=float, ndim=2)
        Contour matrix in 'y-direction'.

    Nagel2019: Nagel, J. Finite-difference simulation of eddy currents
    in nonmagnetic sheets via electric vector potential. IEEE
    Transactions on Magnetics. DOI: 10.1109/TMAG.2019.2940204

    """
    N_col = (Nx-1)*(Ny-1)
    N_row = Nx*Ny
    Cx = np.zeros((N_row, N_col), dtype='cfloat')
    Cy = np.zeros((N_row, N_col), dtype='cfloat')
    for i in range(1, Nx-1):
        for j in range(1, Ny-1):
            Cx[j*Nx+i, (j-1)*(Nx-1)+i] = 1j
            Cx[j*Nx+i, j*(Nx-1)+i] = -1j
            Cx[j*Nx+i, j*(Nx-1)+i-1] = -1j
            Cx[j*Nx+i, (j-1)*(Nx-1)+i-1] = 1j
            Cy[j*Nx+i, (j-1)*(Nx-1)+i] = 1j
            Cy[j*Nx+i, j*(Nx-1)+i] = 1j
            Cy[j*Nx+i, j*(Nx-1)+i-1] = -1j
            Cy[j*Nx+i, (j-1)*(Nx-1)+i-1] = -1j
    return omega*dx*Cx/2, omega*dy*Cy/2


def volume_int_I(dx: float, dy: float, t: float) -> tuple[float, float]:
    """Return volume integral needed for Biot Savart matrices.

    Following the work of Nagel [cite:Nagel2019], the volume integral
    is defined as:

        I = 8 * Int 1/sqrt(x^2 + y^2 + z^2) dxdydz,

    with the bounds:
        x : (0, dx/2)
        y : (0, dy/2)
        z : (0, t/2)

    The value of I is obtained through numerical integration.

    Parameters
    ----------
    dx : float
        Grid spacing in x-direction.
    dy : float
        Grid spacing in y-direction.
    t : float
        Plate thickness.

    Returns
    -------
    I : float
        Volume integral.
    err : float
        Estimate of absolute error.

    Nagel2019: Nagel, J. Finite-difference simulation of eddy currents
    in nonmagnetic sheets via electric vector potential. IEEE
    Transactions on Magnetics. DOI: 10.1109/TMAG.2019.2940204

    """
    def f(x, y, z):
        return 1/np.sqrt(x**2 + y**2 + z**2)
    res = nquad(f, [[0, dx/2],
                    [0, dy/2],
                    [0, t/2]])
    I = res[0]
    err = res[1]
    return 8*I, err


def biot_savart_matrix(x: ArrayFloat, y: ArrayFloat, t: float,
                       mu_0=4*np.pi*1E-7) -> ArrayFloat:
    """Return Biot-Savart matrix.

    The Biot-Savart matrix relates the eddy current density J to the
    magnetic vector A as:

        {a_x} = [N]{J_x},

        {a_y} = [N]{J_y},

    with a_x and a_y denoting the x- and y-components of the magnetic
    vector potential A, and J_x and J_y denoting the x- and
    y-components of the eddy current density.

    Parameters
    ----------
    x : nd.array(dtype=float, dim=1)
        Array X-coordinates where the electric vector potentials are
        defined.
    y : nd.array(dtype=float, dim=1)
        Array Y-coordinates where the electric vector potentials are
        defined.
    t : float
        Plate thickness.
    mu_0 : float, optional
        Magnatic permeability (the default is the magnetic
        permeability of free space).


    Returns
    -------
    I : float
        Volume integral.
    err : float
        Estimate of absolute error.

    Nagel2019: Nagel, J. Finite-difference simulation of eddy currents
    in nonmagnetic sheets via electric vector potential. IEEE
    Transactions on Magnetics. DOI: 10.1109/TMAG.2019.2940204

    """
    dx = x[1]-x[0]
    dy = y[1]-y[0]
    xc, yc = current_coordinates(x, y)
    pos = np.array([np.array([x, y, 0]) for y in yc for x in xc])
    N = np.zeros((len(xc)*len(yc), len(xc)*len(yc)))
    dV = dx*dy*t
    I, _ = volume_int_I(dx, dy, t)
    for i, loc in enumerate(pos):
        with np.errstate(divide='ignore'):
            N[i] = dV/np.linalg.norm(loc-pos, axis=1)
        N[i, i] = I
    return (mu_0/4/np.pi)*N


def phase_shift(J: ArrayCompl, phi: float) -> ArrayCompl:
    """Apply phase shift to array with complex numbers.

    Parameters
    ----------
    J : nd.array(dtype=complex, dim=1)
        Array with complex numbers to which the phase shift is applied.
    phi : float
        Phase shift angle (in radians).

    Returns
    -------
    Js : nd.array(dtype=complex, dim=1)
        Phase shifted numbers.

    """
    r = np.abs(J)
    theta = np.angle(J) + phi
    Js = np.array([cmath.rect(r[i], theta[i]) for i in range(len(r))])
    return Js


# ----------------------------------------------------------------------
# Input file parser
#

def parse_line(line: str) -> tuple[str, list[str]]:
    """Return command and arguments from line from input file."""
    ln = line.split(";", 1)[0].split(",")
    if bool(ln[0].split()):
        command = ln[0].split()[0]
        args = [arg.split()[0] for arg in ln[1:]]
    else:
        command = 'pass'
        args = []
    return command, args


def parse_file(fn: str) -> tuple[dict, dict]:
    """Return coil and plate dictionary from input file.

    Parameters
    ----------
    fn : string
        Path to file.

    Returns
    -------
    coil : dict
        Dictionary with coil information. Has the following keys:
        freq, amplitiude, points, lines, circles, arcs, esize.
    plate : dict
        Dictionary with plate information. Has the following keys: Lx,
        Ly, dx, dy, thickness, cond.

    """
    points = []
    lines = []
    arcs = []
    circles = []
    freq = amplitude = esize = 0.0
    Lx = Ly = thickness = dx = dy = cond = 0.0

    with open(fn) as f:
        lns = f.readlines()
        for ln in lns:
            command, args = parse_line(ln)
            if command == 'p':
                points.append(args)
            elif command == 'line':
                lines.append(args)
            elif command == 'arc':
                arcs.append(args)
            elif command == 'circle':
                circles.append(args)
            elif command == 'esize':
                esize = float(args[0])
            elif command == 'freq':
                freq = float(args[0])
            elif command == 'amplitude':
                amplitude = float(args[0])
            elif command == 'lenx':
                Lx = float(args[0])
            elif command == 'leny':
                Ly = float(args[0])
            elif command == 'thickness':
                thickness = float(args[0])
            elif command == 'dx':
                dx = float(args[0])
            elif command == 'dy':
                dy = float(args[0])
            elif command == 'cond':
                cond = float(args[0])

    coil = {"freq": freq,
            "amplitude": amplitude,
            "points": np.array(points, dtype=float),
            "lines": np.array(lines, dtype=int),
            "circles": np.array(circles, dtype=int),
            "arcs": np.array(arcs, dtype=int),
            "esize": esize}
    plate = {"Lx": Lx, "Ly": Ly,
             "dx": dx, "dy": dy,
             "thickness": thickness, "cond": cond}
    return coil, plate


def run_input_file(fn: str) -> None:
    """Run simulation based on input file.

    Parameters
    ----------
    fn : string
        Path to file.

    """
    print("Welcome to eppy's eddy current calculator.")
    print("------------------------------------------")

    coil, plate = parse_file(fn)
    print("Input file loaded.")

    # plate dimensions and cell size
    Lx, Ly = plate["Lx"], plate["Ly"]
    dx, dy = plate["dx"], plate["dy"]
    t = plate["thickness"]

    # position vector for points on XY plane
    Nx = int(np.ceil(Lx/dx + 1))
    Ny = int(np.ceil(Ly/dy + 1))
    X = np.linspace(-Lx/2, Lx/2, Nx)
    Y = np.linspace(-Ly/2, Ly/2, Ny)
    pos = np.array([np.array([x, y, 0]) for y in Y for x in X])

    # resistivity
    rho = 1/plate["cond"]

    # coil excitation frequency and current
    omega = 2*np.pi*coil["freq"]
    current = coil["amplitude"]

    # system matrix
    M = system_matrix(rho, dx, dy, Nx, Ny)
    N = biot_savart_matrix(X, Y, t)
    Cx, Cy = contour_matrices(dx, dy, Nx, Ny, omega)
    Dx, Dy = derivative_matrices(dx, dy, Nx, Ny)
    K = M + Cx@N@Dy - Cy@N@Dx

    # unknown electric vector potential and mask
    T = np.zeros(Nx*Ny, dtype=complex)
    print("System matrices generated.")

    # coil
    R, dl = coil_segments(coil["points"], coil["esize"],
                          lines=coil["lines"],
                          circles=coil["circles"],
                          arcs=coil["arcs"])

    # magnetic field
    B = biot_savart(dl, R, pos, current)
    Bz = B[:, 2]
    flux = rhs(Bz, omega, dx, dy)
    print("Magnetic field and flux calculated.")

    # solve system
    mask = mask_bc(Nx, Ny)
    T[mask] = np.linalg.solve(K[:, mask][mask, :], flux[mask])

    # calculate currents
    Jx = np.dot(Dy, T)
    Jy = -np.dot(Dx, T)
    print("Currents calculated.")

    # plot Z-component of coil magnetic field and eddy current distr.
    fig, ax = plt.subplots(nrows=1, ncols=2, squeeze=True, figsize=(12, 6))
    _, _ = plot_mf(X, Y, Bz, label='z', levels=10, ax=ax[0]) # type: ignore
    _, cs_I = plot_current_density(X, Y, Jx, Jy, label="mag", ax=ax[1]) # type: ignore
    _, _ = plot_current_streamlines(X, Y, Jx, Jy, ax=ax[1]) # type: ignore

    # labels
    ax[0].set_title("Z-component of magnetic field (coil)")
    ax[0].set_xlabel("x [m]")
    ax[0].set_ylabel("y [m]")
    ax[1].set_title("Eddy current distribution [A/m^2]")
    ax[1].set_xlabel("x [m]")
    ax[1].set_ylabel("y [m]")

    # add color bar
    fig.subplots_adjust(right=0.8)
    cbar_ax = fig.add_axes([0.85, 0.25, 0.05, 0.5])
    fig.colorbar(cs_I, cax=cbar_ax)

    # limits
    ax[0].set_xlim([-Lx/2, Lx/2])
    ax[0].set_ylim([-Ly/2, Ly/2])
    ax[1].set_xlim([-Lx/2, Lx/2])
    ax[1].set_ylim([-Ly/2, Ly/2])

    # show plot
    fig_fn = os.path.splitext(fn)[0] + ".png"
    plt.savefig(fig_fn, dpi=300)
    plt.show()
    print("All done! Figure is saved as '{}'.".format(fig_fn))


if __name__ == "__main__":
    fn = sys.argv[1]
    run_input_file(fn)