def run_simulation(wdir, arp=True, **kwargs):
    """
    Example
    -------
    100-nm silica sphere with 10-nm thick Ag coating,
    embedded in water; arprec 20 digits; illuminated with YAG (1064nm);
    scan xz plane (21x21, +-200nm)

    bhfield-arp-db.exe mpdigit wl r_core r_coat
                       n_grid_x xspan_min xspan_max
                       n_grid_y yspan_min yspan_max
                       n_grid_z zspan_min zspan_max
                       case  Kreibig
                       [n_med n_core k_core n_coat k_coat (case=other)]

    bhfield-arp-db.exe 20 1.064 0.050 0.060
                       21 -0.2 0.2
                       1 0 0
                       21 -0.2 0.2
                       other 0
                       1.3205 1.53413 0 0.565838 7.23262


    Explanation of parameters
    -------------------------
    mpdigit:
        arprec's number of precision digits;
        increase it to overcome round-off errors
    wl[um]:
        light wavelength in vacuum
    r_core[um], r_coat[um]:
        core & coat radii
    n_grid_x xspan_min[um] xspan_max[um]:
        number & span of grid points for field computation; x span
    n_grid_y yspan_min[um] yspan_max[um]:
        y span
    n_grid_z zspan_min[um] zspan_max[um]:
        z span
    Kreibig:
        Kreibig mean free path correction for Ag (0.0 - 1.0)
    case:
        nanoshell/liposome/HPC/barber/other
    n_med n_core k_core n_coat k_coat (case=other only):
        refractive indices of medium (real), core & coat (n, k)


    If `case=other`, complex refractive indices
    (n, k at the particular wavelength) must be specified.
    Otherwise (case = nanoshell etc) the medium/core/coat
    materials are predefined and the n,k values
    are taken from the data file (Ag_palik.nk etc).
    The latter reflects our own interest and is intended
    for use in our lab, so general users may not find it useful :-)
    """
    wdir = pathlib.Path(wdir)
    cmd = "{pathbhfield} {mpdigit} {wl:f} {r_core:f} {r_coat:f} " \
          + "{n_grid_x:d} {xspan_min:f} {xspan_max:f} " \
          + "{n_grid_y:d} {yspan_min:f} {yspan_max:f} " \
          + "{n_grid_z:d} {zspan_min:f} {zspan_max:f} " \
          + "{case} {Kreibig:f} {n_med:f} {n_core:f} {k_core:f} " \
          + "{n_coat:f} {k_coat:f}"
    old_dir = pathlib.Path.cwd()
    os.chdir(wdir)

    kwargs["pathbhfield"] = get_binary(arp=arp)

    if arp:
        kwargs["mpdigit"] = 16
    else:
        kwargs["mpdigit"] = ""

    # run simulation with kwargs
    sp.check_output(cmd.format(**kwargs), shell=True)

    # Go back to orgignal directory before checking (checking might fail)
    os.chdir(old_dir)

    # Check bhdebug.txt to make sure that you specify enough digits to
    # overcome roundoff errors.
    check_simulation(wdir)