def METHOD_NAME(J, n_voxels, affine, calc_res=2., mask=False):
    ''' Calculates the magnetic field caused by a current density distribution J
    Parameters
    ----------
    J: simnibs.msh.ElementData
        ElementData field with the current density field
    n_voxels: list or tuple
        number of voxels in x, y, and z directions for the output B
    affine: ndarray
        A 4x4 matrix specifying the transformation from voxels to xyz for the output B
    calc_res: float (optional)
        Resolution of the space used for the FFT calculations, in mm. Default: 2
    mask: bool (optional)
        Wether to mask the resul using the mesh. Default: False
    Returns
    --------
    B: ndarray
        numpy array of size n_voxesx3 with the magnetic field        
    Reference
    ----------
    Yazdanian, H., Saturnino, G. B., Thielscher, A., & Knudsen, K. (2020).
    Fast evaluation of the Biot-Savart integral using FFT for electrical conductivity imaging.
    Journal of Computational Physics, 109408.
    https://doi.org/10.1016/j.jcp.2020.109408
    Example
    ---------
    Load a simulation result with a 'J' field can calculate B in an ROI
    Defined by a nifti file
    >>> import simnibs
    >>> import nibabel
    >>> simulation_J = simnibs.read_msh('simulation_result_with_J.msh')
    >>> reference_nifti = nibabel.load('reference.nii)
    >>> B = simnibs.simulation.calc_B(simulation_J.field['J'], reference_nifti.shape, reference_nifti.affine)
    '''
    mesh = J.mesh
    # Computational Domain (before zero-padding)
    domain = _comp_domain(mesh, n_voxels, affine)
    # Voxelize J
    J_vol, affine_vol = _voxelize(J, domain, calc_res)
    # Calculate B in the whole domain
    B_vol = _bs_ft(J_vol, calc_res)
    # Interpolate B to the domain of interest
    B = transformations.volumetric_affine(
        (B_vol, affine_vol), np.eye(4), affine, n_voxels,
        intorder=1, keep_vector_length=False
    )
    if mask:
        msk = mesh_io.ElementData(np.ones(mesh.elm.nr), mesh=mesh)
        msk = msk.interpolate_to_grid(n_voxels, affine, method='assign')
        B *= msk[..., None]
    return B