def METHOD_NAME(
    L, R, xc, yc, zc, dincl, ddecl, S, ht, f, xmin, xmax, dx, showDataPts=False
):
    L = np.array(L, dtype=float)
    R = np.array(R, dtype=float)
    xc = np.array(xc, dtype=float)
    yc = np.array(yc, dtype=float)
    zc = np.array(zc, dtype=float)
    dincl = np.array(dincl, dtype=float)
    ddecl = np.array(ddecl, dtype=float)
    S = np.array(S, dtype=float)
    ht = np.array(ht, dtype=float)
    f = np.array(f, dtype=float)
    dx = np.array(dx, dtype=float)
    ymin = xmin
    ymax = xmax
    dely = dx
    # generate the grid
    xp = np.arange(xmin, xmax, dx)
    yp = np.arange(ymin, ymax, dely)
    [y, x] = np.meshgrid(yp, xp)
    z = 0.0 * x - ht
    # set up the response arrays
    real_response = 0.0 * x
    imag_response = 0.0 * x
    # frequency characteristics
    alpha = 2.0 * np.pi * f * L / R
    f_factor = (alpha ** 2.0 + 1j * alpha) / (1 + alpha ** 2.0)
    # amin = 0.01
    # amax = 100.0
    da = 4.0 / 40.0
    alf = np.arange(-2.0, 2.0, da)
    alf = 10.0 ** alf
    fre = alf ** 2.0 / (1.0 + alf ** 2.0)
    fim = alf / (1.0 + alf ** 2.0)
    # simulate anomalies
    yt = y - S / 2.0
    yr = y + S / 2.0
    dm = -S / 2.0
    dp = S / 2.0
    M13 = mind(0.0, dm, 0.0, 90.0, 0.0, 0.0, dp, 0.0, 90.0, 0.0)
    M12 = L * mind(x, yt, z, 90.0, 0.0, xc, yc, zc, dincl, ddecl)
    M23 = L * mind(xc, yc, zc, dincl, ddecl, x, yr, z, 90.0, 0.0)
    c_response = -M12 * M23 * f_factor / (M13 * L)
    # scaled to simulate a net volumetric effect
    if np.logical_and(dincl == 0.0, ddecl == 0.0):
        real_response = np.real(c_response) * 0.0
        imag_response = np.imag(c_response) * 0.0
    else:
        real_response = np.real(c_response) * 1000.0
        imag_response = np.imag(c_response) * 1000.0
    fig, ax = plt.subplots(2, 2, figsize=(14, 8))
    ax[0][0].semilogx(alf, fre, ".-b")
    ax[0][0].semilogx(alf, fim, ".--g")
    ax[0][0].plot([alpha, alpha], [0.0, 1.0], "-k")
    ax[0][0].legend(["Real", "Imag"], loc=2)
    ax[0][0].set_xlabel("$\\alpha = \\omega L /R$")
    ax[0][0].set_ylabel("Frequency Response")
    ax[0][0].set_title("Plot 1: EM responses of loop")
    ax[0][0].grid(which="major", color="0.6", linestyle="-", linewidth="0.5")
    ax[0][0].grid(which="minor", color="0.6", linestyle="-", linewidth="0.5")
    kx = int(np.ceil(xp.size / 2.0))
    ax[0][1].plot(y[kx, :], real_response[kx, :], ".-b")  # kx
    ax[0][1].plot(y[kx, :], imag_response[kx, :], ".--g")
    # ax[0][1].legend(['Real','Imag'],loc=2)
    ax[0][1].set_xlabel("Easting")
    ax[0][1].set_ylabel("H$_s$/H$_p$")
    ax[0][1].set_title("Plot 2: EW cross section along Northing = %1.1f" % (x[kx, 0]))
    ax[0][1].grid(which="major", color="0.6", linestyle="-", linewidth="0.5")
    ax[0][1].grid(which="minor", color="0.6", linestyle="-", linewidth="0.5")
    ax[0][1].set_xlim(np.r_[xmin, xmax])
    vminR = real_response.min()
    vmaxR = real_response.max()
    ax[1][0].plot(np.r_[xp.min(), xp.max()], np.zeros(2), "k--", lw=1)
    clb = plt.colorbar(
        ax[1][0].imshow(
            real_response,
            extent=[xp.min(), xp.max(), yp.min(), yp.max()],
            vmin=vminR,
            vmax=vmaxR,
        ),
        ax=ax[1][0],
    )
    ax[1][0].set_xlim(np.r_[xmin, xmax])
    ax[1][0].set_ylim(np.r_[xmin, xmax])
    ax[1][0].set_xlabel("Easting (m)")
    ax[1][0].set_ylabel("Northing (m)")
    ax[1][0].set_title("Plot 3: Real Component")
    # ax[1][0].colorbar()
    clb.set_label("H$_s$/H$_p$")
    if showDataPts:
        XP, YP = np.meshgrid(xp, yp)
        ax[1][0].plot(XP, YP, ".", color=[0.2, 0.2, 0.2])
    vminI = imag_response.min()
    vmaxI = imag_response.max()
    ax[1][1].plot(np.r_[xp.min(), xp.max()], np.zeros(2), "k--", lw=1)
    clb = plt.colorbar(
        ax[1][1].imshow(
            imag_response,
            extent=[xp.min(), xp.max(), yp.min(), yp.max()],
            vmin=vminI,
            vmax=vmaxI,
        ),
        ax=ax[1][1],
    )
    ax[1][1].set_xlim(np.r_[xmin, xmax])
    ax[1][1].set_ylim(np.r_[xmin, xmax])
    ax[1][1].set_xlabel("Easting (m)")
    ax[1][1].set_ylabel("Northing (m)")
    ax[1][1].set_title("Plot 4: Imag Component")
    clb.set_label("H$_s$/H$_p$")
    if showDataPts:
        ax[1][1].plot(XP, YP, ".", color=[0.2, 0.2, 0.2])
    plt.tight_layout()
    plt.show()