Repository: magpylib/magpylib (language: python)

Issue:
review edgesize setting
this is another useless feature.

The field expressions should be stabilized in such a way that singularities just happen at corners / edges etc.


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d928af0a3..acf5199c1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,17 +3,19 @@ All notable changes to magpylib are documented here.
 
 # Releases
 
-## [Unreleased]
+## [Unreleased v4]
 This is a major update that includes
 
 - API changes
 - New features
 - Improved internal workings
 ### Added
+- New computaion core: ([#376](https://github.com/magpylib/magpylib/issues/376))
+  - top level subpackage `magpylib.core` where all field expression implementations can be accessed directly.
 - New sources classes:
-  - New `CylinderSegment` computation with sector angles `(r1,r2,h,phi1,phi2)`. ([#386](https://github.com/magpylib/magpylib/issues/386), [#385](https://github.com/magpylib/magpylib/issues/385))
+  - New `CylinderSegment` computation with dimension `(r1,r2,h,phi1,phi2)` with the inner radius `r1`, the outer radius `r2` the height `h` and the cylinder section angles `phi1 < phi2`. ([#386](https://github.com/magpylib/magpylib/issues/386), [#385](https://github.com/magpylib/magpylib/issues/385))
   - New `CustomSource` class for user defined field functions ([#349](https://github.com/magpylib/magpylib/issues/349), [#409](https://github.com/magpylib/magpylib/issues/409), [#411](https://github.com/magpylib/magpylib/pull/411))
-- New display features:
+- New `display` (now `show`) features:
   - 🚀 New `plotly` plotting backend: ([#396](https://github.com/magpylib/magpylib/pull/396), [#353](https://github.com/magpylib/magpylib/issues/353))
     - `plotly` remains a optional dependency, needs extra installation ([#395](https://github.com/magpylib/magpylib/issues/395))
     - Interactive `animation` feature with `animation_time` and `animation_fps` `display` arguments
@@ -34,11 +36,11 @@ This is a major update that includes
 ### Changed
 - `Box` class renamed to `Cuboid`. ([#350](https://github.com/magpylib/magpylib/issues/350))
 - `Circular` class renamed to `Loop`. ([#402](https://github.com/magpylib/magpylib/pull/402))
-- `magpylib.defaults.edgesize` set to `1e-8` by default to avoid problems in `Cuboid` corners.
 - Magpylib objects can now be _initialized_ **without excitation** and **without dimension** attributes.
 - `magpylib.Config` parameters are now in `magpylib.defaults` ([#387](https://github.com/magpylib/magpylib/issues/387))
 - `getBv` and `getHv` are now integrated into `getB` and `getH` ([#449](https://github.com/magpylib/magpylib/pull/449))
-- Renamed `display` arguments
+- `display` is renamed to `show` ([#453](https://github.com/magpylib/magpylib/pull/453), [#451](https://github.com/magpylib/magpylib/issues/451))
+- Renamed `show` (prev. `display`) arguments
   - `axis` ➡️ `canvas`
   - `show_direction` ➡️ `style_magnetization_show`
   - `show_path` ➡️ `path`
@@ -46,14 +48,21 @@ This is a major update that includes
   - `size_direction` ➡️ `style_magnetization_size`
 - `move` and `rotate` methods have `start='auto'` instead of `start=-1`. Apply to full path if input is a scalar, append/merge if input is a vector ([#438](https://github.com/magpylib/magpylib/discussions/438), [#444](https://github.com/magpylib/magpylib/issues/444))
 - `show_path` in the `display` function is integrated in styling options with `style_path_show` ([#453](https://github.com/magpylib/magpylib/pull/453))
-- `display` is renamed to `show` ([#453](https://github.com/magpylib/magpylib/pull/453), [#451](https://github.com/magpylib/magpylib/issues/451))
 
 ### Updated
 - Computation:
-  - Updated `Cylinder` computation with a new exact closed form implementation without iteration, from a new paper (F.Slanovc, preprint-2021). Computation times are around 50-100 µs, with some performance improvement planned in the future.
-  - Added `__len__` dunder for `Collection` ([#383](https://github.com/magpylib/magpylib/issues/383))
+  - Improoved field computations: ([#374](https://github.com/magpylib/magpylib/issues/374))
+    - negative dimension input taken as absolute when only positive dimensions are allowed
+    - scale invariant field evaluations
+    - special cases caught within 1e-15 rtol and atol to account for numerical imprecision with positioning (e.g. object rotation).
+    - supress numpy divide/invalid warnings. return np.nan as (0,0,0) (e.g. on magnet edges or on line currents) and allow return of np.inf.
+  - Updated `Cylinder` diametral magnetization computation with novel closed form implementation. 
+    - Much faster (100-1000x)
+    - Numerically stable for small `r`.
   - Improved performance of `getB` and `getH` functions and methods.
-  - Improved numerical stability ([#374](https://github.com/magpylib/magpylib/issues/374))
+  - Improved numerical stability of current Loop field. [#374](https://github.com/magpylib/magpylib/issues/374) 
+- Other
+  - Added `__len__` dunder for `Collection` ([#383](https://github.com/magpylib/magpylib/issues/383))
 - Docs:
   - General docs improvements ([#399](https://github.com/magpylib/magpylib/issues/399), [#294](https://github.com/magpylib/magpylib/issues/294))
   - [MATLAB example](https://magpylib.readthedocs.io/en/latest/_pages/3_MATLAB/#example) updated ([#346](https://github.com/magpylib/magpylib/issues/346), [#366](https://github.com/magpylib/magpylib/pull/366))
diff --git a/docs/_pages/0_documentation.md b/docs/_pages/0_documentation.md
index 288508e26..d532dff46 100644
--- a/docs/_pages/0_documentation.md
+++ b/docs/_pages/0_documentation.md
@@ -354,7 +354,7 @@ mag = np.array([(100,0,0)]*5)
 dim = np.array([(1,2,45,90,-1,1)]*5)
 obs_pos = np.array([(0,0,0)]*5)
 
-B = magpy.lib.magnet_cyl_tile_H_Slanovc2021(mag, dim, obs_pos)
+B = magpy.lib.magnet_cylinder_section_core(mag, dim, obs_pos)
 print(B)
 
 # out: [[   0.           0.        -186.1347833]
diff --git a/magpylib/__init__.py b/magpylib/__init__.py
index a81d1a8a9..27a95f3fc 100644
--- a/magpylib/__init__.py
+++ b/magpylib/__init__.py
@@ -297,7 +297,7 @@
 >>> dim = np.array([(1,2,45,90,-1,1)]*5)
 >>> poso = np.array([(0,0,0)]*5)
 >>>
->>> B = magpy.lib.magnet_cyl_tile_H_Slanovc2021(mag, dim, poso)
+>>> B = magpy.lib.magnet_cylinder_segment_core(mag, dim, poso)
 >>> print(B)
 [[   0.           0.        -186.1347833]
  [   0.           0.        -186.1347833]
@@ -315,14 +315,11 @@
 __credits__ = 'Silicon Austria Labs - Sensor Systems'
 __all__ = ['magnet', 'current', 'misc', 'getB', 'getH',
     'Sensor', 'Collection', 'show', 'display', 'defaults', '__version__',
-    '__author__', '__credits__']
+    '__author__', '__credits__', 'core']
 
 # create interface to outside of package
+from magpylib import magnet, current, misc, display, core
 from magpylib._src.defaults.defaults_classes import default_settings as defaults
-from magpylib import magnet
-from magpylib import current
-from magpylib import misc
-from magpylib import display
 from magpylib._src.fields import getB, getH
 from magpylib._src.obj_classes import Sensor
 from magpylib._src.obj_classes.class_Collection import Collection
diff --git a/magpylib/_src/defaults/defaults_values.py b/magpylib/_src/defaults/defaults_values.py
index 4912b5140..197908278 100644
--- a/magpylib/_src/defaults/defaults_values.py
+++ b/magpylib/_src/defaults/defaults_values.py
@@ -2,8 +2,6 @@
 
 DEFAULTS = {
     "checkinputs": True,
-    "edgesize": 1e-8,
-    "itercylinder": 50,
     "display": {
         "autosizefactor": 10,
         "animation": {"fps":20, "maxfps": 30, "maxframes": 200, "time": 5, "slider": True},
diff --git a/magpylib/_src/display/display_matplotlib.py b/magpylib/_src/display/display_matplotlib.py
index 40404733e..aec8671dc 100644
--- a/magpylib/_src/display/display_matplotlib.py
+++ b/magpylib/_src/display/display_matplotlib.py
@@ -14,7 +14,7 @@
     faces_cylinder,
     system_size,
     faces_sphere,
-    faces_cylinder_section,
+    faces_cylinder_segment,
 )
 from magpylib._src.input_checks import check_excitations
 from magpylib._src.style import get_style
@@ -394,7 +394,7 @@ def display_matplotlib(
                 elif obj._object_type == "Cylinder":
                     faces = faces_cylinder(obj, path_frames)
                 elif obj._object_type == "CylinderSegment":
-                    faces = faces_cylinder_section(obj, path_frames)
+                    faces = faces_cylinder_segment(obj, path_frames)
                 elif obj._object_type == "Sphere":
                     faces = faces_sphere(obj, path_frames)
                 elif obj._object_type == "Line":
diff --git a/magpylib/_src/display/display_utility.py b/magpylib/_src/display/display_utility.py
index 0c5e77a1e..194ac9959 100644
--- a/magpylib/_src/display/display_utility.py
+++ b/magpylib/_src/display/display_utility.py
@@ -264,7 +264,7 @@ def faces_cylinder(src, show_path):
     return all_faces
 
 
-def faces_cylinder_section(src, show_path):
+def faces_cylinder_segment(src, show_path):
     """
     Compute vertices and faces of CylinderSegment for plotting.
 
@@ -280,7 +280,7 @@ def faces_cylinder_section(src, show_path):
     # pylint: disable=protected-access
     res = 15  # surface discretization
 
-    # generate cylinder section faces
+    # generate cylinder segment faces
     r1, r2, h, phi1, phi2 = src.dimension
     res_tile = (
         int((phi2 - phi1) / 360 * 2 * res) + 2
diff --git a/magpylib/_src/display/plotly/plotly_display.py b/magpylib/_src/display/plotly/plotly_display.py
index 39cebcc91..98e972fa4 100644
--- a/magpylib/_src/display/plotly/plotly_display.py
+++ b/magpylib/_src/display/plotly/plotly_display.py
@@ -519,9 +519,7 @@ def get_plotly_traces(
             )
             make_func = make_Cuboid
         elif isinstance(input_obj, Cylinder):
-            base_vertices = min(
-                50, Config.itercylinder
-            )  # no need to render more than 50 vertices
+            base_vertices = 50
             kwargs.update(
                 mag=input_obj.magnetization,
                 diameter=input_obj.dimension[0],
@@ -530,9 +528,7 @@ def get_plotly_traces(
             )
             make_func = make_Cylinder
         elif isinstance(input_obj, CylinderSegment):
-            vert = min(
-                50, Config.itercylinder
-            )  # no need to render more than 50 vertices
+            vert = 50
             kwargs.update(
                 mag=input_obj.magnetization, dimension=input_obj.dimension, vert=vert,
             )
diff --git a/magpylib/_src/fields/__init__.py b/magpylib/_src/fields/__init__.py
index 0adb00cd0..560776530 100644
--- a/magpylib/_src/fields/__init__.py
+++ b/magpylib/_src/fields/__init__.py
@@ -1,7 +1,6 @@
 """_src.fields"""
 
-__all__ = ['getB', 'getB_dict', 'getH', 'getH_dict']
+__all__ = ['getB', 'getH']
 
 # create interface to outside of package
 from magpylib._src.fields.field_wrap_BH_level3 import getB, getH
-from magpylib._src.fields.field_wrap_BH_level2_dict import getB_dict, getH_dict
diff --git a/magpylib/_src/fields/field_BH_cuboid.py b/magpylib/_src/fields/field_BH_cuboid.py
index b8139d35a..1639b5ae3 100644
--- a/magpylib/_src/fields/field_BH_cuboid.py
+++ b/magpylib/_src/fields/field_BH_cuboid.py
@@ -4,111 +4,42 @@
 """
 
 import numpy as np
-from magpylib._src.defaults.defaults_classes import default_settings as Config
-
-
-def field_BH_cuboid(
-        bh: bool,
-        mag: np.ndarray,
-        dim: np.ndarray,
-        pos_obs: np.ndarray
-        ) -> np.ndarray:
-    """ setting up the Cuboid field computation
-    - separate mag=0 cases (returning 0)
-    - separate edge/corner cases (returning 0)
-    - call field computation for general cases
-    - select B or H
-    - transform B->H (inside check)
-
-    ### Args:
-    - bh (boolean): True=B, False=H
-    - mag (ndarray Nx3): homogeneous magnetization vector in units of mT
-    - dim (ndarray Nx3): dimension of Cuboid side lengths in units of mm
-    - pos_obs (ndarray Nx3): position of observer in units of mm
-
-    ### Returns:
-    - B/H-field (ndarray Nx3): magnetic field vectors at pos_obs in units of mT / kA/m
-    """
-
-    edgesize = Config.edgesize
-
-    # allocate field vectors ----------------------
-    B = np.zeros((len(mag),3))
-
-    # special case mag = 0 ------------------------
-    mask0 = (mag[:,0]==0) * (mag[:,1]==0) * (mag[:,2]==0)
-
-    # special cases for edge/corner fields --------
-    x, y, z = np.copy(pos_obs.T)
-    a, b, c = dim.T/2
-
-    mx1 = (abs(abs(x)-a) < edgesize)
-    my1 = (abs(abs(y)-b) < edgesize)
-    mz1 = (abs(abs(z)-c) < edgesize)
-
-    mx2 = (abs(x)-a < edgesize) # within actual edge
-    my2 = (abs(y)-b < edgesize)
-    mz2 = (abs(z)-c < edgesize)
-
-    mask_xedge = my1 & mz1 & mx2
-    mask_yedge = mx1 & mz1 & my2
-    mask_zedge = mx1 & my1 & mz2
-    mask_edge = mask_xedge | mask_yedge | mask_zedge
-
-    # not a special case --------------------------
-    mask_gen = ~mask_edge & ~mask0
-
-    # compute field -------------------------------
-    if np.any(mask_gen):
-        B[mask_gen] = magnet_cuboid_B_Yang1999(mag[mask_gen], dim[mask_gen], pos_obs[mask_gen])
-
-    # return B or compute and return H -------------
-    if bh:
-        return B
-
-    # if inside magnet subtract magnetization vector
-    poso_abs = np.abs(pos_obs)
-    case1 = poso_abs[:,0]<=dim[:,0]/2-edgesize
-    case2 = poso_abs[:,1]<=dim[:,1]/2-edgesize
-    case3 = poso_abs[:,2]<=dim[:,2]/2-edgesize
-    mask_inside = case1*case2*case3
-    B[mask_inside] -= mag[mask_inside]
-    # transform units mT -> kA/m
-    H = B*10/4/np.pi
-    return H
+from magpylib._src.input_checks import check_field_input
 
 
 # ON INTERFACE
-def magnet_cuboid_B_Yang1999(
-    mag: np.ndarray,
-    dim: np.ndarray,
-    pos_obs: np.ndarray) -> np.ndarray:
+def magnet_cuboid_field(
+    magnetization: np.ndarray,
+    dimension: np.ndarray,
+    observer: np.ndarray,
+    field='B') -> np.ndarray:
     """
-    B-field in Cartesian CS of Cuboid magnet with homogenous magnetization.
-    The Cuboid sides are parallel to the CS axes.
-    The geometric center of the Cuboid lies in the origin.
-
-    Implementation from [Yang1999], [Engel-Herbert2005], [Camacho2013], [Cichon2019].
+    Computes the magnetic field of a homogeneously magnetized cuboid magnet
+    in Cartesian coordinates . The cuboid sides are parallel to the coordinate axes.
+    The geometric center of the cuboid lies in the origin.
 
     Parameters
     ----------
-    mag: ndarray, shape (n,3)
-        Homogeneous magnetization vector in units of [mT]
+    magnetization: ndarray, shape (n,3)
+        Homogeneous magnetization vector in units of [mT].
+
+    dimension: ndarray, shape (n,3)
+        Cuboid side lengths in units of [mm]. Positive input expected.
 
-    dim: ndarray, shape (n,3)
-        Cuboid side lengths in units of [mm]
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
 
-    pos_obs: ndarray, shape (n,3)
-        position of observer in units of [mm]
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
 
     Returns
     -------
-    B-field: ndarray
-        B-field of Cuboid in Cartesian CS, shape (n,3) in units of [mT].
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of magnet in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
 
     Examples
     --------
-    Compute the field of three instances.
+    Compute the field of three different instances.
 
     >>> import numpy as np
     >>> import magpylib as magpy
@@ -126,7 +57,7 @@ def magnet_cuboid_B_Yang1999(
     Field computations via magnetic surface charge density. Published
     several times with similar expressions:
 
-    Camacho: Revista Mexicana de F´ısica E 59 (2013) 8–17
+    Camacho: Revista Mexicana de Fisica E 59 (2013) 8–17
 
     Engel-Herbert: Journal of Applied Physics 97(7):074504 - 074504-4 (2005)
 
@@ -144,105 +75,160 @@ def magnet_cuboid_B_Yang1999(
     positions to their bottQ4 counterparts. see also
 
     Cichon: IEEE Sensors Journal, vol. 19, no. 7, April 1, 2019, p.2509
-
-    Numerical instabilities: see GitHub discussion
     """
     # pylint: disable=too-many-statements
 
-    magx, magy, magz = mag.T
-    x, y, z = pos_obs.T
-    a, b, c = dim.T/2
-    n = len(magx)
-
-    # avoid indeterminate forms by evaluating in bottQ4 only --------
-    # basic masks
-    maskx = x<0
-    masky = y>0
-    maskz = z>0
-
-    # change all positions to their bottQ4 counterparts
-    x[maskx] = x[maskx]*-1
-    y[masky] = y[masky]*-1
-    z[maskz] = z[maskz]*-1
-
-    # create sign flips for position changes
-    qsigns = np.ones((n,3,3))
-    qs_flipx = np.array([[ 1,-1,-1],[-1, 1, 1],[-1, 1, 1]])
-    qs_flipy = np.array([[ 1,-1, 1],[-1, 1,-1],[ 1,-1, 1]])
-    qs_flipz = np.array([[ 1, 1,-1],[ 1, 1,-1],[-1,-1, 1]])
-    # signs flips can be applied subsequently
-    qsigns[maskx] = qsigns[maskx]*qs_flipx
-    qsigns[masky] = qsigns[masky]*qs_flipy
-    qsigns[maskz] = qsigns[maskz]*qs_flipz
-
-    # field computations --------------------------------------------
-    # Note: in principle the computation for all three mag-components can be
-    #   vectorized itself using symmetries. However, tiling the three
-    #   components will cost more than is gained by the vectorized evaluation
-
-    # Note: making the following computation steps is not necessary
-    #   as mkl will cache such small computations
-    xma, xpa = x-a, x+a
-    ymb, ypb = y-b, y+b
-    zmc, zpc = z-c, z+c
-
-    xma2, xpa2 = xma**2, xpa**2
-    ymb2, ypb2 = ymb**2, ypb**2
-    zmc2, zpc2 = zmc**2, zpc**2
-
-    mmm = np.sqrt(xma2 + ymb2 + zmc2)
-    pmp = np.sqrt(xpa2 + ymb2 + zpc2)
-    pmm = np.sqrt(xpa2 + ymb2 + zmc2)
-    mmp = np.sqrt(xma2 + ymb2 + zpc2)
-    mpm = np.sqrt(xma2 + ypb2 + zmc2)
-    ppp = np.sqrt(xpa2 + ypb2 + zpc2)
-    ppm = np.sqrt(xpa2 + ypb2 + zmc2)
-    mpp = np.sqrt(xma2 + ypb2 + zpc2)
-
-    ff2x = np.log((xma+mmm) * (xpa+ppm) * (xpa+pmp) * (xma+mpp))     \
-            -np.log((xpa+pmm) * (xma+mpm) * (xma+mmp) * (xpa+ppp))
-
-    ff2y = np.log((-ymb+mmm) * (-ypb+ppm) * (-ymb+pmp) * (-ypb+mpp)) \
-           -np.log((-ymb+pmm) * (-ypb+mpm) * ( ymb-mmp) * ( ypb-ppp))
-
-    ff2z = np.log((-zmc+mmm) * (-zmc+ppm) * (-zpc+pmp) * (-zpc+mpp)) \
-           -np.log((-zmc+pmm) * ( zmc-mpm) * (-zpc+mmp) * ( zpc-ppp))
-
-    ff1x = (np.arctan2((ymb*zmc),(xma*mmm)) - np.arctan2((ymb*zmc),(xpa*pmm))
-            - np.arctan2((ypb*zmc),(xma*mpm)) + np.arctan2((ypb*zmc),(xpa*ppm))
-            - np.arctan2((ymb*zpc),(xma*mmp)) + np.arctan2((ymb*zpc),(xpa*pmp))
-            + np.arctan2((ypb*zpc),(xma*mpp)) - np.arctan2((ypb*zpc),(xpa*ppp)))
-
-    ff1y = (np.arctan2((xma*zmc),(ymb*mmm)) - np.arctan2((xpa*zmc),(ymb*pmm))
-            - np.arctan2((xma*zmc),(ypb*mpm)) + np.arctan2((xpa*zmc),(ypb*ppm))
-            - np.arctan2((xma*zpc),(ymb*mmp)) + np.arctan2((xpa*zpc),(ymb*pmp))
-            + np.arctan2((xma*zpc),(ypb*mpp)) - np.arctan2((xpa*zpc),(ypb*ppp)))
-
-    ff1z = (np.arctan2((xma*ymb),(zmc*mmm)) - np.arctan2((xpa*ymb),(zmc*pmm))
-            - np.arctan2((xma*ypb),(zmc*mpm)) + np.arctan2((xpa*ypb),(zmc*ppm))
-            - np.arctan2((xma*ymb),(zpc*mmp)) + np.arctan2((xpa*ymb),(zpc*pmp))
-            + np.arctan2((xma*ypb),(zpc*mpp)) - np.arctan2((xpa*ypb),(zpc*ppp)))
-
-    # contributions from x-magnetization
-    bx_magx = magx * ff1x * qsigns[:,0,0]  # the 'missing' third sign is hidden in ff1x
-    by_magx = magx * ff2z * qsigns[:,0,1]
-    bz_magx = magx * ff2y * qsigns[:,0,2]
-    # contributions from y-magnetization
-    bx_magy =  magy * ff2z * qsigns[:,1,0]
-    by_magy =  magy * ff1y * qsigns[:,1,1]
-    bz_magy = -magy * ff2x * qsigns[:,1,2]
-    # contributions from z-magnetization
-    bx_magz =  magz * ff2y * qsigns[:,2,0]
-    by_magz = -magz * ff2x * qsigns[:,2,1]
-    bz_magz =  magz * ff1z * qsigns[:,2,2]
-
-    # summing all contributions
-    bx_tot = bx_magx + bx_magy + bx_magz
-    by_tot = by_magx + by_magy + by_magz
-    bz_tot = bz_magx + bz_magy + bz_magz
-
-    # combine with special edge/corner cases
-    # B = np.c_[bx_tot, by_tot, bz_tot]      # faster for 10^5 and more evaluations
-    B = np.concatenate(((bx_tot,),(by_tot,),(bz_tot,)), axis=0).T
-
-    return B / (4*np.pi)
+    bh = check_field_input(field, 'magnet_cuboid_field()')
+
+    magx, magy, magz = magnetization.T
+    a, b, c = np.abs(dimension.T)/2
+    x, y, z = observer.T
+
+    # This implementation is completely scale invariant as only observer/dimension
+    # ratios appear in equations below.
+
+    # dealing with special cases -----------------------------------
+
+    # allocate B with zeros
+    B_all = np.zeros((len(magx),3))
+
+    # SPECIAL CASE 1: mag = (0,0,0)
+    mask1 = (magx==0) * (magy==0) * (magz==0) # 2x faster than np.all()
+
+    # SPECIAL CASE 2: 0 in dimension
+    mask2 = (a*b*c).astype(bool)
+
+    # SPECIAL CASE 3: observer lies on-edge/corner
+    # -> 1e-15 to account for numerical inprecision when e.g. rotating
+    # -> /a /b /c to account for the "missing" scaling (1e-15 is large when
+    #    a is e.g. 1e-15 itself)
+
+    mx1 = abs(abs(x)-a) < 1e-15 * a  # on surface
+    my1 = abs(abs(y)-b) < 1e-15 * b  # on surface
+    mz1 = abs(abs(z)-c) < 1e-15 * c  # on surface
+
+    mx2 = (abs(x)-a) < 1e-15 * a # within cuboid dimension
+    my2 = (abs(y)-b) < 1e-15 * b # within cuboid dimension
+    mz2 = (abs(z)-c) < 1e-15 * c # within cuboid dimension
+
+    mask_xedge = my1 & mz1 & mx2
+    mask_yedge = mx1 & mz1 & my2
+    mask_zedge = mx1 & my1 & mz2
+    mask3 = mask_xedge | mask_yedge | mask_zedge
+
+    # on-wall is not a special case
+
+    # continue only with general cases ----------------------------
+    mask_gen = ~mask1 & mask2 & ~mask3
+
+    if np.any(mask_gen):
+        magx, magy, magz = magnetization[mask_gen].T
+        a, b, c = dimension[mask_gen].T/2
+        x, y, z = np.copy(observer[mask_gen]).T
+
+
+        # avoid indeterminate forms by evaluating in bottQ4 only --------
+        # basic masks
+        maskx = x<0
+        masky = y>0
+        maskz = z>0
+
+        # change all positions to their bottQ4 counterparts
+        x[maskx] = x[maskx]*-1
+        y[masky] = y[masky]*-1
+        z[maskz] = z[maskz]*-1
+
+        # create sign flips for position changes
+        qsigns = np.ones((len(magx),3,3))
+        qs_flipx = np.array([[ 1,-1,-1],[-1, 1, 1],[-1, 1, 1]])
+        qs_flipy = np.array([[ 1,-1, 1],[-1, 1,-1],[ 1,-1, 1]])
+        qs_flipz = np.array([[ 1, 1,-1],[ 1, 1,-1],[-1,-1, 1]])
+        # signs flips can be applied subsequently
+        qsigns[maskx] = qsigns[maskx]*qs_flipx
+        qsigns[masky] = qsigns[masky]*qs_flipy
+        qsigns[maskz] = qsigns[maskz]*qs_flipz
+
+        # field computations --------------------------------------------
+        # Note: in principle the computation for all three mag-components can be
+        #   vectorized itself using symmetries. However, tiling the three
+        #   components will cost more than is gained by the vectorized evaluation
+
+        # Note: making the following computation steps is not necessary
+        #   as mkl will cache such small computations
+        xma, xpa = x-a, x+a
+        ymb, ypb = y-b, y+b
+        zmc, zpc = z-c, z+c
+
+        xma2, xpa2 = xma**2, xpa**2
+        ymb2, ypb2 = ymb**2, ypb**2
+        zmc2, zpc2 = zmc**2, zpc**2
+
+        mmm = np.sqrt(xma2 + ymb2 + zmc2)
+        pmp = np.sqrt(xpa2 + ymb2 + zpc2)
+        pmm = np.sqrt(xpa2 + ymb2 + zmc2)
+        mmp = np.sqrt(xma2 + ymb2 + zpc2)
+        mpm = np.sqrt(xma2 + ypb2 + zmc2)
+        ppp = np.sqrt(xpa2 + ypb2 + zpc2)
+        ppm = np.sqrt(xpa2 + ypb2 + zmc2)
+        mpp = np.sqrt(xma2 + ypb2 + zpc2)
+
+        with np.errstate(divide='ignore', invalid='ignore'):
+            ff2x = np.log((xma+mmm) * (xpa+ppm) * (xpa+pmp) * (xma+mpp))     \
+                    -np.log((xpa+pmm) * (xma+mpm) * (xma+mmp) * (xpa+ppp))
+
+            ff2y = np.log((-ymb+mmm) * (-ypb+ppm) * (-ymb+pmp) * (-ypb+mpp)) \
+                -np.log((-ymb+pmm) * (-ypb+mpm) * ( ymb-mmp) * ( ypb-ppp))
+
+            ff2z = np.log((-zmc+mmm) * (-zmc+ppm) * (-zpc+pmp) * (-zpc+mpp)) \
+                -np.log((-zmc+pmm) * ( zmc-mpm) * (-zpc+mmp) * ( zpc-ppp))
+
+        ff1x = (np.arctan2((ymb*zmc),(xma*mmm)) - np.arctan2((ymb*zmc),(xpa*pmm))
+                - np.arctan2((ypb*zmc),(xma*mpm)) + np.arctan2((ypb*zmc),(xpa*ppm))
+                - np.arctan2((ymb*zpc),(xma*mmp)) + np.arctan2((ymb*zpc),(xpa*pmp))
+                + np.arctan2((ypb*zpc),(xma*mpp)) - np.arctan2((ypb*zpc),(xpa*ppp)))
+
+        ff1y = (np.arctan2((xma*zmc),(ymb*mmm)) - np.arctan2((xpa*zmc),(ymb*pmm))
+                - np.arctan2((xma*zmc),(ypb*mpm)) + np.arctan2((xpa*zmc),(ypb*ppm))
+                - np.arctan2((xma*zpc),(ymb*mmp)) + np.arctan2((xpa*zpc),(ymb*pmp))
+                + np.arctan2((xma*zpc),(ypb*mpp)) - np.arctan2((xpa*zpc),(ypb*ppp)))
+
+        ff1z = (np.arctan2((xma*ymb),(zmc*mmm)) - np.arctan2((xpa*ymb),(zmc*pmm))
+                - np.arctan2((xma*ypb),(zmc*mpm)) + np.arctan2((xpa*ypb),(zmc*ppm))
+                - np.arctan2((xma*ymb),(zpc*mmp)) + np.arctan2((xpa*ymb),(zpc*pmp))
+                + np.arctan2((xma*ypb),(zpc*mpp)) - np.arctan2((xpa*ypb),(zpc*ppp)))
+
+        # contributions from x-magnetization
+        bx_magx = magx * ff1x * qsigns[:,0,0]  # the 'missing' third sign is hidden in ff1x
+        by_magx = magx * ff2z * qsigns[:,0,1]
+        bz_magx = magx * ff2y * qsigns[:,0,2]
+        # contributions from y-magnetization
+        bx_magy =  magy * ff2z * qsigns[:,1,0]
+        by_magy =  magy * ff1y * qsigns[:,1,1]
+        bz_magy = -magy * ff2x * qsigns[:,1,2]
+        # contributions from z-magnetization
+        bx_magz =  magz * ff2y * qsigns[:,2,0]
+        by_magz = -magz * ff2x * qsigns[:,2,1]
+        bz_magz =  magz * ff1z * qsigns[:,2,2]
+
+        # summing all contributions
+        bx_tot = bx_magx + bx_magy + bx_magz
+        by_tot = by_magx + by_magy + by_magz
+        bz_tot = bz_magx + bz_magy + bz_magz
+
+        # B = np.c_[bx_tot, by_tot, bz_tot]      # faster for 10^5 and more evaluations
+        B = np.concatenate(((bx_tot,),(by_tot,),(bz_tot,)), axis=0).T
+
+        # combine with special edge/corner cases
+        B_all[mask_gen] = B
+
+    B = B_all / (4*np.pi)
+
+    # return B or compute and return H -------------
+    if bh:
+        return B
+
+    # if inside magnet subtract magnetization vector
+    mask_inside = mx2 & my2 & mz2
+    B[mask_inside] -= magnetization[mask_inside]
+    H = B*10/4/np.pi # mT -> kA/m
+    return H
diff --git a/magpylib/_src/fields/field_BH_cylinder.py b/magpylib/_src/fields/field_BH_cylinder.py
index 93ae99f11..d5816504e 100644
--- a/magpylib/_src/fields/field_BH_cylinder.py
+++ b/magpylib/_src/fields/field_BH_cylinder.py
@@ -7,13 +7,14 @@
 import numpy as np
 from scipy.special import ellipk, ellipe
 from magpylib._src.fields.special_cel import cel
-from magpylib._src.utility import close
+from magpylib._src.utility import cyl_field_to_cart, cart_to_cyl_coordinates
+from magpylib._src.input_checks import check_field_input
 
 
-# ON INTERFACE
-def magnet_cyl_axial_B_Derby2009(
-    dim: np.ndarray,
-    pos_obs: np.ndarray) -> list:
+def fieldB_cylinder_axial(
+    z0: np.ndarray,
+    r: np.ndarray,
+    z: np.ndarray) -> list:
     """
     B-field in Cylindrical CS of Cylinder magnet with homogenous axial unit
     magnetization. The Cylinder axis coincides with the z-axis of the
@@ -32,34 +33,12 @@ def magnet_cyl_axial_B_Derby2009(
     -------
     B-field: ndarray
         B-field array of shape (n,2) in cylindrical coordinates (Br,Bz) in units of [mT].
-
-    Examples
-    --------
-    Compute the field of three instances.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> dim = np.array([(2,2), (2,3), (3,4)])
-    >>> obs = np.array([(.1,2), (2,3), (4,4)])
-    >>> B = magpy.lib.magnet_cyl_axial_B_Derby2009(dim, obs)
-    >>> print(B)
-    [[0.0080244  0.12016141]
-     [0.02622945 0.01635106]
-     [0.01931265 0.00475101]]
-
-    Notes
-    -----
-    Field computed from the current picture (perfect solenoid). Integrals reduced to cel
-    function.
     """
-
-    d,h = dim.T / 2       # d/h are now radius and h/2
-    r,z = pos_obs.T
-    n = len(d)
+    n = len(z0)
 
     # some important quantities
-    zph, zmh = z+h, z-h
-    dpr, dmr = d+r, d-r
+    zph, zmh = z+z0, z-z0
+    dpr, dmr = 1+r, 1-r
 
     sq0 = np.sqrt(zmh**2+dpr**2)
     sq1 = np.sqrt(zph**2+dpr**2)
@@ -70,20 +49,20 @@ def magnet_cyl_axial_B_Derby2009(
     one = np.ones(n)
 
     # radial field (unit magnetization)
-    Br = d*(cel(k1, one, one, -one)/sq1 - cel(k0, one, one, -one)/sq0)/np.pi
+    Br = (cel(k1, one, one, -one)/sq1 - cel(k0, one, one, -one)/sq0)/np.pi
 
     # axial field (unit magnetization)
-    Bz = d/dpr*(zph*cel(k1, gamma**2, one, gamma)/sq1
+    Bz = 1/dpr*(zph*cel(k1, gamma**2, one, gamma)/sq1
               - zmh*cel(k0, gamma**2, one, gamma)/sq0)/np.pi
 
-    return np.array([Br, Bz]).T
+    return Br, Bz
 
 
-# ON INTERFACE
-def magnet_cyl_dia_H_Rauber2021(
-        tetta: np.ndarray,
-        dim: np.ndarray,
-        pos_obs: np.ndarray,
+def fieldH_cylinder_diametral(
+        z0: np.ndarray,
+        r: np.ndarray,
+        phi: np.ndarray,
+        z: np.ndarray,
         ) -> np.ndarray:
     """
     H-field in Cylindrical CS of Cylinder magnet with homogenous
@@ -92,6 +71,11 @@ def magnet_cyl_dia_H_Rauber2021(
 
     Implementation from [Rauber2021].
 
+    H-Field computed analytically via the magnetic scalar potential. Final integration
+    reduced to complete elliptic integrals.
+
+    Numerical Instabilities: See discussion on GitHub.
+
     Parameters
     ----------
     dim: ndarray, shape (n,2)
@@ -105,98 +89,245 @@ def magnet_cyl_dia_H_Rauber2021(
     -------
     H-field: ndarray
         H-field array of shape (n,3) in cylindrical coordinates (Hr, Hphi, Hz) in units of [kA/m].
+    """
+    # pylint: disable=too-many-statements
+
+    n = len(z0)
+
+    # allocate to treat small r special cases
+    Hr, Hphi, Hz = np.empty((3,n))
+
+    # compute repeated quantities for all cases
+    zp = z+z0
+    zm = z-z0
+
+    zp2 = zp**2
+    zm2 = zm**2
+    r2 = r**2
+
+    # case small_r: numerical instability of general solution
+    mask_small_r = r<.05
+    mask_general = ~mask_small_r
+    if np.any(mask_small_r):
+        phiX = phi[mask_small_r]
+        zpX, zmX = zp[mask_small_r], zm[mask_small_r]
+        zp2X, zm2X = zp2[mask_small_r], zm2[mask_small_r]
+        rX, r2X = r[mask_small_r], r2[mask_small_r]
+
+        # taylor series for small r
+        zpp = zp2X + 1
+        zmm = zm2X + 1
+        sqrt_p = np.sqrt(zpp)
+        sqrt_m = np.sqrt(zmm)
+
+        frac1 = zpX/sqrt_p
+        frac2 = zmX/sqrt_m
+
+        r3X = r2X*rX
+        r4X = r3X*rX
+        r5X = r4X*rX
+
+        term1 =  (frac1 - frac2)
+        term2 = (frac1/zpp**2 - frac2/zmm**2)*r2X/8
+        term3 = ((3-4*zp2X)*frac1/zpp**4 - (3-4*zm2X)*frac2/zmm**4)/64*r4X
+
+        Hr[mask_small_r] = -np.cos(phiX)/4*( term1 + 9*term2 + 25*term3)
+
+        Hphi[mask_small_r] = np.sin(phiX)/4*( term1 + 3*term2 + 5*term3)
+
+        Hz[mask_small_r] = -np.cos(phiX)/4*(
+            rX*(1/zpp/sqrt_p - 1/zmm/sqrt_m) +
+            3/8*r3X*((1 - 4*zp2X)/zpp**3/sqrt_p - (1 - 4*zm2X)/zmm**3/sqrt_m) +
+            15/64*r5X*((1-12*zp2X+8*zp2X**2)/zpp**5/sqrt_p - (1-12*zm2X+8*zm2X**2)/zmm**5/sqrt_m))
+
+        # if there are small_r, select the general/case variables
+        # when there are no small_r cases it is not necessary to slice with [True, True, Tue,...]
+        phi = phi[mask_general]
+        n = len(phi)
+        zp, zm = zp[mask_general], zm[mask_general]
+        zp2, zm2 = zp2[mask_general], zm2[mask_general]
+        r, r2 = r[mask_general], r2[mask_general]
+
+
+    if np.any(mask_general):
+        rp = r+1
+        rm = r-1
+        rp2 = rp**2
+        rm2 = rm**2
+
+        ap2 = zp2+rm**2
+        am2 = zm2+rm**2
+        ap = np.sqrt(ap2)
+        am = np.sqrt(am2)
+
+        argp = -4*r/ap2
+        argm = -4*r/am2
+
+        # special case r=r0 : indefinite form
+        #   result is numerically stable in the vicinity of of r=r0
+        #   so only the special case must be caught (not the surroundings)
+        mask_special = rm==0
+        argc = np.ones(n)*1e16      # should be np.Inf but leads to 1/0 problems in cel
+        argc[~mask_special] = -4*r[~mask_special]/rm2[~mask_special]
+        # special case 1/rm
+        one_over_rm = np.zeros(n)
+        one_over_rm[~mask_special] = 1/rm[~mask_special]
+
+        elle_p = ellipe(argp)
+        elle_m = ellipe(argm)
+        ellk_p = ellipk(argp)
+        ellk_m = ellipk(argm)
+        onez = np.ones(n)
+        ellpi_p = cel(np.sqrt(1-argp), 1-argc, onez, onez) # elliptic_Pi
+        ellpi_m = cel(np.sqrt(1-argm), 1-argc, onez, onez) # elliptic_Pi
+
+        # compute fields
+        Hr[mask_general] = - np.cos(phi)/(4*np.pi*r2)*(
+            - zm*am             * elle_m   +  zp*ap             * elle_p
+            + zm/am*(2+zm2) * ellk_m   -  zp/ap*(2+zp2) * ellk_p
+            + (zm/am* ellpi_m  -  zp/ap * ellpi_p) * rp*(r2+1) * one_over_rm)
+
+        Hphi[mask_general] = np.sin(phi)/(4*np.pi*r2)*(
+            + zm*am              * elle_m   -  zp*ap              * elle_p
+            - zm/am*(2+zm2+2*r2) * ellk_m   +  zp/ap*(2+zp2+2*r2) * ellk_p
+            + zm/am*rp2          * ellpi_m  -  zp/ap*rp2          * ellpi_p
+            )
+
+        Hz[mask_general] = - np.cos(phi)/(2*np.pi*r)*(
+            + am              * elle_m  -  ap              * elle_p
+            - (1+zm2+r2)/am * ellk_m  +  (1+zp2+r2)/ap * ellk_p
+            )
+
+    return Hr, Hphi, Hz
+
+
+# ON INTERFACE
+def magnet_cylinder_field(
+    magnetization: np.ndarray,
+    dimension: np.ndarray,
+    observer: np.ndarray,
+    field='B',
+    ) -> np.ndarray:
+    """
+    Computes the magnetic field of a homogeneously magnetized cylinder
+    in Cartesian coordinates. The cylinder axis coincides with the z-axis and
+    the geometric center of the cylinder lies in the origin.
+
+    Parameters
+    ----------
+    magnetization: ndarray, shape (n,3)
+        Homogeneous magnetization vector in units of [mT].
+
+    dimension: ndarray, shape (n,2)
+        Cylinder dimension (d,h) with diameter d and height h in units of [mm].
+
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
+
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of magnet in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
 
     Examples
     --------
-    Compute field for three instances.
+    Compute the B-field of two different cylinder magnets at position (1,2,3).
 
     >>> import numpy as np
     >>> import magpylib as magpy
-    >>> tetta = np.zeros(3)
-    >>> dim = np.array([(2,2), (2,3), (3,4)])
-    >>> obs = np.array([(.1,0,2), (2,0.12,3), (4,0.2,1)])
-    >>> B = magpy.lib.magnet_cyl_dia_H_Rauber2021(tetta, dim, obs)
+    >>> mag = np.array([(0,0,1000), (100,0,100)])
+    >>> dim = np.array([(1,1), (2,3)])
+    >>> obs = np.array([(1,2,3), (1,2,3)])
+    >>> B = magpy.lib.magnet_cylinder_field(mag, dim, obs)
     >>> print(B)
-    [[-0.05992403  0.          0.0080244 ]
-     [ 0.00193283  0.00219048  0.02604082]
-     [ 0.05270086  0.00606112  0.01546927]]
+    [[ 0.77141782  1.54283565  1.10384481]
+     [-0.15185713  2.90352915  2.23601722]]
 
     Notes
     -----
-    H-Field computed analytically via the magnetic scalar potential. Final integration
-    reduced to complete elliptic integrals.
-
-    Numerical Instabilities: See discussion on GitHub.
+    Axial implementation based on [Derby].
+    Diametral implementation based on [caciagli2018] and [rauber/leitner/ortner 2022 wip].
     """
 
-    r0, z0 = dim.T/2
-    r, phi, z = pos_obs.T
-    n = len(r0)
+    bh = check_field_input(field, 'magnet_cylinder_field()')
 
-    # phi is now relative between mag and pos_obs
-    phi = phi-tetta
+    # transform to Cy CS --------------------------------------------
+    r, phi, z = cart_to_cyl_coordinates(observer)
+    r0,z0 = dimension.T/2
 
-    # implementation of Rauber2021
-    # compute repeated quantities
-    zp = z+z0
-    zm = z-z0
-    rp = r+r0
-    rm = r-r0
+    # scale invariance (make dimensionless)
+    r  = np.copy(r/r0)
+    z  = np.copy(z/r0)
+    z0 = np.copy(z0/r0)
 
-    zp2 = zp**2
-    zm2 = zm**2
-    rp2 = rp**2
-    rm2 = rm**2
-    r2 = r**2
-    r02 = r0**2
-    rr0 = r*r0
-
-    ap2 = zp2+rm**2
-    am2 = zm2+rm**2
-    ap = np.sqrt(ap2)
-    am = np.sqrt(am2)
-
-    argp = -4*rr0/ap2
-    argm = -4*rr0/am2
-
-    # special case r=r0 : indefinite form
-    #   result is numerically stable in the vicinity of of r=r0
-    #   so only the special case must be caught (not the surroundings)
-    mask_special = rm==0
-    argc = np.ones(n)*1e16      # should be np.Inf but leads to 1/0 problems in cel
-    argc[~mask_special] = -4*rr0[~mask_special]/rm2[~mask_special]
-
-    elle_p = ellipe(argp)
-    elle_m = ellipe(argm)
-    ellk_p = ellipk(argp)
-    ellk_m = ellipk(argm)
-    onez = np.ones(n)
-    ellpi_p = cel(np.sqrt(1-argp), 1-argc, onez, onez) # elliptic_Pi
-    ellpi_m = cel(np.sqrt(1-argm), 1-argc, onez, onez) # elliptic_Pi
-
-    # compute fields
-    Hphi = np.sin(phi)/(4*np.pi*r2)*(
-        + zm*am                  * elle_m   -  zp*ap                  * elle_p
-        - zm/am*(2*r02+zm2+2*r2) * ellk_m   +  zp/ap*(2*r02+zp2+2*r2) * ellk_p
-        + zm/am*rp2              * ellpi_m  -  zp/ap*rp2              * ellpi_p
-        )
-
-    Hz = - np.cos(phi)/(2*np.pi*r)*(
-        + am              * elle_m  -  ap              * elle_p
-        - (r02+zm2+r2)/am * ellk_m  +  (r02+zp2+r2)/ap * ellk_p
-        )
-
-    # special case 1/rm
-    one_over_rm = np.zeros(n)
-    one_over_rm[~mask_special] = 1/rm[~mask_special]
-
-    Hr = - np.cos(phi)/(4*np.pi*r2)*(
-        - zm*am             * elle_m   +  zp*ap             * elle_p
-        + zm/am*(2*r02+zm2) * ellk_m   -  zp/ap*(2*r02+zp2) * ellk_p
-        + (zm/am* ellpi_m  -  zp/ap * ellpi_p) * rp*(r2+r02) * one_over_rm)
-
-    return np.array([Hr, Hphi, Hz]).T
+    # allocate field vectors ----------------------------------------
+    Br, Bphi, Bz = np.zeros((3,len(r)))
+
+    # create masks to distinguish between cases ---------------------
+    m0 = np.isclose(r, 1, rtol=1e-15, atol=0)         # on Cylinder hull plane
+    m1 = np.isclose(abs(z), z0, rtol=1e-15, atol=0)   # on top or bottom plane
+    m2 = np.abs(z)<=z0      # in-between top and bottom plane
+    m3 = r<=1               # inside Cylinder hull plane
+
+    # special case: mag = 0
+    mask0 = (np.linalg.norm(magnetization, axis=1)==0)
+
+    # special case: on Cylinder edge
+    mask_edge = (m0 & m1)
 
+    # general case
+    mask_gen = ~mask0 & ~mask_edge
+
+    # axial/transv magnetization cases
+    magx, magy, magz = magnetization.T
+    mask_tv = (magx != 0) | (magy != 0)
+    mask_ax = (magz != 0)
+
+    # inside/outside
+    mask_inside = m2 & m3
+
+    # general case masks
+    mask_tv = mask_tv & mask_gen
+    mask_ax = mask_ax & mask_gen
+    mask_inside = mask_inside & mask_gen
+
+    # transversal magnetization contributions -----------------------
+    if any(mask_tv):
+        magxy = np.sqrt(magx**2 + magy**2)[mask_tv]
+        tetta = np.arctan2(magy[mask_tv], magx[mask_tv])
+        br_tv, bphi_tv, bz_tv = fieldH_cylinder_diametral(z0[mask_tv], r[mask_tv],
+            phi[mask_tv]-tetta, z[mask_tv])
+
+        # add to H-field (inside magxy is missing for B !!!)
+        Br[mask_tv]   += magxy*br_tv
+        Bphi[mask_tv] += magxy*bphi_tv
+        Bz[mask_tv]   += magxy*bz_tv
+
+    # axial magnetization contributions -----------------------------
+    if any(mask_ax):
+        br_ax, bz_ax = fieldB_cylinder_axial(z0[mask_ax], r[mask_ax], z[mask_ax])
+        Br[mask_ax] += magz[mask_ax]*br_ax
+        Bz[mask_ax] += magz[mask_ax]*bz_ax
+
+    # transform field to cartesian CS -------------------------------
+    Bx, By = cyl_field_to_cart(phi, Br, Bphi)
+
+    # add/subtract Mag when inside for B/H --------------------------
+    if bh:
+        if any(mask_tv): # tv computes H-field
+            Bx[mask_tv*mask_inside] += magx[mask_tv*mask_inside]
+            By[mask_tv*mask_inside] += magy[mask_tv*mask_inside]
+        return np.concatenate(((Bx,),(By,),(Bz,)),axis=0).T
+
+    if any(mask_ax): # ax computes B-field
+        Bz[mask_tv*mask_inside] -= magz[mask_tv*mask_inside]
+    return np.concatenate(((Bx,),(By,),(Bz,)),axis=0).T*10/4/np.pi
+
+
+# old iterative solution by Furlani
 
 # def magnet_cyl_dia_H_Furlani1994(
 #         tetta: np.ndarray,
@@ -298,109 +429,3 @@ def magnet_cyl_dia_H_Rauber2021(
 #     Bz = r0/2/niter*np.sum(summandz, axis=0)
 
 #     return np.array([Br, Bphi, Bz]).T
-
-
-
-def field_BH_cylinder(
-        bh: bool,
-        mag: np.ndarray,
-        dim: np.ndarray,
-        pos_obs: np.ndarray
-        ) -> np.ndarray:
-    """ setting up the Cylinder field computation
-    - transform to Cylindrical CS
-    - separate mag=0 cases (returning 0)
-    - separate edge/corner cases (returning 0)
-    - separate magz and magxy
-    - call field computation for general cases
-    - select B or H
-    - transform B<-->H (inside check)
-
-    ### Args:
-    - bh (boolean): True=B, False=H
-    - mag (ndarray Nx3): homogeneous magnetization vector in units of mT
-    - dim (ndarray Nx2): dimension of Cylinder DxH in units of mm
-    - pos_obs (ndarray Nx3): position of observer in units of mm
-
-    ### Returns:
-    - B/H-field (ndarray Nx3): magnetic field vectors at pos_obs in units of mT / kA/m
-    """
-
-    # transform to Cy CS --------------------------------------------
-    x, y, z = pos_obs.T
-    r, phi = np.sqrt(x**2+y**2), np.arctan2(y, x)
-    pos_obs_cy = np.concatenate(((r,),(phi,),(z,)),axis=0).T
-
-    # allocate field vectors ----------------------------------------
-    Br, Bphi, Bz = np.zeros((3,len(x)))
-
-    # create masks to distinguish between cases ---------------------
-
-    r0,z0 = dim.T/2
-    m0 = close(r, r0)        # on Cylinder plane
-    m1 = close(abs(z), z0)   # on top or bottom plane
-    m2 = np.abs(z)<=z0       # inside Cylinder plane
-    m3 = r<=r0               # in-between top and bottom plane
-
-    # special case: mag = 0
-    mask0 = (np.linalg.norm(mag,axis=1)==0)
-
-    # special case: on Cylinder surface
-    mask_edge = (m0 & m2) | (m1 & m3)
-
-    # general case
-    mask_gen = ~mask0 & ~mask_edge
-
-    # axial/transv magnetization cases
-    magx, magy, magz = mag.T
-    mask_tv = (magx != 0) | (magy != 0)
-    mask_ax = (magz != 0)
-
-    # inside/outside
-    mask_inside = m2 & m3
-
-    # general case masks
-    mask_tv = mask_tv & mask_gen
-    mask_ax = mask_ax & mask_gen
-    mask_inside = mask_inside & mask_gen
-
-    # transversal magnetization contributions -----------------------
-    if any(mask_tv):
-        # select non-zero tv parts
-        magxy = np.sqrt(magx**2 + magy**2)[mask_tv]
-        tetta = np.arctan2(magy[mask_tv], magx[mask_tv])
-        pos_obs_tv = pos_obs_cy[mask_tv]
-        dim_tv = dim[mask_tv]
-        # compute H-field (in mT)
-        br_tv, bphi_tv, bz_tv = magnet_cyl_dia_H_Rauber2021(tetta, dim_tv, pos_obs_tv).T
-        # add to H-field (inside magxy is missing for B)
-        Br[mask_tv]   += magxy*br_tv
-        Bphi[mask_tv] += magxy*bphi_tv
-        Bz[mask_tv]   += magxy*bz_tv
-
-    # axial magnetization contributions -----------------------------
-    if any(mask_ax):
-        # select non-zero ax parts
-        pos_obs_ax = pos_obs_cy[mask_ax]
-        magz_ax = magz[mask_ax]
-        dim_ax = dim[mask_ax]
-        # compute B-field
-        br_ax, bz_ax = magnet_cyl_axial_B_Derby2009(dim_ax, pos_obs_ax[:,(0,2)]).T
-        # add to B-field
-        Br[mask_ax] += magz_ax*br_ax
-        Bz[mask_ax] += magz_ax*bz_ax
-
-    # transform field to cartesian CS -------------------------------
-    Bx = Br*np.cos(phi) - Bphi*np.sin(phi)
-    By = Br*np.sin(phi) + Bphi*np.cos(phi)
-
-    # add/subtract Mag when inside for B/H --------------------------
-    if bh:
-        if any(mask_tv): # tv computes H-field
-            Bx[mask_tv*mask_inside] += magx[mask_tv*mask_inside]
-            By[mask_tv*mask_inside] += magy[mask_tv*mask_inside]
-        return np.concatenate(((Bx,),(By,),(Bz,)),axis=0).T
-
-    if any(mask_ax): # ax computes B-field
-        Bz[mask_tv*mask_inside] -= magz[mask_tv*mask_inside]
-    return np.concatenate(((Bx,),(By,),(Bz,)),axis=0).T*10/4/np.pi
diff --git a/magpylib/_src/fields/field_BH_cylinder_tile.py b/magpylib/_src/fields/field_BH_cylinder_segment.py
similarity index 95%
rename from magpylib/_src/fields/field_BH_cylinder_tile.py
rename to magpylib/_src/fields/field_BH_cylinder_segment.py
index 452eb4f94..b0f0831a7 100644
--- a/magpylib/_src/fields/field_BH_cylinder_tile.py
+++ b/magpylib/_src/fields/field_BH_cylinder_segment.py
@@ -5,9 +5,7 @@
 import numpy as np
 from scipy.special import ellipeinc, ellipkinc
 from magpylib._src.fields.special_el3 import el3_angle
-from magpylib._src.utility import close
-from magpylib._src.defaults.defaults_classes import default_settings as Config
-
+from magpylib._src.input_checks import check_field_input
 
 def arctan_k_tan_2(k, phi):
     """
@@ -26,6 +24,15 @@ def arctan_k_tan_2(k, phi):
     return np.where(np.abs(phi_red) < np.pi, result + np.arctan(k * np.tan(phi_red / 2.0)), result + phi_red / 2.0)
 
 
+def close(arg1: np.ndarray, arg2: np.ndarray) -> np.ndarray:
+    """
+    determine if arg1 and arg2 lie close to each other
+    input: ndarray, shape (n,) or numpy-interpretable scalar
+    output: ndarray, dtype=bool
+    """
+    return np.isclose(arg1, arg2, rtol=1e-14, atol=1e-14)
+
+
 def determine_cases(r, phi, z, r1, phi1, z1):
     """
     Determine case of input parameter set.
@@ -1243,13 +1250,13 @@ def case235(r, r_i, r_bar_i, phi_bar_j, phi_bar_M, phi_bar_Mj, theta_M, z_bar_k)
     return results
 
 # ON INTERFACE
-def magnet_cyl_tile_H_Slanovc2021(
+def magnet_cylinder_segment_core(
     mag: np.ndarray,
     dim: np.ndarray,
     obs_pos: np.ndarray
     ) ->np.ndarray:
     """
-    H-field of Cylinder Tile magnet with homogenous magnetization.
+    H-field of Cylinder segment magnet with homogenous magnetization.
     The full cylinder axis coincides with the z-axis of the CS.
     The geometric center of the full Cylinder is in the origin.
 
@@ -1262,7 +1269,7 @@ def magnet_cyl_tile_H_Slanovc2021(
     obs_pos : ndarray, shape (n,3)
         observer positions (r,phi,z) in cy CS, units: [mm] [rad]
     dim: ndarray, shape (n,6)
-        section dimensions (r1,r2,phi1,phi2,z1,z2) in cy CS , units: [mm] [rad]
+        segment dimensions (r1,r2,phi1,phi2,z1,z2) in cy CS , units: [mm] [rad]
 
     Returns
     -------
@@ -1279,7 +1286,7 @@ def magnet_cyl_tile_H_Slanovc2021(
     >>> mag = np.array([(100,0,0), (200,.1,pi/4)])
     >>> dim = np.array([(1,2,0,pi/2,-1,1), (.1,3,-.3,pi,0,1)])
     >>> obs = np.array([(.1,0,3), (1,pi,3)])
-    >>> B = magpy.lib.magnet_cyl_tile_H_Slanovc2021(mag, dim, obs)
+    >>> B = magpy.lib.magnet_cylinder_segment_core(mag, dim, obs)
     >>> print(B)
     [[-0.84506541 -0.9207606   1.48474874]
      [ 3.95719801  3.59131966  3.11703698]]
@@ -1342,35 +1349,72 @@ def magnet_cyl_tile_H_Slanovc2021(
     return result.T
 
 
-def field_BH_cylinder_tile(
-        bh: bool,
-        mag: np.ndarray,
-        dim: np.ndarray,
-        pos_obs: np.ndarray
-        ) -> np.ndarray:
+def magnet_cylinder_segment_field(
+    magnetization: np.ndarray,
+    dimension: np.ndarray,
+    observer: np.ndarray,
+    field='B',
+    ) -> np.ndarray:
     """
-    ### Args:
-    - bh (boolean): True=B, False=H
-    - mag (ndarray Nx3): homogeneous magnetization vector in cartesian CS units of [mT]
-    - dim (ndarray Nx6): dimension of Cylinder (r1,r2,phi1,phi2,z1,z2) in units of [mm], [deg]
-        expects 0<=r1<r2, -2pi<phi1<phi2<2pi, z1<z2
-    - pos_obs (ndarray Nx3): position of observer in units of [mm]
-    - return 0 when on surface (analytic code yields unphysical solutions on edges and surfaces)
-
-    ### Returns:
-    - B/H-field (ndarray Nx3): magnetic field vectors in cartesian CS at pos_obs
-        in units of [mT] or [kA/m]
+    Computes the magnetic field of a homogeneously magnetized cylinder segment in
+    Cartesian coordinates. The full cylinder axis coincides with the z-axis of the coordinate
+    system. The geometric center of the full cylinder is in the origin.
+
+    Parameters
+    ----------
+    magnetization: ndarray, shape (n,3)
+        Homogeneous magnetization vector in units of [mT].
+
+    dimension: ndarray, shape (n,5)
+        Cylinder segment dimensions (r1,r2,h,phi1,phi2) with inner radius r1, outer radius r2,
+        height h in units of [mm] and the two segment angles phi1 and phi2 in units of [deg].
+
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
+
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of magnet in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
+
+    Examples
+    --------
+    Compute the field of two different cylinder segment magnets at position (1,1,1).
+
+    >>> import numpy as np
+    >>> import magpylib as magpy
+    >>> mag = np.array([(0,0,100), (50,50,0)])
+    >>> dim = np.array([(0,1,2,0,90), (1,2,4,35,125)])
+    >>> obs = np.array([(1,1,1), (1,1,1)])
+    >>> B = magpy.lib.magnet_cylinder_segment_field(mag, dim, obs)
+    >>> print(B)
+    [[ 6.27410168  6.27410168 -1.20044166]
+     [29.84602335 20.75731598  0.34961733]]
+
+    Notes
+    -----
+    Implementation based on [Slanovc2022].
     """
-    BHfinal = np.zeros((len(mag),3))
+    bh = check_field_input(field, 'magnet_cylinder_segment_field()')
+
+    BHfinal = np.zeros((len(magnetization),3))
+
+    r1,r2,h,phi1,phi2 = dimension.T
+    r1 = abs(r1)
+    r2 = abs(r2)
+    h = abs(h)
+    z1, z2 = -h/2, h/2
 
     # transform dim deg->rad
-    r1, r2, phi1, phi2, z1, z2 = dim.T
     phi1 = phi1/180*np.pi
     phi2 = phi2/180*np.pi
     dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
 
     # transform obs_pos to Cy CS --------------------------------------------
-    x, y, z = pos_obs.T
+    x, y, z = observer.T
     r, phi = np.sqrt(x**2+y**2), np.arctan2(y, x)
     pos_obs_cy = np.concatenate(((r,),(phi,),(z,)),axis=0).T
 
@@ -1384,10 +1428,10 @@ def field_BH_cylinder_tile(
     mask_phi1 = close(phio1, phi1) | close(phio2, phi1)
     mask_phi2 = close(phio1, phi2) | close(phio2, phi2)
 
-    # r, phi ,z lies in-between
-    mask_r_in = (r1-Config.edgesize<r) & (r<r2+Config.edgesize)
+    # r, phi ,z lies in-between, avoid numerical fluctuations (e.g. due to rotations) by including 1e-14
+    mask_r_in = (r1-1e-14<r) & (r<r2+1e-14)
     mask_phi_in = (np.sign(phio1-phi1)!=np.sign(phio1-phi2)) | (np.sign(phio2-phi1)!=np.sign(phio2-phi2))
-    mask_z_in = (z1-Config.edgesize<z) & (z<z2+Config.edgesize)
+    mask_z_in = (z1-1e-14<z) & (z<z2+1e-14)
 
     # on surface
     mask_surf_z = (close(z, z1) | close(z, z2)) & mask_phi_in & mask_r_in # top / bottom
@@ -1404,7 +1448,7 @@ def field_BH_cylinder_tile(
         return BHfinal
 
     # redefine input if there are some surface-points -------------------------
-    magg = mag[mask_not_on_surf]
+    magg = magnetization[mask_not_on_surf]
     dim = dim[mask_not_on_surf]
     pos_obs_cy = pos_obs_cy[mask_not_on_surf]
     phi = phi[mask_not_on_surf]
@@ -1416,7 +1460,7 @@ def field_BH_cylinder_tile(
     mag_sph = np.concatenate(((m,),(phi_m,),(th_m,)),axis=0).T
 
     # compute H and transform to cart CS -------------------------------------
-    H_cy = magnet_cyl_tile_H_Slanovc2021(mag_sph, dim, pos_obs_cy)
+    H_cy = magnet_cylinder_segment_core(mag_sph, dim, pos_obs_cy)
     Hr, Hphi, Hz = H_cy.T
     Hx = Hr*np.cos(phi) - Hphi*np.sin(phi)
     Hy = Hr*np.sin(phi) + Hphi*np.cos(phi)
@@ -1430,5 +1474,5 @@ def field_BH_cylinder_tile(
     B = H/(10/4/np.pi) # kA/m -> mT
     BHfinal[mask_not_on_surf] = B
     maskX = mask_inside*mask_not_on_surf
-    BHfinal[maskX] += mag[maskX]
+    BHfinal[maskX] += magnetization[maskX]
     return BHfinal
diff --git a/magpylib/_src/fields/field_BH_dipole.py b/magpylib/_src/fields/field_BH_dipole.py
index be0029605..533af1396 100644
--- a/magpylib/_src/fields/field_BH_dipole.py
+++ b/magpylib/_src/fields/field_BH_dipole.py
@@ -3,36 +3,66 @@
 """
 
 import numpy as np
+from magpylib._src.input_checks import check_field_input
 
-
-def field_BH_dipole(
-        bh: bool,
-        moment: np.ndarray,
-        pos_obs: np.ndarray
-        ) -> np.ndarray:
+def dipole_field(
+    moment: np.ndarray,
+    observer: np.ndarray,
+    field='B'
+    ) -> np.ndarray:
     """
-    The B-field of a magnetic dipole moment. The solution is set up so that
-    the moment is given in mT*mm^3, i.e. the field is similar to the outside-field
-    of a spherical magnet with Volume = 1 mm^3.
+    Computes the magnetic field of a magnetic dipole moment in Cartesian coordinates.
 
-    Paramters:
+    Parameters
     ----------
-    - bh (boolean): True=B, False=H
-    - moment (ndarray Nx3): dipole moment vector in units of mT*mm^3
-    - pos_obs (ndarray Nx3): position of observer in units of mm
+    moment: ndarray, shape (n,3)
+        Dipole moment vector in units of [mT*mm^3].
+
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
 
-    Returns:
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of dipole in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
+
+    Examples
     --------
-    B/H-field (ndarray Nx3): magnetic field vectors at pos_obs in units of mT / kA/m
+    >>> import numpy as np
+    >>> import magpylib as magpy
+    >>> mom = np.array([(1,2,3), (0,0,1)])
+    >>> obs = np.array([(1,1,1), (0,0,2)])
+    >>> B = magpy.lib.dipole_Bfield(mom, obs)
+    >>> print(B)
+    [[0.07657346 0.06125877 0.04594407]
+     [0.         0.         0.01989437]]
+
+    Notes
+    -----
+    The field is similar to the outside-field of a spherical magnet with Volume = 1 [mm^3].
     """
+    bh = check_field_input(field, 'dipole_field()')
 
-    x, y, z = pos_obs.T
+    x, y, z = observer.T
     r = np.sqrt(x**2+y**2+z**2)   # faster than np.linalg.norm
-    B = (3*(np.sum(moment*pos_obs,axis=1)*pos_obs.T)/r**5 - moment.T/r**3).T/4/np.pi
+    with np.errstate(divide='ignore', invalid='ignore'):
+        # 0/0 produces invalid warn and results in np.nan
+        # x/0 produces divide warn and results in np.inf
+        B = (3*np.sum(moment*observer,axis=1)*observer.T/r**5 - moment.T/r**3).T/4/np.pi
+
+    # when r=0 return np.inf in all non-zero moment directions
+    mask1 = r==0
+    if np.any(mask1):
+        with np.errstate(divide='ignore', invalid='ignore'):
+            B[mask1] = moment[mask1]/0.
+            np.nan_to_num(B, copy=False, posinf=np.inf, neginf=np.NINF)
 
+    # return B or H
     if bh:
         return B
 
-    # adjust and return H
     H = B*10/4/np.pi
     return H
diff --git a/magpylib/_src/fields/field_BH_line.py b/magpylib/_src/fields/field_BH_line.py
index f73221cb3..909b0187f 100644
--- a/magpylib/_src/fields/field_BH_line.py
+++ b/magpylib/_src/fields/field_BH_line.py
@@ -4,14 +4,14 @@
 
 import numpy as np
 from numpy.linalg import norm
-from magpylib._src.defaults.defaults_classes import default_settings as Config
+from magpylib._src.input_checks import check_field_input
 
 
 def field_BH_line_from_vert(
-    bh: bool,
     current: np.ndarray,
     vertex_sets: list,  # list of mix3 ndarrays
-    pos_obs: np.ndarray
+    pos_obs: np.ndarray,
+    field: str,
     ) -> np.ndarray:
     """
     This function accepts n (mi,3) shaped vertex-sets, creates a single long
@@ -48,7 +48,7 @@ def field_BH_line_from_vert(
     pos_end = np.repeat(pos_end, npp, axis=0)
 
     # compute field
-    field = field_BH_line(bh, curr_tile, pos_start, pos_end, pos_obs)
+    field = current_line_field(curr_tile, pos_start, pos_end, pos_obs, field=field)
     field = np.reshape(field, (nseg, npp, 3))
 
     # sum for each vertex set
@@ -58,76 +58,108 @@ def field_BH_line_from_vert(
     return np.reshape(field_sum, (-1,3))
 
 
-def field_BH_line(
-    bh: bool,
+# ON INTERFACE
+def current_line_field(
     current: np.ndarray,
-    pos_start: np.ndarray,
-    pos_end: np.ndarray,
-    pos_obs: np.ndarray
+    start: np.ndarray,
+    end: np.ndarray,
+    observer: np.ndarray,
+    field='B'
     ) -> np.ndarray:
     """
-    ### Args:
-    - bh (boolean): True=B, False=H
-    - current (float): current on line in units of [A]
-    - pos_start (ndarray nx3) start position of line segments
-    - pos_end (ndarray nx3) end positions of line segments
-    - pos_obs (ndarray nx3): n observer positions in units of [mm]
-
-    ### Returns:
-    - B-field (ndarray nx3): B-field vectors at pos_obs in units of mT
-
-    ### init_state:
-    Line current flowing in a straight line from pos_start to pos_end.
-
-    ### Computation info:
+    Computes the magnetic field of line current segments in Cartesian coordinates.
+
+    The field set to (0,0,0) on line.
+
+    Parameters
+    ----------
+    current: ndarray, shape (n,)
+        Electrical current in units of [A].
+
+    start: ndarray, shape (n,3)
+        Line start positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    end: ndarray, shape (n,3)
+        Line end positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
+
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of current in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
+
+    Examples
+    --------
+    Compute the field of two segments. The 2nd observer lies on the segment
+    so that [0  0  0] is returned.
+
+    >>> import numpy as np
+    >>> import magpylib as magpy
+    >>> curr = np.array([1,2])
+    >>> start = np.array([(-1,0,0), (-1,0,0)])
+    >>> end   = np.array([( 1,0,0), ( 2,0,0)])
+    >>> obs   = np.array([( 0,0,1), ( 0,0,0)])
+    >>> B = magpy.lib.current_line_Bfield(curr, start, end, obs)
+    >>> print(B)
+    [[ 0.         -0.14142136  0.        ]
+     [ 0.          0.          0.        ]]
+
+    Notes
+    -----
     Field computation via law of Biot Savart. See also countless online ressources.
     eg. http://www.phys.uri.edu/gerhard/PHY204/tsl216.pdf
-
-    ### Numerical instabilities:
-        - singularity at r=0, B set to 0 within Config.edgesize
     """
     # pylint: disable=too-many-statements
 
+    bh = check_field_input(field, 'current_line_field()')
+
+    # allocate for special case treatment
+    ntot = len(current)
+    field_all = np.zeros((ntot,3))
+
     # Check for zero-length segments
-    mask0 = np.all(pos_start==pos_end, axis=1)
+    mask0 = np.all(start==end, axis=1)
     if np.all(mask0):
-        n0 = len(pos_obs)
-        return np.zeros((n0,3))
-
-    any_zero_segments = np.any(mask0)
+        return field_all
 
     # continue only with non-zero segments
-    if any_zero_segments:
+    if np.any(mask0):
         not_mask0 = ~mask0     # avoid multiple computation of ~mask
-        p1 = pos_start[not_mask0]
-        p2 = pos_end[not_mask0]
-        po = pos_obs[not_mask0]
         current = current[not_mask0]
-    else:
-        p1 = pos_start        # just renaming
-        p2 = pos_end
-        po = pos_obs
+        start = start[not_mask0]
+        end = end[not_mask0]
+        observer = observer[not_mask0]
+
+    # rename
+    p1,p2,po = start, end, observer
+
+    # make dimensionless (avoid all large/small input problems) by introducing
+    # the segment length as characteristic length scale.
+    norm_12 = norm(p1-p2, axis=1)
+    p1 = (p1.T/norm_12).T
+    p2 = (p2.T/norm_12).T
+    po = (po.T/norm_12).T
 
     # p4 = projection of pos_obs onto line p1-p2
-    p1p2 = p1-p2
-    norm_12 = norm(p1p2, axis=1)
-    t = np.sum((po-p1)*p1p2, axis=1) / norm_12**2
-    p4 = p1 + (t*p1p2.T).T
+    t = np.sum((po-p1)*(p1-p2), axis=1)
+    p4 = p1 + (t*(p1-p2).T).T
 
     # distance of observer from line
     norm_o4 = norm(po - p4, axis=1)
 
-    # on-line cases (set B=0)
-    mask1 = norm_o4 < Config.edgesize
+    # separate on-line cases (-> B=0)
+    mask1 = norm_o4<1e-15 # account for numerical issues
     if np.all(mask1):
-        n0 = len(pos_obs)
-        return np.zeros((n0,3))
+        return field_all
 
-    any_on_line_cases = np.any(mask1)
-
-    # redefine to avoid tons of slices
-    if any_on_line_cases:
-        not_mask1 = ~mask1     # avoid multiple computation of ~mask
+    # continue only with general off-line cases
+    if np.any(mask1):
+        not_mask1 = ~mask1
         po = po[not_mask1]
         p1 = p1[not_mask1]
         p2 = p2[not_mask1]
@@ -142,7 +174,7 @@ def field_BH_line(
     eB = (cros.T/norm_cros).T
 
     # compute angles
-    norm_o1 = norm(po-p1, axis=1)
+    norm_o1 = norm(po-p1, axis=1)   # improve performance by computing all norms at once
     norm_o2 = norm(po-p2, axis=1)
     norm_41 = norm(p4-p1, axis=1)
     norm_42 = norm(p4-p2, axis=1)
@@ -152,33 +184,24 @@ def field_BH_line(
 
     # determine how p1,p2,p4 are sorted on the line (to get sinTH signs)
     # both points below
-    mask2 = ((norm_41>norm_12) * (norm_41>norm_42))
+    mask2 = ((norm_41>1) * (norm_41>norm_42))
     deltaSin[mask2] = abs(sinTh1[mask2]-sinTh2[mask2])
     # both points above
-    mask3 = ((norm_42>norm_12) * (norm_42>norm_41))
+    mask3 = ((norm_42>1) * (norm_42>norm_41))
     deltaSin[mask3] = abs(sinTh2[mask3]-sinTh1[mask3])
     # one above one below or one equals p4
     mask4 = ~mask2 * ~mask3
     deltaSin[mask4] = abs(sinTh1[mask4]+sinTh2[mask4])
 
-    field = (deltaSin/norm_o4*eB.T* current/10).T # m->mm, T->mT
-
-    if any_on_line_cases: # broadcast into np.zeros
-        n1 = len(mask1)
-        fields1 = np.zeros((n1,3))
-        fields1[not_mask1] = field
-        field = fields1
+    field = (deltaSin/norm_o4*eB.T / norm_12 * current / 10).T # m->mm, T->mT
 
-    if any_zero_segments: # broadcast into np.zeros
-        n0 = len(mask0)
-        fields0 = np.zeros((n0,3))
-        fields0[not_mask0] = field
-        field = fields0
+    # broadcast general case results into allocated vector
+    mask0[~mask0] = mask1
+    field_all[~mask0] = field
 
-    # return B
+    # return B or H
     if bh:
-        return field
+        return field_all
 
-    # return H (mT -> kA/m)
-    H = field*10/4/np.pi
-    return H
+    # H: mT -> kA/m
+    return field_all*10/4/np.pi
diff --git a/magpylib/_src/fields/field_BH_loop.py b/magpylib/_src/fields/field_BH_loop.py
index 5db1bb1d6..90e063ca9 100644
--- a/magpylib/_src/fields/field_BH_loop.py
+++ b/magpylib/_src/fields/field_BH_loop.py
@@ -2,89 +2,43 @@
 Implementations of analytical expressions for the magnetic field of
 a circular current loop. Computation details in function docstrings.
 """
-# pylint: disable=no-name-in-module
+# pylint: disable=no-name-in-module # reason scipy.special.ellipe import
 
 import numpy as np
 from scipy.special import ellipe
 from magpylib._src.fields.special_cel import cel_loop_stable
-
-
-def field_BH_loop(
-    bh: bool,
-    current: np.ndarray,
-    dia: np.ndarray,
-    pos_obs: np.ndarray
-    ) -> list:
-    """
-    Field of circular current loop.
-    - wraps fundamental implementation Smythe1950
-    - selects B or H
-    - sets singularity at wire to 0
-    - Cylinder CS <-> Cartesian CS
-    """
-    # pylint: disable=too-many-locals
-
-    x, y, z = pos_obs.T
-    r0 = dia/2
-    n = len(x)
-
-    # cylindrical coordinates ----------------------------------------------
-    r, phi = np.sqrt(x**2+y**2), np.arctan2(y, x)
-    pos_obs_cy = np.concatenate(((r,),(z,)),axis=0).T
-
-    # allocate fields as zeros ----------------------------------------------
-    #    singularities will become zero, other will be overwritten
-    Br_all = np.zeros(n)
-    Bz_all = np.zeros(n)
-
-    # determine singularity positions (pos_obs on wire) ------------------------
-    mask0 = (r==r0)*(z==0)
-
-    # forward only non-singularity observer positions for computation--------
-    Br, Bz = current_loop_B_Leitner2021(r0[~mask0], pos_obs_cy[~mask0]).T
-
-    # insert non-singular computations into total vectors----------------------
-    Br_all[~mask0] = Br
-    Bz_all[~mask0] = Bz
-
-    # transform field to cartesian CS -----------------------------------------
-    Bx = Br_all*np.cos(phi)
-    By = Br_all*np.sin(phi)
-    B_cart = np.concatenate(((Bx,),(By,),(Bz_all,)),axis=0) # ugly but fast
-
-    B_cart *= current/10
-
-    # B or H field
-    if bh:
-        return B_cart.T
-    return B_cart.T*10/4/np.pi
-
+from magpylib._src.utility import cart_to_cyl_coordinates, cyl_field_to_cart
+from magpylib._src.input_checks import check_field_input
 
 # ON INTERFACE
-def current_loop_B_Leitner2021(
-    radius: np.ndarray,
-    pos_obs: np.ndarray
+def current_loop_field(
+    current: np.ndarray,
+    diameter: np.ndarray,
+    observer: np.ndarray,
+    field='B'
     ) -> np.ndarray:
     """
-    B-field in cylindrical CS of circular line-current loop. The
-    current loop lies in the z=0 plane with the origin at its center.
-    The field is computed with 12-14 digits precision. Implementation
-    from [Leitner2021].
+    Computes the magnetic field of a circular line-current loop in Cartesian coordinates.
+    The loop lies in the z=0 plane with the coordinate origin at its center.
 
     Parameters
     ----------
-    radius: ndarray, shape (n,)
-        radius of current loop in units of [mm].
+    current: ndarray, shape (n,)
+        Electrical current in units of [A].
+
+    diameter: ndarray, shape (n,)
+        Diameter of loop in units of [mm].
 
-    pos_obs: ndarray, shape (n,2)
-        position of observer in cylindrical coordinates (r, z)
-        in units of [mm].
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
 
     Returns
     -------
-    B-field: ndarray
-        B-field of current loop in cylindrical coordinates (Br, Bz),
-        shape (n,2) in units of [mT].
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of current in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
 
     Examples
     --------
@@ -92,13 +46,14 @@ def current_loop_B_Leitner2021(
 
     >>> import numpy as np
     >>> import magpylib as magpy
-    >>> rad = np.array([1,2,3])
+    >>> cur = np.array([1,1,2])
+    >>> dia = np.array([2,4,6])
     >>> obs = np.array([(1,1), (2,2), (3,3)])
-    >>> B = magpy.lib.current_loop_B_Smythe1950(rad, obs)
+    >>> B = magpy.lib.current_loop_field(cur, dia, obs)
     >>> print(B)
-    [[1.14331448 0.96483239]
-     [0.57165724 0.48241619]
-     [0.38110483 0.3216108 ]]
+    [[0.11433145 0.09648324]
+     [0.05716572 0.04824162]
+     [0.07622097 0.06432216]]
 
     Notes
     -----
@@ -111,36 +66,71 @@ def current_loop_B_Leitner2021(
     loop," 2001.
 
     Ortner, "Feedback of Eddy Currents in Layered Materials for Magnetic Speed Sensing",
-    IEEE Transactions on Magnetics ( Volume: 53, Issue: 8, Aug. 2017)
+    IEEE Transactions on Magnetics ( Volume: 53, Issue: 8, Aug. 2017).
 
-    Leitner, "work in progress"
+    New numerically stable implementation based on [Ortner/Leitner wip]. Based
+    thereon the field is computed with >12 digits precision everywhere. On the
+    loop the result is set to (0,0,0).
     """
 
-    # inputs   -----------------------------------------------------------
-    r0 = radius
-    r, z = pos_obs.T
-    n = len(r0)
+    bh = check_field_input(field, 'current_loop_field()')
+
+    r, phi, z = cart_to_cyl_coordinates(observer)
+    rad = np.abs(diameter/2)
+    n = len(rad)
+
+    # define masks for special cases
+    mask_radius0 = rad==0
+    # rel pos deviation by 1e-15 to account for num errors (e.g. when rotating)
+    mask_on_loop = np.logical_and(abs(r-rad)<1e-15*rad, z==0)
+    mask_general = ~np.logical_or(mask_radius0, mask_on_loop)
+
+    # collect general case inputs
+    r = r[mask_general]
+    z = z[mask_general]
+    rad = rad[mask_general]
+    nX = len(rad)
 
-    # make dimensionless
-    rb = r/r0
-    zb = z/r0
+    # express through ratios (make dimensionless, avoid large/small input values)
+    rb = r/rad
+    zb = z/rad
 
     # pre-compute small quantities that might not be cached
     z2 = zb**2
     brack = (z2+(rb+1)**2)
     k2 = 4*rb/brack
-    pf = 1/np.sqrt(brack)/(1-k2)/r0
     xi = cel_loop_stable(k2)
 
-    # rb=0 requires special treatment
+    # rb=0 (on z-axis) requires special treatment because ellipe(x)/x and
+    # cel(x)/x for x->0 both appear in the expressions, both are finite at x=0
+    # but evaluation gives 0/0=nan
+    # To avoid treating x=0 as a special case (with masks), one would have to find
+    # designated algorithms for these expressions (which require additional evaluation
+    # and computation time because xi must be evaluated in any case).
+
     mask1 = rb==0
-    z_over_r = np.zeros(n)
-    k2_over_rb = np.ones(n)*4/(z2+1)
+    z_over_r = np.zeros(nX)
+    k2_over_rb = np.ones(nX)*4/(z2+1)
     z_over_r[~mask1] = zb[~mask1]/rb[~mask1]    # will be zero when r=0
-    k2_over_rb[~mask1] = k2[~mask1]/rb[~mask1] # will be zero when r=0
+    k2_over_rb[~mask1] = k2[~mask1]/rb[~mask1]  # will be zero when r=0
 
     # field components
+    pf = 1/np.sqrt(brack)/(1-k2)/rad
     Br = pf * z_over_r * xi
     Bz = pf * (k2_over_rb*ellipe(k2) - xi)
 
-    return np.array([Br, Bz]).T
+    # current and [mT] unit
+    Br_tot, Bz_tot = np.zeros((2,n))
+    Br_tot[mask_general] = Br
+    Bz_tot[mask_general] = Bz
+
+
+    # transform field to cartesian CS
+    Bx_tot, By_tot = cyl_field_to_cart(phi, Br_tot)
+    B_cart = (np.concatenate(((Bx_tot,),(By_tot,),(Bz_tot,)), axis=0)*current).T # ugly but fast
+
+    # B or H field
+    if bh:
+        return B_cart/10
+
+    return B_cart/4/np.pi
diff --git a/magpylib/_src/fields/field_BH_sphere.py b/magpylib/_src/fields/field_BH_sphere.py
index e0582a8e6..570147d2b 100644
--- a/magpylib/_src/fields/field_BH_sphere.py
+++ b/magpylib/_src/fields/field_BH_sphere.py
@@ -4,57 +4,85 @@
 """
 
 import numpy as np
+from magpylib._src.input_checks import check_field_input
 
 
-def field_BH_sphere(
-        bh: bool,
-        mag: np.ndarray,
-        dim: np.ndarray,
-        pos_obs: np.ndarray
-        ) -> np.ndarray:
+def magnet_sphere_field(
+    magnetization: np.ndarray,
+    diameter: np.ndarray,
+    observer: np.ndarray,
+    field='B'
+    )->np.ndarray:
     """
-    The B-field of a homogeneously magnetized spherical magnet corresponds to a dipole
-    field on the outside and is 2/3*mag in the inside (see e.g. "Theoretical Physics, Bertelmann")
+    Computes the magnetic field of a homogeneously magnetized spherical magnet
+    in Cartesian coordinates.
 
-    - separate mag=0 cases (returning 0)
-    - select B or H
-
-    Paramters:
+    Parameters
     ----------
-    - bh (boolean): True=B, False=H
-    - mag (ndarray Nx3): homogeneous magnetization vector in units of mT
-    - dim (ndarray Nx3): Sphere diameter
-    - pos_obs (ndarray Nx3): position of observer in units of mm
+    magnetization: ndarray, shape (n,3)
+        Homogeneous magnetization vector in units of [mT].
+
+    diameter: ndarray, shape (n,3)
+        Sphere diameter in units of [mm].
+
+    observer: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    field: str, default='B'
+        If 'B' return B-field in units of [mT], if 'H' return H-field in units of [kA/m].
+
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of magnet in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
 
-    Returns:
+    Examples
     --------
-    B/H-field (ndarray Nx3): magnetic field vectors at pos_obs in units of mT / kA/m
+    Compute the field of two different spherical magnets at position (1,1,1),
+    inside and outside.
+
+    >>> import numpy as np
+    >>> import magpylib as magpy
+    >>> dia = np.array([1,5])
+    >>> obs = np.array([(1,1,1), (1,1,1)])
+    >>> mag = np.array([(1,2,3), (0,0,3)])
+    >>> B = magpy.lib.magnet_sphere_Bfield(mag, dia, obs)
+    >>> print(B)
+    [[0.04009377 0.03207501 0.02405626]
+     [0.         0.         2.        ]]
+
+    Notes
+    -----
+    The field corresponds to a dipole field on the outside and is 2/3*mag
+    in the inside (see e.g. "Theoretical Physics, Bertelmann").
     """
 
-    x, y, z = np.copy(pos_obs.T)
+    bh = check_field_input(field, 'magnet_sphere_field()')
+
+    # all special cases r0=0 and mag=0 automatically covered
+
+    x, y, z = np.copy(observer.T)
     r = np.sqrt(x**2+y**2+z**2)   # faster than np.linalg.norm
-    radius = dim/2
+    r0 = abs(diameter)/2
 
     # inside field & allocate
-    B = mag*2/3
+    B = magnetization*2/3
 
-    # overwrite outside field
-    # special case mag=0 automatically reflected in field formulas
-    mask_out = r>=radius
+    # overwrite outside field entries
+    mask_out = (r>=r0)
 
-    mag1 = mag[mask_out]
-    poso1 = pos_obs[mask_out]
+    mag1 = magnetization[mask_out]
+    obs1 = observer[mask_out]
     r1 = r[mask_out]
-    dim1 = dim[mask_out]
-    val = (3*(np.sum(mag1*poso1,axis=1)*poso1.T)/r1**5 - mag1.T/r1**3)*dim1**3/24
-    B[mask_out] = val.T
+    r01 = r0[mask_out]
+
+    field_out = (3*(np.sum(mag1*obs1,axis=1)*obs1.T)/r1**5 - mag1.T/r1**3)*r01**3/3
+    B[mask_out] = field_out.T
 
-    # return B
     if bh:
         return B
 
     # adjust and return H
-    mask_in = ~mask_out
-    B[mask_in] = -mag[mask_in]/3
+    B[~mask_out] = -magnetization[~mask_out]/3
     H = B*10/4/np.pi
     return H
diff --git a/magpylib/_src/fields/field_wrap_BH_level1.py b/magpylib/_src/fields/field_wrap_BH_level1.py
index e7b487418..3b1afdc66 100644
--- a/magpylib/_src/fields/field_wrap_BH_level1.py
+++ b/magpylib/_src/fields/field_wrap_BH_level1.py
@@ -1,11 +1,11 @@
 import numpy as np
-from magpylib._src.fields.field_BH_cuboid import field_BH_cuboid
-from magpylib._src.fields.field_BH_cylinder import field_BH_cylinder
-from magpylib._src.fields.field_BH_cylinder_tile import field_BH_cylinder_tile
-from magpylib._src.fields.field_BH_sphere import field_BH_sphere
-from magpylib._src.fields.field_BH_dipole import field_BH_dipole
-from magpylib._src.fields.field_BH_loop import field_BH_loop
-from magpylib._src.fields.field_BH_line import field_BH_line, field_BH_line_from_vert
+from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
+from magpylib._src.fields.field_BH_cylinder import magnet_cylinder_field
+from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_field
+from magpylib._src.fields.field_BH_sphere import magnet_sphere_field
+from magpylib._src.fields.field_BH_dipole import dipole_field
+from magpylib._src.fields.field_BH_loop import current_loop_field
+from magpylib._src.fields.field_BH_line import current_line_field, field_BH_line_from_vert
 from magpylib._src.exceptions import MagpylibInternalError
 
 
@@ -29,7 +29,7 @@ def getBH_level1(**kwargs:dict) -> np.ndarray:
 
     # base inputs of all sources
     src_type = kwargs['source_type']
-    bh = kwargs['bh']      # True=B, False=H
+    field = kwargs['field']      # 'B' or 'H'
 
     rot = kwargs['orientation']    # only rotation object allowed as input
     pos = kwargs['position']
@@ -43,52 +43,49 @@ def getBH_level1(**kwargs:dict) -> np.ndarray:
     if src_type == 'Cuboid':
         mag = kwargs['magnetization']
         dim = kwargs['dimension']
-        B = field_BH_cuboid(bh, mag, dim, pos_rel_rot)
+        B = magnet_cuboid_field(mag, dim, pos_rel_rot, field=field)
 
     elif src_type == 'Cylinder':
         mag = kwargs['magnetization']
         dim = kwargs['dimension']
-        B = field_BH_cylinder(bh, mag, dim, pos_rel_rot)
+        B = magnet_cylinder_field(mag, dim, pos_rel_rot, field=field)
 
     elif src_type == 'CylinderSegment':
         mag = kwargs['magnetization']
         dim = kwargs['dimension']
-        r1,r2,h,phi1,phi2 = dim.T
-        z1, z2 = -h/2, h/2
-        dim = np.array([r1,r2,phi1,phi2,z1,z2]).T
-        B = field_BH_cylinder_tile(bh, mag, dim, pos_rel_rot)
+        B = magnet_cylinder_segment_field(mag, dim, pos_rel_rot, field=field)
 
     elif src_type == 'Sphere':
         mag = kwargs['magnetization']
         dia = kwargs['diameter']
-        B = field_BH_sphere(bh, mag, dia, pos_rel_rot)
+        B = magnet_sphere_field(mag, dia, pos_rel_rot, field=field)
 
     elif src_type == 'Dipole':
         moment = kwargs['moment']
-        B = field_BH_dipole(bh, moment, pos_rel_rot)
+        B = dipole_field(moment, pos_rel_rot, field=field)
 
     elif src_type == 'Loop':
         current = kwargs['current']
         dia = kwargs['diameter']
-        B = field_BH_loop(bh, current, dia, pos_rel_rot)
+        B = current_loop_field(current, dia, pos_rel_rot, field=field)
 
     elif src_type =='Line':
         current = kwargs['current']
         if 'vertices' in kwargs:
             vertices = kwargs['vertices']
-            B = field_BH_line_from_vert(bh, current, vertices, pos_rel_rot)
+            B = field_BH_line_from_vert(current, vertices, pos_rel_rot, field=field)
         else:
             pos_start = kwargs['segment_start']
             pos_end = kwargs['segment_end']
-            B = field_BH_line(bh, current, pos_start, pos_end, pos_rel_rot)
+            B = current_line_field(current, pos_start, pos_end, pos_rel_rot, field=field)
 
     elif src_type == 'CustomSource':
-        bh_key = 'B' if bh else 'H'
-        if kwargs[f'field_{bh_key}_lambda'] is not None:
-            B = kwargs[f'field_{bh_key}_lambda'](pos_rel_rot)
+        #bh_key = 'B' if bh else 'H'
+        if kwargs[f'field_{field}_lambda'] is not None:
+            B = kwargs[f'field_{field}_lambda'](pos_rel_rot)
         else:
             raise MagpylibInternalError(
-                f'{bh_key}-field calculation not implemented for CustomSource class'
+                f'{field}-field calculation not implemented for CustomSource class'
             )
 
     else:
diff --git a/magpylib/_src/fields/field_wrap_BH_level2.py b/magpylib/_src/fields/field_wrap_BH_level2.py
index b97a2814f..2136ee4db 100644
--- a/magpylib/_src/fields/field_wrap_BH_level2.py
+++ b/magpylib/_src/fields/field_wrap_BH_level2.py
@@ -32,8 +32,8 @@ def tile_dim_cylinder(group: list, n_pp: int):
     return dimv
 
 
-def tile_dim_cylinder_section(group: list, n_pp: int):
-    """ tile up cylinder section dimensions.
+def tile_dim_cylinder_segment(group: list, n_pp: int):
+    """ tile up cylinder segment dimensions.
     """
     dims = np.array([src.dimension for src in group])
     dimv = np.tile(dims, n_pp).reshape((-1, 5))
@@ -105,7 +105,7 @@ def get_src_dict(group: list, n_pix: int, n_pp: int, poso: np.ndarray) -> dict:
 
     elif src_type == 'CylinderSegment':
         magv = tile_mag(group, n_pp)
-        dimv = tile_dim_cylinder_section(group, n_pp)
+        dimv = tile_dim_cylinder_segment(group, n_pp)
         kwargs.update({'magnetization':magv,  'dimension':dimv})
 
     elif src_type == 'Dipole':
@@ -135,7 +135,7 @@ def get_src_dict(group: list, n_pix: int, n_pp: int, poso: np.ndarray) -> dict:
     return kwargs
 
 
-def getBH_level2(bh, sources, observers, sumup, squeeze, **kwargs) -> np.ndarray:
+def getBH_level2(sources, observers, **kwargs) -> np.ndarray:
     """...
 
     Parameters
@@ -143,13 +143,14 @@ def getBH_level2(bh, sources, observers, sumup, squeeze, **kwargs) -> np.ndarray
     - bh (bool): True=getB, False=getH
     - sources (src_obj or list): source object or 1D list of L sources/collections with similar
         pathlength M and/or 1.
-    - observers (sens_obj or list or pos_obs): pos_obs or sensor object or 1D list of K sensors with
-        similar pathlength M and/or 1 and sensor pixel of shape (N1,N2,...,3).
-    - sumup (bool): default=False returns [B1,B2,...] for every source, True returns sum(Bi)
+    - observers (sens_obj or list or pos_obs): pos_obs or sensor object or 1D list of K
+        sensors with similar pathlength M and/or 1 and sensor pixel of shape (N1,N2,...,3).
+    kwargs:
+    - 'sumup' (bool): False returns [B1,B2,...] for every source, True returns sum(Bi)
         for all sources.
-    - squeeze (bool): default=True, If True output is squeezed (axes of length 1 are eliminated)
-    - kwargs : if not empty, returns getBH_dict_level2 with source_type=sources and observer=
-        observers
+    - 'squeeze' (bool): True output is squeezed (axes of length 1 are eliminated)
+    - 'field' (str): 'B' computes B field, 'H' computes H-field
+    - getBH_dict inputs
 
     Returns
     -------
@@ -169,22 +170,24 @@ def getBH_level2(bh, sources, observers, sumup, squeeze, **kwargs) -> np.ndarray
     # pylint: disable=too-many-branches
     # pylint: disable=too-many-statements
 
-
     # CHECK AND FORMAT INPUT ---------------------------------------------------
     if isinstance(sources, str):
         return getBH_dict_level2(
-            bh = bh,
             source_type=sources,
             observer=observers,
-            sumup=sumup,
-            squeeze=squeeze,
             **kwargs
         )
-    if kwargs:
+
+    # bad user inputs mixing getBH_dict kwargs with object oriented interface
+    kwargs_check = kwargs.copy()
+    for popit in ['field', 'sumup', 'squeeze']:
+        kwargs_check.pop(popit)
+    if kwargs_check:
         raise MagpylibBadUserInput(
-            f"Keyword arguments {tuple(kwargs.keys())} are only allowed when the source is "
-            "defined by a string (e.g. sources='Cylinder')"
+            f"Keyword arguments {tuple(kwargs_check.keys())} are only allowed when the source "
+            "is defined by a string (e.g. sources='Cylinder')"
         )
+
     # format sources input:
     #   input: allow only bare src objects or 1D lists/tuple of src and col
     #   out: sources = ordered list of sources
@@ -265,7 +268,7 @@ def getBH_level2(bh, sources, observers, sumup, squeeze, **kwargs) -> np.ndarray
     groups = {}
     for ind,src in enumerate(src_list):
         if src._object_type=='CustomSource':
-            group_key = src.field_B_lambda if bh else src.field_H_lambda
+            group_key = src.field_B_lambda if kwargs['field']=='B' else src.field_H_lambda
         else:
             group_key = src._object_type
         if group_key not in groups:
@@ -279,7 +282,7 @@ def getBH_level2(bh, sources, observers, sumup, squeeze, **kwargs) -> np.ndarray
         lg = len(group['sources'])
         gr = group['sources']
         src_dict = get_src_dict(gr, n_pix, n_pp, poso)  # compute array dict for level1
-        B_group = getBH_level1(bh=bh, **src_dict)       # compute field
+        B_group = getBH_level1(field=kwargs['field'], **src_dict)       # compute field
         B_group = B_group.reshape((lg,m,n_pix,3))       # reshape (2% slower for large arrays)
         for i in range(lg):                             # put into dedicated positions in B
             B[group['order'][i]] = B_group[i]
@@ -311,19 +314,19 @@ def getBH_level2(bh, sources, observers, sumup, squeeze, **kwargs) -> np.ndarray
                 for j in range(m): # THIS LOOP IS EXTREMELY SLOW !!!! github issue #283
                     Bpart = B[:,j,i*k_pixel:(i+1)*k_pixel]           # select part
                     Bpart_flat = np.reshape(Bpart, (k_pixel*l0,3))   # flatten for rot package
-                    Bpart_flat_rot = sens._orientation[j].inv().apply(Bpart_flat)  # apply rotation
-                    B[:,j,i*k_pixel:(i+1)*k_pixel] = np.reshape(Bpart_flat_rot, (l0,k_pixel,3)) # ov
+                    Bpart_flat_rot = sens._orientation[j].inv().apply(Bpart_flat)  # apply rotat
+                    B[:,j,i*k_pixel:(i+1)*k_pixel] = np.reshape(Bpart_flat_rot, (l0,k_pixel,3))
 
     # rearrange sensor-pixel shape
     sens_px_shape = (k,) + pix_shape
     B = B.reshape((l0,m)+sens_px_shape)
 
     #
-    if sumup:
+    if kwargs['sumup']:
         B = np.sum(B, axis=0, keepdims=True)
 
     # reduce all size-1 levels
-    if squeeze:
+    if kwargs['squeeze']:
         B = np.squeeze(B)
 
     # reset tiled objects
diff --git a/magpylib/_src/fields/field_wrap_BH_level2_dict.py b/magpylib/_src/fields/field_wrap_BH_level2_dict.py
index 484933e9d..c04df62f5 100644
--- a/magpylib/_src/fields/field_wrap_BH_level2_dict.py
+++ b/magpylib/_src/fields/field_wrap_BH_level2_dict.py
@@ -138,205 +138,207 @@ def getBH_dict_level2(**kwargs: dict) -> np.ndarray:
     return B
 
 
-def getB_dict(**kwargs):
-    """
-    B-Field computation in units of [mT] from a dictionary of input vectors of
-    length N.
-
-    This function avoids the object-oriented Magpylib interface and gives direct
-    access to the field implementations. It is the fastest way to compute fields
-    with Magpylib.
-
-    "Static" inputs of shape (x,) will automatically be tiled to shape (N,x) to
-    fit with other inputs.
-
-    Required inputs depend on chosen source_type!
-
-    Parameters
-    ----------
-    source_type: string
-        Source type for computation. Must be either 'Cuboid', 'Cylinder', 'Cylinder_old', 'Sphere',
-        'Dipole', 'Loop' or 'Line'. Expected input parameters depend on source_type.
-
-    position: array_like, shape (3,) or (N,3), default=(0,0,0)
-        Source positions in units of [mm].
-
-    orientation: scipy Rotation object, default=unit rotation
-        Source rotations relative to the initial state (see object docstrings).
-
-    observer: array_like, shape (3,) or (N,3)
-        Observer positions in units of [mm].
+# def getB_dict(**kwargs):
+#     """
+#     B-Field computation in units of [mT] from a dictionary of input vectors of
+#     length N.
+
+#     This function avoids the object-oriented Magpylib interface and gives direct
+#     access to the field implementations. It is the fastest way to compute fields
+#     with Magpylib.
+
+#     "Static" inputs of shape (x,) will automatically be tiled to shape (N,x) to
+#     fit with other inputs.
+
+#     Required inputs depend on chosen source_type!
+
+#     Parameters
+#     ----------
+#     source_type: string
+#         Source type for computation. Must be either 'Cuboid', 'Cylinder', 'Cylinder_old',
+#           'Sphere',
+#         'Dipole', 'Loop' or 'Line'. Expected input parameters depend on source_type.
+
+#     position: array_like, shape (3,) or (N,3), default=(0,0,0)
+#         Source positions in units of [mm].
+
+#     orientation: scipy Rotation object, default=unit rotation
+#         Source rotations relative to the initial state (see object docstrings).
+
+#     observer: array_like, shape (3,) or (N,3)
+#         Observer positions in units of [mm].
+
+#     squeeze: bool, default=True
+#         If True, the output is squeezed, i.e. all axes of length 1 in the output are eliminated.
+
+#     magnetization: array_like, shape (3,) or (N,3)
+#         Only `source_type in ('Cuboid', 'Cylinder', 'Sphere')`! Magnetization vector (mu0*M) or
+#         remanence field of homogeneous magnet magnetization in units of [mT].
 
-    squeeze: bool, default=True
-        If True, the output is squeezed, i.e. all axes of length 1 in the output are eliminated.
+#     moment:  array_like, shape (3,) or (N,3)
+#         Only `source_type = 'Moment'`! Magnetic dipole moment in units of [mT*mm^3]. For
+#         homogeneous magnets the relation is moment = magnetization*volume.
 
-    magnetization: array_like, shape (3,) or (N,3)
-        Only `source_type in ('Cuboid', 'Cylinder', 'Sphere')`! Magnetization vector (mu0*M) or
-        remanence field of homogeneous magnet magnetization in units of [mT].
+#     current: array_like, shape (N,)
+#         Only `source_type in ('Line', 'Loop')`! Current flowing in loop in units of [A].
+
+#     dimension: array_like
+#         Only `source_type in ('Cuboid', 'Cylinder', 'CylinderSegment')`! Magnet dimension
+#         input in units of [mm]. Dimension format of sources similar as in object oriented
+#         interface.
 
-    moment:  array_like, shape (3,) or (N,3)
-        Only `source_type = 'Moment'`! Magnetic dipole moment in units of [mT*mm^3]. For
-        homogeneous magnets the relation is moment = magnetization*volume.
+#     diameter: array_like, shape (N)
+#         Only `source_type in (Sphere, Loop)`! Diameter of source in units of [mm].
+
+#     segment_start: array_like, shape (N,3)
+#         Only `source_type = 'Line'`! Start positions of line current segments in units of [mm].
 
-    current: array_like, shape (N,)
-        Only `source_type in ('Line', 'Loop')`! Current flowing in loop in units of [A].
+#     segment_end: array_like, shape (N,3)
+#         Only `source_type = 'Line'`! End positions of line current segments in units of [mm].
 
-    dimension: array_like
-        Only `source_type in ('Cuboid', 'Cylinder', 'CylinderSegment')`! Magnet dimension
-        input in units of [mm]. Dimension format of sources similar as in object oriented
-        interface.
+#     Returns
+#     -------
+#     B-field: ndarray, shape (N,3)
+#         B-field generated by sources at observer positions in units of [mT].
 
-    diameter: array_like, shape (N)
-        Only `source_type in (Sphere, Loop)`! Diameter of source in units of [mm].
+#     Examples
+#     --------
 
-    segment_start: array_like, shape (N,3)
-        Only `source_type = 'Line'`! Start positions of line current segments in units of [mm].
-
-    segment_end: array_like, shape (N,3)
-        Only `source_type = 'Line'`! End positions of line current segments in units of [mm].
-
-    Returns
-    -------
-    B-field: ndarray, shape (N,3)
-        B-field generated by sources at observer positions in units of [mT].
-
-    Examples
-    --------
-
-    Three-fold evaluation of the dipole field. For each computation the moment is (100,100,100).
-
-    >>> import magpylib as magpy
-    >>> B = magpy.getB_dict(
-    >>>     source_type='Dipole',
-    >>>     position=[(1,2,3), (2,3,4), (3,4,5)],
-    >>>     moment=(100,100,100),
-    >>>     observer=[(1,1,1), (2,2,2), (3,3,3)])
-    >>> print(B)
-    [[-0.71176254  0.56941003  1.85058261]
-     [-0.71176254  0.56941003  1.85058261]
-     [-0.71176254  0.56941003  1.85058261]]
-
-    Six-fold evaluation of a Cuboid magnet field with increasing size and magnetization
-    of the magnet. Position and orientation are by default (0,0,0) and unit-orientation,
-    respectively. The observer position is (1,2,3) for each evaluation.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.getB_dict(
-    >>>     source_type='Cuboid',
-    >>>     magnetization = [(0,0,m) for m in np.linspace(500,1000,6)],
-    >>>     dimension = [(a,a,a) for a in np.linspace(1,2,6)],
-    >>>     observer=(1,2,3))
-    >>> print(B)
-    [[ 0.48818967  0.97689261  0.70605984]
-     [ 1.01203491  2.02636222  1.46575704]
-     [ 1.87397714  3.756164    2.72063422]
-     [ 3.19414311  6.41330652  4.65485356]
-     [ 5.10909461 10.2855981   7.4881383 ]
-     [ 7.76954697 15.70382556 11.48192812]]
-
-
-    """
-    return getBH_dict_level2(bh=True, **kwargs)
-
-
-# ON INTERFACE
-def getH_dict(**kwargs):
-    """
-    H-Field computation in units of [kA/m] from a dictionary of input vectors of
-    length N.
+#     Three-fold evaluation of the dipole field. For each computation the moment is (100,100,100).
 
-    This function avoids the object-oriented Magpylib interface and gives direct
-    access to the field implementations. It is the fastest way to compute fields
-    with Magpylib.
+#     >>> import magpylib as magpy
+#     >>> B = magpy.getB_dict(
+#     >>>     source_type='Dipole',
+#     >>>     position=[(1,2,3), (2,3,4), (3,4,5)],
+#     >>>     moment=(100,100,100),
+#     >>>     observer=[(1,1,1), (2,2,2), (3,3,3)])
+#     >>> print(B)
+#     [[-0.71176254  0.56941003  1.85058261]
+#      [-0.71176254  0.56941003  1.85058261]
+#      [-0.71176254  0.56941003  1.85058261]]
 
-    "Static" inputs of shape (x,) will automatically be tiled to shape (N,x) to
-    fit with other inputs.
+#     Six-fold evaluation of a Cuboid magnet field with increasing size and magnetization
+#     of the magnet. Position and orientation are by default (0,0,0) and unit-orientation,
+#     respectively. The observer position is (1,2,3) for each evaluation.
 
-    Required inputs depend on chosen source_type!
-
-    Parameters
-    ----------
-    source_type: string
-        Source type for computation. Must be either 'Cuboid', 'Cylinder','Cylinder_old', 'Sphere',
-        'Dipole', 'Loop' or 'Line'. Expected input parameters depend on source_type.
-
-    position: array_like, shape (3,) or (N,3), default=(0,0,0)
-        Source positions in units of [mm].
-
-    orientation: scipy Rotation object, default=unit rotation
-        Source rotations relative to the initial state (see object docstrings).
-
-    observer: array_like, shape (3,) or (N,3)
-        Observer positions in units of [mm].
-
-    squeeze: bool, default=True
-        If True, the output is squeezed, i.e. all axes of length 1 in the output are eliminated.
-
-    magnetization: array_like, shape (3,) or (N,3)
-        Only `source_type in ('Cuboid', 'Cylinder', 'Sphere')`! Magnetization vector (mu0*M) or
-        remanence field of homogeneous magnet magnetization in units of [mT].
-
-    moment:  array_like, shape (3,) or (N,3)
-        Only `source_type = 'Moment'`! Magnetic dipole moment in units of [mT*mm^3]. For
-        homogeneous magnets the relation is moment = magnetization*volume.
-
-    current: array_like, shape (N,)
-        Only `source_type in ('Line', 'Loop')`! Current flowing in loop in units of [A].
-
-    dimension: array_like
-        Only `source_type in ('Cuboid', 'Cylinder', 'CylinderSegment')`! Magnet dimension
-        input in units of [mm]. Dimension format of sources similar as in object oriented
-        interface.
-
-    diameter: array_like, shape (N)
-        Only `source_type in (Sphere, Loop)`! Diameter of source in units of [mm].
-
-    segment_start: array_like, shape (N,3)
-        Only `source_type = 'Line'`! Start positions of line current segments in units of [mm].
-
-    segment_end: array_like, shape (N,3)
-        Only `source_type = 'Line'`! End positions of line current segments in units of [mm].
-
-    Returns
-    -------
-    H-field: ndarray, shape (N,3)
-        H-field generated by sources at observer positions in units of [kA/m].
-
-    Examples
-    --------
-
-    Three-fold evaluation of the dipole field. For each computation the moment is (100,100,100).
-
-    >>> import magpylib as magpy
-    >>> H = magpy.getH_dict(
-    >>>     source_type='Dipole',
-    >>>     position=[(1,2,3), (2,3,4), (3,4,5)],
-    >>>     moment=(100,100,100),
-    >>>     observer=[(1,1,1), (2,2,2), (3,3,3)])
-    >>> print(H)
-    [[-0.56640264  0.45312211  1.47264685]
-     [-0.56640264  0.45312211  1.47264685]
-     [-0.56640264  0.45312211  1.47264685]]
-
-    Six-fold evaluation of a Cuboid magnet field with increasing size and magnetization
-    of the magnet. Position and orientation are (0,0,0) and unit-orientation, respectively,
-    by default. The observer position is (1,2,3) for each evaluation.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> H = magpy.getH_dict(
-    >>>     source_type='Cuboid',
-    >>>     magnetization = [(0,0,m) for m in np.linspace(500,1000,6)],
-    >>>     dimension = [(a,a,a) for a in np.linspace(1,2,6)],
-    >>>     observer=(1,2,3))
-    >>> print(H)
-    [[ 0.388489    0.77738644  0.56186457]
-     [ 0.80535179  1.61252782  1.16641239]
-     [ 1.49126363  2.98906034  2.16501192]
-     [ 2.54181833  5.10354717  3.70421476]
-     [ 4.06568831  8.1850189   5.95887113]
-     [ 6.18280903 12.49670731  9.13702808]]
-
-    """
-    return getBH_dict_level2(bh=False, **kwargs)
+#     >>> import numpy as np
+#     >>> import magpylib as magpy
+#     >>> B = magpy.getB_dict(
+#     >>>     source_type='Cuboid',
+#     >>>     magnetization = [(0,0,m) for m in np.linspace(500,1000,6)],
+#     >>>     dimension = [(a,a,a) for a in np.linspace(1,2,6)],
+#     >>>     observer=(1,2,3))
+#     >>> print(B)
+#     [[ 0.48818967  0.97689261  0.70605984]
+#      [ 1.01203491  2.02636222  1.46575704]
+#      [ 1.87397714  3.756164    2.72063422]
+#      [ 3.19414311  6.41330652  4.65485356]
+#      [ 5.10909461 10.2855981   7.4881383 ]
+#      [ 7.76954697 15.70382556 11.48192812]]
+
+
+#     """
+#     return getBH_dict_level2(field='B', **kwargs)
+
+
+# # ON INTERFACE
+# def getH_dict(**kwargs):
+#     """
+#     H-Field computation in units of [kA/m] from a dictionary of input vectors of
+#     length N.
+
+#     This function avoids the object-oriented Magpylib interface and gives direct
+#     access to the field implementations. It is the fastest way to compute fields
+#     with Magpylib.
+
+#     "Static" inputs of shape (x,) will automatically be tiled to shape (N,x) to
+#     fit with other inputs.
+
+#     Required inputs depend on chosen source_type!
+
+#     Parameters
+#     ----------
+#     source_type: string
+#         Source type for computation. Must be either 'Cuboid', 'Cylinder','Cylinder_old',
+#           'Sphere',
+#         'Dipole', 'Loop' or 'Line'. Expected input parameters depend on source_type.
+
+#     position: array_like, shape (3,) or (N,3), default=(0,0,0)
+#         Source positions in units of [mm].
+
+#     orientation: scipy Rotation object, default=unit rotation
+#         Source rotations relative to the initial state (see object docstrings).
+
+#     observer: array_like, shape (3,) or (N,3)
+#         Observer positions in units of [mm].
+
+#     squeeze: bool, default=True
+#         If True, the output is squeezed, i.e. all axes of length 1 in the output are eliminated.
+
+#     magnetization: array_like, shape (3,) or (N,3)
+#         Only `source_type in ('Cuboid', 'Cylinder', 'Sphere')`! Magnetization vector (mu0*M) or
+#         remanence field of homogeneous magnet magnetization in units of [mT].
+
+#     moment:  array_like, shape (3,) or (N,3)
+#         Only `source_type = 'Moment'`! Magnetic dipole moment in units of [mT*mm^3]. For
+#         homogeneous magnets the relation is moment = magnetization*volume.
+
+#     current: array_like, shape (N,)
+#         Only `source_type in ('Line', 'Loop')`! Current flowing in loop in units of [A].
+
+#     dimension: array_like
+#         Only `source_type in ('Cuboid', 'Cylinder', 'CylinderSegment')`! Magnet dimension
+#         input in units of [mm]. Dimension format of sources similar as in object oriented
+#         interface.
+
+#     diameter: array_like, shape (N)
+#         Only `source_type in (Sphere, Loop)`! Diameter of source in units of [mm].
+
+#     segment_start: array_like, shape (N,3)
+#         Only `source_type = 'Line'`! Start positions of line current segments in units of [mm].
+
+#     segment_end: array_like, shape (N,3)
+#         Only `source_type = 'Line'`! End positions of line current segments in units of [mm].
+
+#     Returns
+#     -------
+#     H-field: ndarray, shape (N,3)
+#         H-field generated by sources at observer positions in units of [kA/m].
+
+#     Examples
+#     --------
+
+#     Three-fold evaluation of the dipole field. For each computation the moment is (100,100,100).
+
+#     >>> import magpylib as magpy
+#     >>> H = magpy.getH_dict(
+#     >>>     source_type='Dipole',
+#     >>>     position=[(1,2,3), (2,3,4), (3,4,5)],
+#     >>>     moment=(100,100,100),
+#     >>>     observer=[(1,1,1), (2,2,2), (3,3,3)])
+#     >>> print(H)
+#     [[-0.56640264  0.45312211  1.47264685]
+#      [-0.56640264  0.45312211  1.47264685]
+#      [-0.56640264  0.45312211  1.47264685]]
+
+#     Six-fold evaluation of a Cuboid magnet field with increasing size and magnetization
+#     of the magnet. Position and orientation are (0,0,0) and unit-orientation, respectively,
+#     by default. The observer position is (1,2,3) for each evaluation.
+
+#     >>> import numpy as np
+#     >>> import magpylib as magpy
+#     >>> H = magpy.getH_dict(
+#     >>>     source_type='Cuboid',
+#     >>>     magnetization = [(0,0,m) for m in np.linspace(500,1000,6)],
+#     >>>     dimension = [(a,a,a) for a in np.linspace(1,2,6)],
+#     >>>     observer=(1,2,3))
+#     >>> print(H)
+#     [[ 0.388489    0.77738644  0.56186457]
+#      [ 0.80535179  1.61252782  1.16641239]
+#      [ 1.49126363  2.98906034  2.16501192]
+#      [ 2.54181833  5.10354717  3.70421476]
+#      [ 4.06568831  8.1850189   5.95887113]
+#      [ 6.18280903 12.49670731  9.13702808]]
+
+#     """
+#     return getBH_dict_level2(field='H', **kwargs)
diff --git a/magpylib/_src/fields/field_wrap_BH_level3.py b/magpylib/_src/fields/field_wrap_BH_level3.py
index 2e2f16b01..7b9e620ba 100644
--- a/magpylib/_src/fields/field_wrap_BH_level3.py
+++ b/magpylib/_src/fields/field_wrap_BH_level3.py
@@ -133,7 +133,7 @@ def getB(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
     [[3.99833439 3.99833439 3.46340502]
      [0.09728187 0.34647784 0.52129178]]
     """
-    return getBH_level2(True, sources, observers, sumup, squeeze, **kwargs)
+    return getBH_level2(sources, observers, sumup=sumup, squeeze=squeeze, field='B', **kwargs)
 
 
 # ON INTERFACE
@@ -268,4 +268,4 @@ def getH(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
     [[3.18177341 3.18177341 2.75609015]
      [0.07741445 0.27571831 0.41483082]]
     """
-    return getBH_level2(False, sources, observers, sumup, squeeze, **kwargs)
+    return getBH_level2(sources, observers, sumup=sumup, squeeze=squeeze, field='H', **kwargs)
diff --git a/magpylib/_src/fields/special_cel.py b/magpylib/_src/fields/special_cel.py
index 9618a001c..00011d8bf 100644
--- a/magpylib/_src/fields/special_cel.py
+++ b/magpylib/_src/fields/special_cel.py
@@ -134,36 +134,70 @@ def ellipticPi(x, y):
 def cel_loop_stable(k2):
     """
     numerically stabilized version of the function
-        xi_loop = (2-k**2)E(k**2)-2(1-k**2)K(k**2)
-        needed for the current.loop fields. See paper
-        Leitner2021
+    xi_loop = (2-k**2)E(k**2)-2(1-k**2)K(k**2) = cel(np.sqrt(1-k2), 1, k2, k2*(k2-1))
+    needed for the circular current loop field. This modicfication of the cel()
+    algorithm is numerically stable when k2 becomes small.
     """
-    n = len(k2)
-    result = np.empty(n)
+    # if k2 == 0: return 0    # on axis
 
-    mask1 = (k2 > 0.04)
-
-    if np.any(mask1):
-        k2m = k2[mask1]
-        result[mask1] = cel(np.sqrt(1-k2m), np.ones(np.sum(mask1)), k2m, k2m*(k2m-1))
-    if np.any(~mask1):
-        result[~mask1] = cel_loop_taylor0(k2[~mask1])
-    return result
+    pp = 1-k2          # allocate pp and use temporarily for 1-k2
+    k  = np.sqrt(pp)
+    kk = 1+k           # allocate pp and use temporarily for 1+k
 
+    g  = 1 if isinstance(k2,(float,int)) else np.ones(len(k2))
+    cc = k2**2
+    ss = 2*cc*pp/(kk-k2)
+    pp = kk
+    em = kk
+    kk = k
 
-def cel_loop_taylor0(k2):
-    """
-    taylor expansion of the function xi_loop about k2=0
-        See paper Leitner2021
-    """
-    C2 = 0.5890486225480862
-    C3 = 0.1472621556370216
-    C4 = 0.06902913545485386
-    C5 = 0.04026699568199808
-    C6 = 0.02642521591631124
-    C7 = 0.01868640268367724
-    C8 = 0.01391747699878044
-    C9 = 0.01076947624905629
-    C10 = 0.008581926385966734
-    return (C2*k2**2 + C3*k2**3 + C4*k2**4 + C5*k2**5
-        + C6*k2**6 + C7*k2**7 + C8*k2**8 + C9*k2**9 + C10*k2**10)
+    errtol = .000001
+    while np.any(abs(g-k) > g*errtol):
+        k = 2*np.sqrt(kk)
+        kk = k*em
+        f = cc
+        cc = cc + ss/pp
+        g = kk/pp
+        ss = 2*(ss + f*g)
+        pp = g + pp
+        g = em
+        em = k+em
+    return (np.pi/2)*(ss+cc*em)/(em*(em+pp))
+
+
+# def cel_loop_stable_old(k2):
+#     """
+#     numerically stabilized version of the function
+#         xi_loop = (2-k**2)E(k**2)-2(1-k**2)K(k**2)
+#         needed for the current.loop fields. See paper
+#         Leitner2021
+#     """
+#     n = len(k2)
+#     result = np.empty(n)
+
+#     mask1 = (k2 > 0.04)
+
+#     if np.any(mask1):
+#         k2m = k2[mask1]
+#         result[mask1] = cel(np.sqrt(1-k2m), np.ones(np.sum(mask1)), k2m, k2m*(k2m-1))
+#     if np.any(~mask1):
+#         result[~mask1] = cel_loop_taylor0(k2[~mask1])
+#     return result
+
+
+# def cel_loop_taylor0(k2):
+#     """
+#     taylor expansion of the function xi_loop about k2=0
+#         See paper Leitner2021
+#     """
+#     C2 = 0.5890486225480862
+#     C3 = 0.1472621556370216
+#     C4 = 0.06902913545485386
+#     C5 = 0.04026699568199808
+#     C6 = 0.02642521591631124
+#     C7 = 0.01868640268367724
+#     C8 = 0.01391747699878044
+#     C9 = 0.01076947624905629
+#     C10 = 0.008581926385966734
+#     return (C2*k2**2 + C3*k2**3 + C4*k2**4 + C5*k2**5
+#         + C6*k2**6 + C7*k2**7 + C8*k2**8 + C9*k2**9 + C10*k2**10)
diff --git a/magpylib/_src/input_checks.py b/magpylib/_src/input_checks.py
index c30d44a05..2fb7d4b1d 100644
--- a/magpylib/_src/input_checks.py
+++ b/magpylib/_src/input_checks.py
@@ -9,6 +9,19 @@
 )
 
 
+def check_field_input(inp, origin):
+    """
+    check field input
+    """
+    if isinstance(inp, str):
+        if inp == 'B':
+            return True
+        if inp == 'H':
+            return False
+    msg = f'{origin} input can only be "field=B" or "field=H".'
+    raise MagpylibBadUserInput(msg)
+
+
 def check_dimensions(sources):
     """check if all sources have dimension (or similar) initialized"""
     # pylint: disable=protected-access
@@ -76,14 +89,14 @@ def check_rot_type(inp):
 
 def check_start_type(start):
     """start input must be int or str"""
-    if not (isinstance(start, int) or start == 'auto'):
+    if not (isinstance(start, (int, np.int_)) or start == 'auto'):
         msg = 'start input must be int or str ("auto")'
         raise MagpylibBadUserInput(msg)
 
 
 def check_angle_type(angle):
     """angle input must be scalar or vector"""
-    if not isinstance(angle, (int, float, list, tuple, np.ndarray)):
+    if not isinstance(angle, (int, float, list, tuple, np.ndarray, np.int_, np.float_)):
         msg = (
             "angle input must be scalar (int, float) or vector (list, tuple, ndarray)."
         )
@@ -128,7 +141,7 @@ def check_degree_type(deg):
 
 def check_scalar_type(inp, origin):
     """scalar input must be int or float or nan"""
-    if not (isinstance(inp, (int, float)) or inp is None):
+    if not (isinstance(inp, (int, float, np.int_, np.float_)) or inp is None):
         msg = origin + " input must be scalar (int or float)."
         raise MagpylibBadUserInput(msg)
 
@@ -151,16 +164,6 @@ def check_vector_type(inp, origin):
         raise MagpylibBadUserInput(msg)
 
 
-# def check_vector_init(inp, origin):
-#     """
-#     - check if vector input was initialized (former None vector)
-#     - return error msg with reference to origin
-#     """
-#     if None in inp:
-#         msg = origin + ' must be initialized.'
-#         raise MagpylibBadUserInput(msg)
-
-
 def check_vector_format(inp, shape, origin):
     """
     - check if vector input has correct format
diff --git a/magpylib/_src/obj_classes/class_BaseGetBH.py b/magpylib/_src/obj_classes/class_BaseGetBH.py
index 71b4baba5..edbf42bba 100644
--- a/magpylib/_src/obj_classes/class_BaseGetBH.py
+++ b/magpylib/_src/obj_classes/class_BaseGetBH.py
@@ -48,7 +48,7 @@ def getB(self, *observers, squeeze=True):
         --------
         """
         observers = format_star_input(observers)
-        return getBH_level2(True, self, observers, False, squeeze)
+        return getBH_level2(self, observers, sumup=False, squeeze=squeeze, field='B')
 
     # ------------------------------------------------------------------
     # INTERFACE
@@ -81,4 +81,4 @@ def getH(self, *observers, squeeze=True):
         --------
         """
         observers = format_star_input(observers)
-        return getBH_level2(False, self, observers, False, squeeze)
+        return getBH_level2(self, observers, sumup=False, squeeze=squeeze, field='H')
diff --git a/magpylib/_src/obj_classes/class_Collection.py b/magpylib/_src/obj_classes/class_Collection.py
index 34e0fdfb5..b43929d09 100644
--- a/magpylib/_src/obj_classes/class_Collection.py
+++ b/magpylib/_src/obj_classes/class_Collection.py
@@ -382,7 +382,7 @@ def getB(self, *children, squeeze=True):
 
         sources, sensors = self._validate_getBH_inputs(*children)
 
-        return getBH_level2(True, sources, sensors, sumup=False, squeeze=squeeze)
+        return getBH_level2(sources, sensors, sumup=False, squeeze=squeeze, field='B')
 
     def getH(self, *children, squeeze=True):
         """
@@ -411,7 +411,7 @@ def getH(self, *children, squeeze=True):
 
         sources, sensors = self._validate_getBH_inputs(*children)
 
-        return getBH_level2(False, sources, sensors, sumup=False, squeeze=squeeze)
+        return getBH_level2(sources, sensors, sumup=False, squeeze=squeeze, field='H')
 
 
 class Collection(BaseGeo, BaseCollection):
diff --git a/magpylib/_src/obj_classes/class_Sensor.py b/magpylib/_src/obj_classes/class_Sensor.py
index 85511b0e5..30ad50415 100644
--- a/magpylib/_src/obj_classes/class_Sensor.py
+++ b/magpylib/_src/obj_classes/class_Sensor.py
@@ -213,7 +213,7 @@ def getB(self, *sources, sumup=False, squeeze=True):
 
         """
         sources = format_star_input(sources)
-        return getBH_level2(True, sources, self, sumup, squeeze)
+        return getBH_level2(sources, self, sumup=sumup, squeeze=squeeze, field='B')
 
     def getH(self, *sources, sumup=False, squeeze=True):
         """
@@ -287,4 +287,4 @@ def getH(self, *sources, sumup=False, squeeze=True):
 
         """
         sources = format_star_input(sources)
-        return getBH_level2(False, sources, self, sumup, squeeze)
+        return getBH_level2(sources, self, sumup=sumup, squeeze=squeeze, field='H')
diff --git a/magpylib/_src/obj_classes/class_mag_CylinderSegment.py b/magpylib/_src/obj_classes/class_mag_CylinderSegment.py
index 2f57af65d..764c142a2 100644
--- a/magpylib/_src/obj_classes/class_mag_CylinderSegment.py
+++ b/magpylib/_src/obj_classes/class_mag_CylinderSegment.py
@@ -115,7 +115,7 @@ def dimension(self):
 
     @dimension.setter
     def dimension(self, dim):
-        """Set Cylinder dimension (d1,d2,h,phi1,phi2), shape (5,), [mm, deg]."""
+        """Set Cylinder dimension (r1,r2,h,phi1,phi2), shape (5,), [mm, deg]."""
         # input type check
         if Config.checkinputs:
             check_vector_type(dim, "dimension")
diff --git a/magpylib/_src/obj_classes/class_mag_Sphere.py b/magpylib/_src/obj_classes/class_mag_Sphere.py
index eb732c1cd..d7459f845 100644
--- a/magpylib/_src/obj_classes/class_mag_Sphere.py
+++ b/magpylib/_src/obj_classes/class_mag_Sphere.py
@@ -109,7 +109,7 @@ def diameter(self, dia):
         """Set Sphere diameter, float, [mm]."""
         # input type check
         if Config.checkinputs:
-            check_scalar_type(dia, "diameter")
+            check_scalar_type(dia, "Sphere diameter")
 
         # secure type
         if dia is None:
diff --git a/magpylib/_src/utility.py b/magpylib/_src/utility.py
index 4c8c60dc6..822bccdb0 100644
--- a/magpylib/_src/utility.py
+++ b/magpylib/_src/utility.py
@@ -67,16 +67,6 @@ def wrong_obj_msg(*objs, allow="sources"):
     return msg
 
 
-def close(arg1: np.ndarray, arg2: np.ndarray) -> np.ndarray:
-    """
-    determine if arg1 and arg2 lie close to each other
-    input: ndarray, shape (n,) or numpy-interpretable scalar
-    output: ndarray, dtype=bool
-    """
-    EDGESIZE = Config.edgesize
-    return np.isclose(arg1, arg2, rtol=0, atol=EDGESIZE)
-
-
 def format_star_input(inp):
     """
     *inputs are always wrapped in tuple. Formats *inputs of form "src", "src, src"
@@ -342,3 +332,27 @@ def unit_prefix(number, unit="", precision=3, char_between="") -> str:
         digits = 0
     new_number_str = f"{number / 10 ** digits:.{precision}g}"
     return f"{new_number_str}{char_between}{prefix}{unit}"
+
+
+def cart_to_cyl_coordinates(observer):
+    """
+    cartesian observer positions to cylindrical coordinates
+    observer: ndarray, shape (n,3)
+    """
+    x, y, z = observer.T
+    r, phi = np.sqrt(x**2+y**2), np.arctan2(y, x)
+    return r, phi, z
+
+
+def cyl_field_to_cart(phi, Br, Bphi=None):
+    """
+    transform Br,Bphi to Bx, By
+    """
+    if Bphi is not None:
+        Bx = Br*np.cos(phi) - Bphi*np.sin(phi)
+        By = Br*np.sin(phi) + Bphi*np.cos(phi)
+    else:
+        Bx = Br*np.cos(phi)
+        By = Br*np.sin(phi)
+
+    return Bx,By
diff --git a/magpylib/core/__init__.py b/magpylib/core/__init__.py
new file mode 100644
index 000000000..ca1ea3abf
--- /dev/null
+++ b/magpylib/core/__init__.py
@@ -0,0 +1,21 @@
+"""
+This sub-package gives direct access to our field implementations.
+"""
+
+__all__ = [
+    'dipole_field',
+    'current_loop_field',
+    'current_line_field',
+    'magnet_sphere_field',
+    'magnet_cuboid_field',
+    'magnet_cylinder_field',
+    'magnet_cylinder_segment_field',
+]
+
+from magpylib._src.fields.field_BH_dipole import dipole_field
+from magpylib._src.fields.field_BH_loop import current_loop_field
+from magpylib._src.fields.field_BH_line import current_line_field
+from magpylib._src.fields.field_BH_sphere import magnet_sphere_field
+from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
+from magpylib._src.fields.field_BH_cylinder import magnet_cylinder_field
+from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_field


Test patch:
diff --git a/tests/test_defaults.py b/tests/test_defaults.py
index 449793d29..a6c69c7df 100644
--- a/tests/test_defaults.py
+++ b/tests/test_defaults.py
@@ -9,8 +9,6 @@
 
 bad_inputs = {
     "checkinputs": (-1,),  # bool
-    "edgesize": (0,),  # float>0
-    "itercylinder": (0.1,),  # int>0
     "display_autosizefactor": (0,),  # float>0
     "display_animation_maxfps": (0,),  # int>0
     "display_animation_fps": (0,),  # int>0
@@ -90,8 +88,6 @@ def test_defaults_bad_inputs(key, value, expected_errortype):
 # over the values for each key
 good_inputs = {
     "checkinputs": (True, False),  # bool
-    "edgesize": (1e-9, 2),  # float>0
-    "itercylinder": (10,),  # int>0
     "display_autosizefactor": (1,),  # float>0
     "display_animation_maxfps": (10,),  # int>0
     "display_animation_fps": (10,),  # int>0
@@ -204,7 +200,7 @@ def test_bad_default_classes():
 
 def test_resetting_defaults():
     """test setting and resetting the config"""
-    magpy.defaults.itercylinder = 15
-    assert magpy.defaults.itercylinder == 15, "setting config failed"
+    magpy.defaults.checkinputs = False
+    assert magpy.defaults.checkinputs == False, "setting config failed"
     magpy.defaults.reset()
-    assert magpy.defaults.itercylinder == 50, "resetting config failed"
+    assert magpy.defaults.checkinputs == True, "resetting config failed"
diff --git a/tests/test_exceptions.py b/tests/test_exceptions.py
index 184b88d7c..35c34b916 100644
--- a/tests/test_exceptions.py
+++ b/tests/test_exceptions.py
@@ -9,24 +9,24 @@
     MagpylibBadInputShape)
 from magpylib._src.utility import format_obj_input, format_src_inputs, format_obs_inputs
 from magpylib._src.utility import test_path_format as tpf
-from magpylib._src.fields.field_wrap_BH_level2_dict import getB_dict
 
 
 def getBHv_unknown_source_type():
     """ unknown source type """
-    getB_dict(
+    getBH_dict_level2(
         source_type='badName',
         magnetization=(1,0,0),
         dimension=(0,2,1,0,360),
         position=(0,0,-.5),
-        observer=(1.5,0,-.1))
+        observer=(1.5,0,-.1),
+        field='B')
 
 def getBH_level1_internal_error():
     """ bad source_type input should not happen
     """
     x = np.array([(1,2,3)])
     rot = R.from_quat((0,0,0,1))
-    getBH_level1(bh=True,source_type='woot', magnetization=x, dimension=x, observer=x,
+    getBH_level1(field='B',source_type='woot', magnetization=x, dimension=x, observer=x,
         position=x, orientation=rot)
 
 
@@ -35,7 +35,7 @@ def getBH_level2_bad_input1():
     """
     src = magpy.magnet.Cuboid((1,1,2),(1,1,1))
     sens = magpy.Sensor()
-    getBH_level2(True, [src,sens], (0,0,0), False, True)
+    getBH_level2([src,sens], (0,0,0), sumup=False, squeeze=True, field='B')
 
 
 def getBH_level2_bad_input2():
@@ -244,7 +244,7 @@ def bad_input_shape_dipole_mom():
     magpy.misc.Dipole(moment=vec4)
 
 
-########################################################################
+#####################################################################
 class TestExceptions(unittest.TestCase):
     """ test class for exception testing
     """
diff --git a/tests/test_field_cylinder.py b/tests/test_field_cylinder.py
index 2182c17fd..c850f839a 100644
--- a/tests/test_field_cylinder.py
+++ b/tests/test_field_cylinder.py
@@ -4,9 +4,9 @@
 
 import pytest
 import numpy as np
-from magpylib._src.fields.field_BH_cylinder_tile import (
-    magnet_cyl_tile_H_Slanovc2021,
-    field_BH_cylinder_tile,
+from magpylib._src.fields.field_BH_cylinder_segment import (
+    magnet_cylinder_segment_core,
+    magnet_cylinder_segment_field,
 )
 import magpylib as magpy
 
@@ -14,7 +14,7 @@
 # creating test data
 """ import os
 import numpy as np
-from magpylib._src.fields.field_BH_cylinder_tile import magnet_cyl_tile_H_Slanovc2021
+from magpylib._src.fields.field_BH_cylinder_tile import magnet_cylinder_segment_core
 
 N = 1111
 null = np.zeros(N)
@@ -40,7 +40,7 @@
 r = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [112, 212, 132, 232]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -54,7 +54,7 @@
 r = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [122, 222, 132, 232]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -68,7 +68,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [113, 213, 133, 233, 115, 215, 135, 235]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -83,7 +83,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [123, 223, 133, 233, 125, 225, 135, 235]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -97,7 +97,7 @@
 r = r2
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [125, 225, 135, 235, 124, 224, 134, 234]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -112,7 +112,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [211, 221, 212, 222]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -126,7 +126,7 @@
 r = r1
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [214, 224, 215, 225]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -142,7 +142,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [111, 211, 121, 221, 112, 212, 122, 222]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -157,7 +157,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [111, 211, 131, 231, 112, 212, 132, 232]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -171,7 +171,7 @@
 r = r2
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cyl_tile_H_Slanovc2021(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
 DATA["cases [115, 215, 135, 235, 114, 214, 134, 234]"] = {
     "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
     "H_expected": H1,
@@ -193,7 +193,7 @@
 )
 def test_cylinder_tile_slanovc(inputs, H_expected):
     "testing precomputed cylinder test cases"
-    H = magnet_cyl_tile_H_Slanovc2021(**inputs)
+    H = magnet_cylinder_segment_core(**inputs)
     assert np.allclose(np.nan_to_num(H), np.nan_to_num(H_expected))
 
 
@@ -214,15 +214,15 @@ def test_cylinder_field1():
     implementations
     """
     N = 100
-    magg, dim, poso, H0 = np.load("tests/testdata/testdata_full_cyl.npy")
+    magg, dim, poso, B0 = np.load("tests/testdata/testdata_full_cyl.npy")
 
     nulll = np.zeros(N)
     eins = np.ones(N)
     d, h, _ = dim.T
-    dim6 = np.array([nulll, d / 2, nulll, eins * 360, -h / 2, h / 2]).T
-    H1 = field_BH_cylinder_tile(True, magg, dim6, poso)
+    dim5 = np.array([nulll, d / 2, h, nulll, eins * 360]).T
+    B1 = magnet_cylinder_segment_field(magg, dim5, poso)
 
-    assert np.allclose(H1, H0)
+    assert np.allclose(B1, B0)
 
 
 def test_cylinder_slanovc_field2():
@@ -389,3 +389,74 @@ def test_cylinder_tile_vs_fem():
     assert np.amax((fd2[5:-5, 1:] * 1000 - B2[5:-5]).T / amp2[5:-5]) < 0.05
     assert np.amax((fd3[:, 1:] * 1000 - B3).T / amp3) < 0.05
     assert np.amax((fd4[:, 1:] * 1000 - B4).T / amp4) < 0.05
+
+
+def test_cylinder_corner():
+    """test corner =0 behavior"""
+    a = 1
+    s = magpy.magnet.Cylinder((10,10,1000), (2*a,2*a))
+    B = s.getB([[0,a,a],[0,a,-a],[0,-a,-a],[0,-a,a],
+        [a,0,a],[a,0,-a],[-a,0,-a],[-a,0,a]])
+    np.testing.assert_allclose(B, np.zeros((8,3)))
+
+
+def test_cylinder_corner_scaling():
+    """ test corner=0 scaling"""
+    a = 1
+    obs = [[a,0,a+1e-14], [a+1e-14,0,a]]
+    s = magpy.magnet.Cylinder((10,10,1000), (2*a,2*a))
+    Btest = [[5.12553286e+03, -2.26623480e+00, 2.59910242e+02],
+            [5.12803286e+03, -2.26623480e+00, 9.91024238e+00]]
+    np.testing.assert_allclose(s.getB(obs), Btest)
+
+    a = 1000
+    obs = [[a,0,a+1e-14], [a+1e-14,0,a]]
+    s = magpy.magnet.Cylinder((10,10,1000), (2*a,2*a))
+    np.testing.assert_allclose(s.getB(obs), np.zeros((2,3)))
+
+
+def test_cylinder_scaling_invariance():
+    """test scaling invariance"""
+    obs = np.array([
+        [-0.12788963,  0.14872334, -0.35838915],
+        [-0.17319799,  0.39177646,  0.22413971],
+        [-0.15831916, -0.39768996,  0.41800279],
+        [-0.05762575,  0.19985373,  0.02645361],
+        [ 0.19120126, -0.13021813, -0.21615004],
+        [ 0.39272212,  0.36457661, -0.09758084],
+        [-0.39270581, -0.19805643,  0.36988649],
+        [ 0.28942161,  0.31003054, -0.29558298],
+        [ 0.13083584,  0.31396182, -0.11231319],
+        [-0.04097917,  0.43394138, -0.14109254]])
+
+    a = 1e-6
+    s1 = magpy.magnet.Cylinder((10,10,1000), (2*a,2*a))
+    Btest1 = s1.getB(obs*a)
+
+    a = 1
+    s2 = magpy.magnet.Cylinder((10,10,1000), (2*a,2*a))
+    Btest2 = s2.getB(obs)
+
+    a = 1e7
+    s3 = magpy.magnet.Cylinder((10,10,1000), (2*a,2*a))
+    Btest3 = s3.getB(obs*a)
+
+    np.testing.assert_allclose(Btest1, Btest2)
+    np.testing.assert_allclose(Btest1, Btest3)
+
+def test_cylinder_diametral_small_r():
+    """
+    test if the transition from Taylor series to general case is smooth
+    test if the gneral case fluctuations are small
+    """
+    B = magpy.core.magnet_cylinder_field(
+        np.array([(1,1,0)]*1000),
+        np.array([(2,2)]*1000),
+        np.array([(x,0,3) for x in np.logspace(-1.4,-1.2,1000)]),
+        field='B')
+
+    dB = np.log(abs(B[1:]-B[:-1]))
+    ddB = abs(dB[1:]-dB[:-1])
+    ddB = abs(ddB - np.mean(ddB, axis=0))
+
+    assert np.all(ddB < 0.001)
diff --git a/tests/test_field_functions.py b/tests/test_field_functions.py
index c54ac1b82..f02d06f6c 100644
--- a/tests/test_field_functions.py
+++ b/tests/test_field_functions.py
@@ -1,73 +1,44 @@
-import pickle
-import os
 import numpy as np
-from magpylib._src.fields.field_BH_cuboid import field_BH_cuboid
-from magpylib._src.fields.field_BH_cylinder_tile import field_BH_cylinder_tile
-from magpylib._src.fields.field_BH_sphere import field_BH_sphere
-from magpylib._src.fields.field_BH_dipole import field_BH_dipole
-from magpylib._src.fields.field_BH_loop import field_BH_loop
-from magpylib._src.fields.field_BH_line import field_BH_line, field_BH_line_from_vert
-from magpylib._src.defaults.defaults_classes import default_settings as Config
-
-# # GENERATE TEST DATA
-# n = 500
-# magpy.defaults.edgesize = 1e-14
-
-# # dim general
-# dim_gen = np.random.rand(n,3)
-# # dim edge
-# dim_edge = np.array([(2,2,2)]*n)
-
-# # pure edge positions
-# pos_edge = (np.random.rand(n,3)-0.5)*2
-# pos_edge[:,0]=1 + (np.random.rand(n)-0.5)*1e-14
-# pos_edge[:,1]=1  + (np.random.rand(n)-0.5)*1e-14
-# for i in range(n):
-#     np.random.shuffle(pos_edge[i])
-# # general positions
-# pos_gen = (np.random.rand(n,3)-0.5)*.5 # half inner half outer
-# # mixed positions
-# aa = np.r_[pos_edge, pos_gen]
-# np.random.shuffle(aa)
-# pos_mix = aa[:n]
-
-# # general mag (no special cases at this point)
-# mag = (np.random.rand(n,3)-0.5)*1000
-
-# poss = [pos_edge, pos_gen, pos_mix]
-# dims = [dim_edge, dim_gen, dim_edge]
-
-# Bs = []
-# for dim,pos in zip(dims,poss):
-#     Bs += [field_BH_cuboid(True, mag, dim, pos)]
-# Bs = np.array(Bs)
-
-# print(Bs)
-# pickle.dump([mag, dims, poss, Bs], open('testdata_field_BH_cuboid.p','wb'))
-
-
-def test_field_BH_cuboid():
-    """ test cuboid field
-    """
-    Config.edgesize=1e-14
-    mag, dims, poss, B = pickle.load(open(
-        os.path.abspath('tests/testdata/testdata_field_BH_cuboid.p') ,'rb'))
-    Btest = []
-    for dim,pos in zip(dims,poss):
-        Btest += [field_BH_cuboid(True, mag, dim, pos)]
-    Btest = np.array(Btest)
-    assert np.allclose(B, Btest), 'Cuboid field computation broken'
+from numpy.testing import assert_allclose
+from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
+from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_field
+from magpylib._src.fields.field_BH_sphere import magnet_sphere_field
+from magpylib._src.fields.field_BH_dipole import dipole_field
+from magpylib._src.fields.field_BH_loop import current_loop_field
+from magpylib._src.fields.field_BH_line import current_line_field, field_BH_line_from_vert
 
 
-def test_field_BH_cuboid_mag0():
+def test_magnet_cuboid_Bfield():
+    """ test cuboid field
+    """
+    mag = np.array([(0,0,0), (1,2,3), (1,2,3), (1,2,3), (1,2,3), (2,2,2), (2,2,2), (1,1,1), (1,1,1)])
+    dim = np.array([(1,2,3), (-1,-2,2), (1,2,2), (0,2,2), (1,2,3), (2,2,2), (2,2,2), (2,2,2), (2,2,2)])
+    pos = np.array([(1,2,3), (1,-1,0), (1,-1,0), (1,-1,0), (1,2,3), (1,1+1e-14,0), (1,1,1), (1,-1,2), (1+1e-14,-1,2)])
+    B = magnet_cuboid_field(mag, dim, pos)
+
+    Btest = [
+    [ 0.         , 0.         , 0.        ],
+    [-0.14174376 ,-0.16976459 ,-0.20427478],
+    [-0.14174376 ,-0.16976459 ,-0.20427478],
+    [ 0.         , 0.         , 0.        ],
+    [ 0.02596336 , 0.04530334 , 0.05840059],
+    [     np.inf ,     np.inf ,-0.29516724],
+    [ 0.         , 0.         , 0.        ],
+    [-0.0009913  ,-0.08747071 , 0.04890262],
+    [-0.0009913  ,-0.08747071 , 0.04890262],]
+
+    np.testing.assert_allclose(B, Btest, rtol=1e-5)
+
+
+def test_magnet_cuboid_field_mag0():
     """ test cuboid field magnetization=0
     """
     n = 10
     mag = np.zeros((n,3))
     dim = np.random.rand(n,3)
     pos = np.random.rand(n,3)
-    B = field_BH_cuboid(True, mag, dim, pos)
-    assert np.allclose(mag,B), 'Cuboid magnetization=0 case broken'
+    B = magnet_cuboid_field(mag, dim, pos)
+    assert_allclose(mag,B)
 
 
 def test_field_BH_cylinder_tile_mag0():
@@ -75,14 +46,13 @@ def test_field_BH_cylinder_tile_mag0():
     """
     n = 10
     mag = np.zeros((n,3))
-    r1,r2,phi1,phi2,z1,z2 = np.random.rand(6,n)
+    r1,r2,h,phi1,phi2 = np.random.rand(5,n)
     r2=r1+r2
     phi2=phi1+phi2
-    z2=z1+z2
-    dim = np.array([r1,r2,phi1,phi2,z1,z2]).T
+    dim = np.array([r1,r2,h,phi1,phi2]).T
     pos = np.random.rand(n,3)
-    B = field_BH_cylinder_tile(True, mag, dim, pos)
-    assert np.allclose(mag,B), 'CylinderTile magnetization=0 case broken'
+    B = magnet_cylinder_segment_field(mag, dim, pos)
+    assert_allclose(mag, B)
 
 
 def test_field_sphere_vs_v2():
@@ -101,20 +71,20 @@ def test_field_sphere_vs_v2():
     dim = np.array([1.23]*7)
     mag = np.array([(33,66,99)]*7)
     poso = np.array([(0,0,0),(.2,.2,.2),(.4,.4,.4),(-1,-1,-2),(.1,.1,.1),(1,2,-3),(-3,2,1)])
-    B = field_BH_sphere(True, mag, dim, poso )
+    B = magnet_sphere_field(mag, dim, poso )
 
-    assert np.allclose(result_v2, B), 'vs_v2 failed'
+    np.testing.assert_allclose(result_v2, B, rtol=1e-6)
 
 
-def test_field_BH_sphere_mag0():
+def test_magnet_sphere_field_mag0():
     """ test cuboid field magnetization=0
     """
     n = 10
     mag = np.zeros((n,3))
     dim = np.random.rand(n)
     pos = np.random.rand(n,3)
-    B = field_BH_sphere(True, mag, dim, pos)
-    assert np.allclose(mag,B), 'Cuboid magnetization=0 case broken'
+    B = magnet_sphere_field(mag, dim, pos)
+    assert_allclose(mag,B)
 
 
 def test_field_dipole1():
@@ -122,28 +92,31 @@ def test_field_dipole1():
     """
     poso = np.array([(1,2,3),(-1,2,3)])
     mom = np.array([(2,3,4),(0,-3,-2)])
-    B = field_BH_dipole(True, mom, poso)*np.pi
+    B = dipole_field(mom, poso, field='B')*np.pi
     Btest = np.array([
         (0.01090862,0.02658977,0.04227091),
         (0.0122722,-0.01022683,-0.02727156),
         ])
 
-    assert np.allclose(B,Btest)
+    assert_allclose(B, Btest, rtol=1e-6)
 
 
 def test_field_dipole2():
     """ test nan return when pos_obs=0
     """
-    poso = np.array([(0,0,0)])
-    mom = np.array([(-1,2,3)])
-    np.seterr(all='ignore')
-    B = field_BH_dipole(True, mom, poso)*np.pi
-    np.seterr(all='warn')
+    moment = np.array([(100,200,300)]*2 + [(0,0,0)]*2)
+    observer = np.array([(0,0,0),(1,2,3)]*2)
+    B = dipole_field(moment, observer)
+
+    assert all(np.isinf(B[0]))
+    assert_allclose(B[1:],
+        [[0.3038282,  0.6076564,  0.91148459],
+         [0.,         0.,         0.        ],
+         [0.,         0.,         0.        ]])
 
-    assert all(np.isnan(B[0]))
 
 
-def test_field_circular():
+def test_field_loop():
     """ test if field function gives correct outputs
     """
     # from hyperphysics
@@ -166,36 +139,48 @@ def test_field_circular():
         [-0.27939151,  0.27939151,  0.01220605],
         [ 3.25697271,  0.32569727, -5.49353046]]
 
-    pos_test = np.array(pos_test_hyper + pos_test_mag2 + [[1,0,0]])
-    Btest = np.array(Btest_hyper + Btest_mag2 + [[0,0,0]])
+    pos_test = np.array(pos_test_hyper + pos_test_mag2)
+    Btest = np.array(Btest_hyper + Btest_mag2)
 
-    current = np.array([1,1] + [123]*8 + [123])
-    dim = np.array([2,2] + [2]*8 + [2])
+    current = np.array([1,1] + [123]*8)
+    dim = np.array([2,2] + [2]*8)
 
-    B = field_BH_loop(True, current, dim, pos_test)
+    B = current_loop_field(current, dim, pos_test)
 
-    assert np.allclose(B, Btest)
+    assert_allclose(B, Btest, rtol=1e-6)
 
     Htest = Btest*10/4/np.pi
-    H = field_BH_loop(False, current, dim, pos_test)
-    assert np.allclose(H, Htest)
+    H = current_loop_field(current, dim, pos_test, field='H')
+    assert_allclose(H, Htest, rtol=1e-6)
 
 
-def test_field_circular2():
+def test_field_loop2():
     """ test if field function accepts correct inputs
     """
     curr = np.array([1])
     dim = np.array([2])
     poso = np.array([[0,0,0]])
-    B = field_BH_loop(True, curr, dim, poso)
+    B = current_loop_field(curr, dim, poso)
 
     curr = np.array([1]*2)
     dim = np.array([2]*2)
     poso = np.array([[0,0,0]]*2)
-    B2 = field_BH_loop(True, curr, dim, poso)
+    B2 = current_loop_field(curr, dim, poso)
+
+    assert_allclose(B, (B2[0],))
+    assert_allclose(B, (B2[1],))
 
-    assert np.allclose(B,B2[0])
-    assert np.allclose(B,B2[1])
+
+def test_field_loop_specials():
+    """ test loop special cases
+    """
+    cur = np.array([1,1,1,1,0,2])
+    dia = np.array([2,2,0,0,2,2])
+    obs = np.array([(0,0,0), (1,0,0), (0,0,0), (1,0,0), (1,0,0), (0,0,0)])
+
+    B = current_loop_field(cur, dia, obs)
+    Btest = [[0,0,0.62831853], [0,0,0], [0,0,0], [0,0,0], [0,0,0], [0,0,1.25663706]]
+    assert_allclose(B, Btest)
 
 
 def test_field_line():
@@ -208,50 +193,50 @@ def test_field_line():
     pe1 = np.array([(2,2,2)])
 
     # only normal
-    B1 = field_BH_line(True, c1, ps1, pe1, po1)
+    B1 = current_line_field(c1, ps1, pe1, po1)
     x1 = np.array([[ 0.02672612, -0.05345225, 0.02672612]])
-    assert np.allclose(x1, B1)
+    assert_allclose(x1, B1, rtol=1e-6)
 
     # only on_line
     po1b = np.array([(1,1,1)])
-    B2 = field_BH_line(True, c1, ps1, pe1, po1b)
+    B2 = current_line_field(c1, ps1, pe1, po1b)
     x2 = np.zeros((1,3))
-    assert np.allclose(x2, B2)
+    assert_allclose(x2, B2, rtol=1e-6)
 
     # only zero-segment
-    B3 = field_BH_line(True, c1, ps1, ps1, po1)
+    B3 = current_line_field(c1, ps1, ps1, po1)
     x3 = np.zeros((1,3))
-    assert np.allclose(x3, B3)
+    assert_allclose(x3, B3, rtol=1e-6)
 
     # only on_line and zero_segment
     c2 = np.array([1]*2)
     ps2 = np.array([(0,0,0)]*2)
     pe2 = np.array([(0,0,0),(2,2,2)])
     po2 = np.array([(1,2,3),(1,1,1)])
-    B4 = field_BH_line(True, c2, ps2, pe2, po2)
+    B4 = current_line_field(c2, ps2, pe2, po2)
     x4 = np.zeros((2,3))
-    assert np.allclose(x4, B4)
+    assert_allclose(x4, B4, rtol=1e-6)
 
     # normal + zero_segment
     po2b = np.array([(1,2,3),(1,2,3)])
-    B5 = field_BH_line(True, c2, ps2, pe2, po2b)
+    B5 = current_line_field(c2, ps2, pe2, po2b)
     x5 = np.array([[0,0,0],[ 0.02672612, -0.05345225, 0.02672612]])
-    assert np.allclose(x5, B5)
+    assert_allclose(x5, B5, rtol=1e-6)
 
     # normal + on_line
     pe2b = np.array([(2,2,2)]*2)
-    B6 = field_BH_line(True, c2, ps2, pe2b, po2)
+    B6 = current_line_field(c2, ps2, pe2b, po2)
     x6 = np.array([[0.02672612, -0.05345225, 0.02672612],[0,0,0]])
-    assert np.allclose(x6, B6)
+    assert_allclose(x6, B6, rtol=1e-6)
 
     # normal + zero_segment + on_line
     c4 = np.array([1]*3)
     ps4 = np.array([(0,0,0)]*3)
     pe4 = np.array([(0,0,0),(2,2,2),(2,2,2)])
     po4 = np.array([(1,2,3),(1,2,3),(1,1,1)])
-    B7 = field_BH_line(True, c4, ps4, pe4, po4)
+    B7 = current_line_field(c4, ps4, pe4, po4)
     x7 = np.array([[0,0,0], [0.02672612, -0.05345225, 0.02672612], [0,0,0]])
-    assert np.allclose(x7, B7)
+    assert_allclose(x7, B7, rtol=1e-6)
 
 
 def test_field_line_from_vert():
@@ -265,7 +250,7 @@ def test_field_line_from_vert():
     vert3 = np.array([(1,2,3),(-2,-3,3),(3,2,1),(3,3,3)])
 
     pos_tiled = np.tile(p, (3,1))
-    B_vert = field_BH_line_from_vert(True, curr, [vert1,vert2,vert3], pos_tiled)
+    B_vert = field_BH_line_from_vert(curr, [vert1,vert2,vert3], pos_tiled, field='B')
 
     B = []
     for i,vert in enumerate([vert1,vert2,vert3]):
@@ -274,7 +259,26 @@ def test_field_line_from_vert():
             p2 = vert[1:]
             po = np.array([pos]*(len(vert)-1))
             cu = np.array([curr[i]]*(len(vert)-1))
-            B += [np.sum(field_BH_line(True, cu, p1, p2, po), axis=0)]
+            B += [np.sum(current_line_field(cu, p1, p2, po), axis=0)]
     B = np.array(B)
 
-    assert np.allclose(B_vert, B)
+    assert_allclose(B_vert, B)
+
+
+def test_field_line_v4():
+    """ test current_line_Bfield() for all cases
+    """
+    cur = np.array([1]*7)
+    start = np.array([(-1,0,0)]*7)
+    end = np.array([(1,0,0), (-1,0,0), (1,0,0), (-1,0,0)] + [(1,0,0)]*3)
+    obs = np.array([(0,0,1),(0,0,0), (0,0,0), (0,0,0), (0,0,1e-16), (2,0,1), (-2,0,1)])
+    B = current_line_field(cur, start, end, obs)
+    Btest = np.array(
+        [[0, -0.14142136, 0],
+         [0,  0.        , 0],
+         [0,  0.        , 0],
+         [0,  0.        , 0],
+         [0,  0.        , 0],
+         [0, -0.02415765, 0],
+         [0, -0.02415765, 0]])
+    np.testing.assert_allclose(B, Btest)
diff --git a/tests/test_getBH_dict.py b/tests/test_getBH_dict.py
index 9c77d78bd..800fd8288 100644
--- a/tests/test_getBH_dict.py
+++ b/tests/test_getBH_dict.py
@@ -4,9 +4,9 @@
 import magpylib
 from magpylib.magnet import Cuboid, Cylinder, Sphere
 from magpylib import getB, getH
-from magpylib._src.fields.field_wrap_BH_level2_dict import getB_dict, getH_dict
 from magpylib._src.exceptions import MagpylibBadUserInput
 
+
 def test_getB_dict1():
     """test field wrapper functions
     """
@@ -24,14 +24,12 @@ def test_getB_dict1():
     rot = pm.orientation
 
     dic = {
-        'source_type': 'Cylinder',
-        'observer': pos_obs,
         'magnetization': mag,
         'dimension': dim,
         'position': pos,
         'orientation':rot
         }
-    B1 = getB_dict(**dic)
+    B1 = getB('Cylinder', pos_obs, **dic)
 
     assert np.allclose(B1, B2, rtol=1e-12, atol=1e-12)
 
@@ -45,13 +43,11 @@ def test_getB_dict2():
     pos = [(1,1,1),(2,2,2),(3,3,3),(5,5,5)]
 
     dic = {
-        'source_type': 'Cylinder',
-        'observer': pos_obs,
         'magnetization': mag,
         'dimension': dim,
         'position': pos
         }
-    B1 = getB_dict(**dic)
+    B1 = getB('Cylinder', pos_obs, **dic)
 
     pm = Cylinder(mag, dim, position=pos)
     B2 = getB([pm],pos_obs)
@@ -67,12 +63,10 @@ def test_getH_dict1():
     dim = [3,3]
 
     dic = {
-        'source_type': 'Cylinder',
-        'observer': pos_obs,
         'magnetization': mag,
         'dimension': dim,
         }
-    B1 = getH_dict(**dic)
+    B1 = getH('Cylinder', pos_obs, **dic)
 
     pm = Cylinder(mag, dim)
     B2 = pm.getH(pos_obs)
@@ -91,14 +85,12 @@ def test_getB_dict3():
     rot = R.from_quat([(t,.2,.3,.4) for t in np.linspace(0,.1,n)])
 
     dic = {
-        'source_type': 'Cuboid',
-        'observer': pos_obs,
         'magnetization': mag,
         'dimension': dim,
         'position': pos,
         'orientation': rot
         }
-    B1 = getB_dict(**dic)
+    B1 = getB('Cuboid', pos_obs, **dic)
 
     B2 = []
     for i in range(n):
@@ -117,12 +109,10 @@ def test_getH_dict3():
     dim = 3
 
     dic = {
-        'source_type': 'Sphere',
-        'observer': pos_obs,
         'magnetization': mag,
         'diameter': dim
         }
-    B1 = getH_dict(**dic)
+    B1 = getH('Sphere', pos_obs, **dic)
 
     B2 = []
     for i in range(3):
@@ -144,14 +134,12 @@ def test_getB_dict4():
     rot = R.from_quat([(t,.2,.3,.4) for t in np.linspace(0,.1,n)])
 
     dic = {
-        'source_type': 'Sphere',
-        'observer': pos_obs,
         'magnetization': mag,
         'diameter': dim,
         'position': pos,
         'orientation': rot
         }
-    B1 = getB_dict(**dic)
+    B1 = getB('Sphere', pos_obs, **dic)
 
     B2 = []
     for i in range(n):
@@ -165,11 +153,11 @@ def test_getB_dict4():
 def test_geBHv_dipole():
     """ test if Dipole implementation gives correct output
     """
-    B = getB_dict(source_type='Dipole', moment=(1,2,3), observer=(1,1,1))
+    B = getB('Dipole', (1,1,1), moment=(1,2,3))
     Btest = np.array([0.07657346,0.06125877,0.04594407])
     assert np.allclose(B,Btest)
 
-    H = getH_dict(source_type='Dipole', moment=(1,2,3), observer=(1,1,1))
+    H = getH('Dipole', (1,1,1), moment=(1,2,3))
     Htest = np.array([0.06093522,0.04874818,0.03656113])
     assert np.allclose(H,Htest)
 
@@ -177,11 +165,11 @@ def test_geBHv_dipole():
 def test_geBHv_circular():
     """ test if Loop implementation gives correct output
     """
-    B = getB_dict(source_type='Loop', current=1, diameter=2, observer=(0,0,0))
+    B = getB('Loop', (0,0,0), current=1, diameter=2)
     Btest = np.array([0,0,0.6283185307179586])
     assert np.allclose(B,Btest)
 
-    H = getH_dict(source_type='Loop', current=1, diameter=2, observer=(0,0,0))
+    H = getH('Loop', (0,0,0), current=1, diameter=2)
     Htest = np.array([0,0,0.6283185307179586*10/4/np.pi])
     assert np.allclose(H,Htest)
 
@@ -189,10 +177,10 @@ def test_geBHv_circular():
 def test_getBHv_squeeze():
     """ test if squeeze works
     """
-    B1 = getB_dict(source_type='Loop', current=1, diameter=2, observer=(0,0,0))
-    B2 = getB_dict(source_type='Loop', current=1, diameter=2, observer=[(0,0,0)])
-    B3 = getB_dict(source_type='Loop', current=1, diameter=2, observer=[(0,0,0)], squeeze=False)
-    B4 = getB_dict(source_type='Loop', current=1, diameter=2, observer=[(0,0,0)]*2)
+    B1 = getB('Loop', (0,0,0)    , current=1, diameter=2)
+    B2 = getB('Loop', [(0,0,0)]  , current=1, diameter=2)
+    B3 = getB('Loop', [(0,0,0)]  , current=1, diameter=2, squeeze=False)
+    B4 = getB('Loop', [(0,0,0)]*2, current=1, diameter=2)
 
     assert B1.ndim == 1
     assert B2.ndim == 1
@@ -203,9 +191,9 @@ def test_getBHv_squeeze():
 def test_getBHv_line():
     """ test getBHv with Line
     """
-    H = getH_dict(
-        source_type='Line',
-        observer=[(1,1,1),(1,2,3),(2,2,2)],
+    H = getH(
+        'Line',
+        [(1,1,1),(1,2,3),(2,2,2)],
         current=1,
         segment_start=(0,0,0),
         segment_end=[(0,0,0),(2,2,2),(2,2,2)])
@@ -219,30 +207,30 @@ def test_getBHv_line2():
     x = 0.14142136
 
     # z-line on x=1
-    B1 = getB_dict(source_type='Line', observer=(0,0,0), current=1,
+    B1 = getB('Line', (0,0,0), current=1,
         segment_start=(1,0,-1), segment_end=(1,0,1))
     assert np.allclose(B1, np.array([0,-x,0]))
 
     # move z-line to x=-1
-    B2 = getB_dict(source_type='Line', position=(-2,0,0), observer=(0,0,0), current=1,
+    B2 = getB('Line', (0,0,0), position=(-2,0,0), current=1,
         segment_start=(1,0,-1), segment_end=(1,0,1))
     assert np.allclose(B2, np.array([0,x,0]))
 
     # rotate 1
     rot = R.from_euler('z', 90, degrees=True)
-    B3 = getB_dict(source_type='Line', orientation=rot, observer=(0,0,0), current=1,
+    B3 = getB('Line', (0,0,0), orientation=rot, current=1,
         segment_start=(1,0,-1), segment_end=(1,0,1))
     assert np.allclose(B3, np.array([x,0,0]))
 
     # rotate 2
     rot = R.from_euler('x', 90, degrees=True)
-    B4 = getB_dict(source_type='Line', orientation=rot, observer=(0,0,0), current=1,
+    B4 = getB('Line', (0,0,0), orientation=rot, current=1,
         segment_start=(1,0,-1), segment_end=(1,0,1))
     assert np.allclose(B4, np.array([0,0,-x]))
 
     # rotate 3
     rot = R.from_euler('y', 90, degrees=True)
-    B5 = getB_dict(source_type='Line', orientation=rot, observer=(0,0,0), current=1,
+    B5 = getB('Line', (0,0,0), orientation=rot, current=1,
         segment_start=(1,0,-1), segment_end=(1,0,1))
     assert np.allclose(B5, np.array([0,-x,0]))
 
@@ -286,12 +274,13 @@ def test_BHv_Cylinder_FEM():
         (0.0121250819870128,0.0104894041620816,0.00303690098080925)])
 
     # compare against FEM
-    B = getB_dict(
-        source_type='CylinderSegment',
+    B = getB(
+        'CylinderSegment',
+        obsp,
         dimension=(1,2,1,90,360),
         magnetization=np.array((1,2,3))*1000/np.sqrt(14),
         position = (0,0,0.5),
-        observer=obsp)
+    )
 
     err = np.linalg.norm(B-Bfem*1000, axis=1)/np.linalg.norm(B,axis=1)
     assert np.amax(err)<0.01
@@ -300,26 +289,29 @@ def test_BHv_Cylinder_FEM():
 def test_BHv_solid_cylinder():
     """compare multiple solid-cylinder solutions against each other"""
     # combine multiple slices to one big Cylinder
-    B1 = getB_dict(
-        source_type='CylinderSegment',
+    B1 = getB(
+        'CylinderSegment',
+        (1,2,3),
         dimension=[(0,1,2,20,120), (0,1,2,120,220), (0,1,2,220,380)],
         magnetization=(22,33,44),
-        observer=(1,2,3))
+    )
     B1 = np.sum(B1,axis=0)
 
     # one big cylinder
-    B2 = getB_dict(
-        source_type='CylinderSegment',
+    B2 = getB(
+        'CylinderSegment',
+        (1,2,3),
         dimension=(0,1,2,0,360),
         magnetization=(22,33,44),
-        observer=(1,2,3))
+    )
 
     # compute with solid cylinder code
-    B3 = getB_dict(
-        source_type='Cylinder',
+    B3 = getB(
+        'Cylinder',
+        (1,2,3),
         dimension=(2,2),
         magnetization=(22,33,44),
-        observer=(1,2,3))
+    )
 
     assert np.allclose(B1,B2)
     assert np.allclose(B1,B3)
diff --git a/tests/test_getBH_interfaces.py b/tests/test_getBH_interfaces.py
index be0bad83e..63c48e980 100644
--- a/tests/test_getBH_interfaces.py
+++ b/tests/test_getBH_interfaces.py
@@ -1,6 +1,5 @@
 import numpy as np
 import magpylib as magpy
-from magpylib._src.fields.field_wrap_BH_level2_dict import getB_dict, getH_dict
 
 
 def test_getB_interfaces1():
@@ -10,30 +9,26 @@ def test_getB_interfaces1():
     src.move(np.linspace((.1,.2,.3), (1,2,3), 10), start=-1)
     poso = [[(-1,-1,-1)]*2]*2
     sens = magpy.Sensor(pixel=poso)
-    B = getB_dict(
-        source_type='Cuboid',
+    B = magpy.getB(
+        'Cuboid',
+        (-1,-1,-1),
         position=src.position,
         magnetization=(1,2,3),
         dimension=(1,2,3),
-        observer=(-1,-1,-1))
+        )
     B1 = np.tile(B,(2,2,1,1))
     B1 = np.swapaxes(B1,0,2)
 
     B_test = magpy.getB(src, sens)
-    assert B_test.shape == B1.shape, 'magpy.getB1 shape failed'
-    assert np.allclose(B1, B_test), 'magpy.getB1 value failed'
+    np.testing.assert_allclose(B1, B_test)
 
     B_test = src.getB(poso)
-    assert B_test.shape == B1.shape, 'src.getB shape failed'
-    assert np.allclose(B1, B_test), 'src.getB value failed'
+    np.testing.assert_allclose(B1, B_test)
 
     B_test = src.getB(sens)
-    assert B_test.shape == B1.shape, 'src.getB3 shape failed'
-    assert np.allclose(B1, B_test), 'src.getB3 value failed'
+    np.testing.assert_allclose(B1, B_test)
 
     B_test = sens.getB(src)
-    assert B_test.shape == B1.shape, 'sens.getB1 shape failed'
-    assert np.allclose(B1, B_test), 'sens.getB1 value failed'
 
 
 def test_getB_interfaces2():
@@ -43,23 +38,22 @@ def test_getB_interfaces2():
     src.move(np.linspace((.1,.2,.3), (1,2,3), 10), start=-1)
     poso = [[(-1,-1,-1)]*2]*2
     sens = magpy.Sensor(pixel=poso)
-    B = getB_dict(
-        source_type='Cuboid',
+    B = magpy.getB(
+        'Cuboid',
+        (-1,-1,-1),
         position=src.position,
         magnetization=(1,2,3),
         dimension=(1,2,3),
-        observer=(-1,-1,-1))
+        )
 
     B2 = np.tile(B,(2,2,2,1,1))
     B2 = np.swapaxes(B2,1,3)
 
     B_test = magpy.getB([src,src], sens)
-    assert B_test.shape == B2.shape, 'magpy.getB2 shape failed'
-    assert np.allclose(B2, B_test), 'magpy.getB2 value failed'
+    np.testing.assert_allclose(B2, B_test)
 
     B_test = sens.getB([src,src])
-    assert B_test.shape == B2.shape, 'sens.getB2 shape failed'
-    assert np.allclose(B2, B_test), 'sens.getB2 value failed'
+    np.testing.assert_allclose(B2, B_test)
 
 
 def test_getB_interfaces3():
@@ -69,27 +63,24 @@ def test_getB_interfaces3():
     src.move(np.linspace((.1,.2,.3), (1,2,3), 10), start=-1)
     poso = [[(-1,-1,-1)]*2]*2
     sens = magpy.Sensor(pixel=poso)
-    B = getB_dict(
-        source_type='Cuboid',
+    B = magpy.getB(
+        'Cuboid',
+        (-1,-1,-1),
         position=src.position,
         magnetization=(1,2,3),
-        dimension=(1,2,3),
-        observer=(-1,-1,-1))
+        dimension=(1,2,3),)
 
     B3 = np.tile(B,(2,2,2,1,1))
     B3 = np.swapaxes(B3,0,3)
 
     B_test = magpy.getB(src, [sens,sens])
-    assert B_test.shape == B3.shape, 'magpy.getB3 shape failed'
-    assert np.allclose(B3, B_test), 'magpy.getB3 value failed'
+    np.testing.assert_allclose(B3, B_test)
 
     B_test = src.getB([poso,poso])
-    assert B_test.shape == B3.shape, 'src.getB2 shape failed'
-    assert np.allclose(B3, B_test), 'src.getB2 value failed'
+    np.testing.assert_allclose(B3, B_test)
 
     B_test = src.getB([sens,sens])
-    assert B_test.shape == B3.shape, 'src.getB4 shape failed'
-    assert np.allclose(B3, B_test), 'src.getB4 value failed'
+    np.testing.assert_allclose(B3, B_test)
 
 
 def test_getH_interfaces1():
@@ -103,30 +94,27 @@ def test_getH_interfaces1():
     poso = [[(-1,-2,-3)]*2]*2
     sens = magpy.Sensor(pixel=poso)
 
-    H = getH_dict(
-        source_type='Cuboid',
+    H = magpy.getH(
+        'Cuboid',
+        (-1,-2,-3),
         position=src.position,
         magnetization=mag,
         dimension=dim,
-        observer=(-1,-2,-3))
+        )
     H1 = np.tile(H,(2,2,1,1))
     H1 = np.swapaxes(H1,0,2)
 
     H_test = magpy.getH(src, sens)
-    assert H_test.shape == H1.shape, 'magpy.getH1 shape failed'
-    assert np.allclose(H1, H_test), 'magpy.getH1 value failed'
+    np.testing.assert_allclose(H1, H_test)
 
     H_test = src.getH(poso)
-    assert H_test.shape == H1.shape, 'src.getH shape failed'
-    assert np.allclose(H1, H_test), 'src.getH value failed'
+    np.testing.assert_allclose(H1, H_test)
 
     H_test = src.getH(sens)
-    assert H_test.shape == H1.shape, 'src.getH3 shape failed'
-    assert np.allclose(H1, H_test), 'src.getH3 value failed'
+    np.testing.assert_allclose(H1, H_test)
 
     H_test = sens.getH(src)
-    assert H_test.shape == H1.shape, 'sens.getH1 shape failed'
-    assert np.allclose(H1, H_test), 'sens.getH1 value failed'
+    np.testing.assert_allclose(H1, H_test)
 
 
 def test_getH_interfaces2():
@@ -140,23 +128,21 @@ def test_getH_interfaces2():
     poso = [[(-1,-2,-3)]*2]*2
     sens = magpy.Sensor(pixel=poso)
 
-    H = getH_dict(
-        source_type='Cuboid',
+    H = magpy.getH(
+        'Cuboid',
+        (-1,-2,-3),
         position=src.position,
         magnetization=mag,
-        dimension=dim,
-        observer=(-1,-2,-3))
+        dimension=dim,)
 
     H2 = np.tile(H,(2,2,2,1,1))
     H2 = np.swapaxes(H2,1,3)
 
     H_test = magpy.getH([src,src], sens)
-    assert H_test.shape == H2.shape, 'magpy.getH2 shape failed'
-    assert np.allclose(H2, H_test), 'magpy.getH2 value failed'
+    np.testing.assert_allclose(H2, H_test)
 
     H_test = sens.getH([src,src])
-    assert H_test.shape == H2.shape, 'sens.getH2 shape failed'
-    assert np.allclose(H2, H_test), 'sens.getH2 value failed'
+    np.testing.assert_allclose(H2, H_test)
 
 
 def test_getH_interfaces3():
@@ -170,24 +156,22 @@ def test_getH_interfaces3():
     poso = [[(-1,-2,-3)]*2]*2
     sens = magpy.Sensor(pixel=poso)
 
-    H = getH_dict(
-        source_type='Cuboid',
+    H = magpy.getH(
+        'Cuboid',
+        (-1,-2,-3),
         position=src.position,
         magnetization=mag,
         dimension=dim,
-        observer=(-1,-2,-3))
+    )
 
     H3 = np.tile(H,(2,2,2,1,1))
     H3 = np.swapaxes(H3,0,3)
 
     H_test = magpy.getH(src, [sens,sens])
-    assert H_test.shape == H3.shape, 'magpy.getH3 shape failed'
-    assert np.allclose(H3, H_test), 'magpy.getH3 value failed'
+    np.testing.assert_allclose(H3, H_test)
 
     H_test = src.getH([poso,poso])
-    assert H_test.shape == H3.shape, 'src.getH2 shape failed'
-    assert np.allclose(H3, H_test), 'src.getH2 value failed'
+    np.testing.assert_allclose(H3, H_test)
 
     H_test = src.getH([sens,sens])
-    assert H_test.shape == H3.shape, 'src.getH4 shape failed'
-    assert np.allclose(H3, H_test), 'src.getH4 value failed'
+    np.testing.assert_allclose(H3, H_test)
diff --git a/tests/test_input_checks.py b/tests/test_input_checks.py
index 47ab9dd19..0843ff330 100644
--- a/tests/test_input_checks.py
+++ b/tests/test_input_checks.py
@@ -215,16 +215,16 @@ def bad_dim_input9():
     magpy.current.Loop(current=1, diameter=(1,1))
 
 def bad_dim_input10():
-    """cylinder section bad dim shape"""
+    """cylinder segment bad dim shape"""
     magpy.magnet.CylinderSegment(magnetization=a3, dimension=(1,2,0,45))
 def bad_dim_input11():
-    """cylinder section bad dim d1>d2"""
+    """cylinder segment bad dim d1>d2"""
     magpy.magnet.CylinderSegment(magnetization=a3, dimension=(3,2,2,0,45))
 def bad_dim_input12():
-    """cylinder section bad dim phi1>phi2"""
+    """cylinder segment bad dim phi1>phi2"""
     magpy.magnet.CylinderSegment(magnetization=a3, dimension=(1,2,2,100,45))
 def bad_dim_input13():
-    """cylinder section bad dim phi2-phi1>360"""
+    """cylinder segment bad dim phi2-phi1>360"""
     magpy.magnet.CylinderSegment(magnetization=a3, dimension=(1,2,2,0,1145))
 
 
@@ -269,8 +269,66 @@ def bad_observer_input4():
     magpy.getB(src, 123)
 
 
+# test bad level0 inputs
+def bad_level0_field_input_0():
+    """sphere"""
+    mag = np.array([(100,200,300)])
+    dim1 = np.array((2,))
+    obs = np.array([(1,2,3)])
+    magpy.core.magnet_sphere_field(mag, dim1, obs, field='x')
+
+def bad_level0_field_input_1():
+    """cylinder"""
+    mag = np.array([(100,200,300)])
+    dim2 = np.array([(2,2)])
+    obs = np.array([(1,2,3)])
+    magpy.core.magnet_cylinder_field(mag, dim2, obs, field='y')
+
+def bad_level0_field_input_2():
+    """cylinder segment"""
+    mag = np.array([(100,200,300)])
+    dim5 = np.array([(0,1,2,0,360)])
+    obs = np.array([(1,2,3)])
+    magpy.core.magnet_cylinder_segment_field(mag, dim5, obs, field=123)
+
+def bad_level0_field_input_3():
+    """dipole"""
+    mag = np.array([(100,200,300)])
+    obs = np.array([(1,2,3)])
+    magpy.core.dipole_field(mag, obs, field=mag)
+
+def bad_level0_field_input_4():
+    """loop"""
+    cur = np.array((20,))
+    dim1 = np.array((2,))
+    obs = np.array([(1,2,3)])
+    magpy.core.current_loop_field(cur, dim1, obs, field=None)
+
+def bad_level0_field_input_5():
+    """line"""
+    cur = np.array((20,))
+    start = np.array([(1,1,1)])
+    end = np.array([(3,3,3)])
+    obs = np.array([(1,2,3)])
+    magpy.core.current_line_field(cur, start, end, obs, field=(1,2,3))
+
+def bad_level0_field_input_6():
+    """cuboid"""
+    mag = np.array([(100,200,300)])
+    dim3 = np.array([(2,2,2)])
+    obs = np.array([(1,2,3)])
+    magpy.core.magnet_cuboid_field(mag, dim3, obs, field='F')
+
 class TestExceptions(unittest.TestCase):
     """ test class for exception testing """
+    def test_level0_field_inputs(self):
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_0)
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_1)
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_2)
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_3)
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_4)
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_5)
+        self.assertRaises(MagpylibBadUserInput, bad_level0_field_input_6)
 
     def test_init(self):
         """ missing inputs when calling display and getB"""
diff --git a/tests/test_missing_plotly.py b/tests/test_missing_plotly.py
index eaf8e7819..80ef5480d 100644
--- a/tests/test_missing_plotly.py
+++ b/tests/test_missing_plotly.py
@@ -5,12 +5,17 @@
 
 
 def test_with_plotly_present():
+    """ no docstring
+    """
+    # pylint: disable=assignment-from-no-return
     src = magpy.magnet.Cylinder((1, 2, 3), (1, 2))
     x = src.show(renderer="json", backend="plotly")
     assert x is None, "display with plotly backend failed"
 
 
 def test_with_plotly_is_missing():
+    """ no docstring
+    """
     with mock.patch.dict(sys.modules, {"plotly": None}):
         src = magpy.magnet.Cylinder((1, 2, 3), (1, 2))
         with pytest.raises(ModuleNotFoundError):
diff --git a/tests/test_obj_Collection.py b/tests/test_obj_Collection.py
index 37c71d64d..42bf18b63 100644
--- a/tests/test_obj_Collection.py
+++ b/tests/test_obj_Collection.py
@@ -6,37 +6,38 @@
 import magpylib as magpy
 from magpylib._src.exceptions import MagpylibBadUserInput
 
-# # GENERATE TESTDATA
-# N = 5
-# mags = (np.random.rand(N,6,3)-0.5)*1000
-# dims3 = (np.random.rand(N,3,3)-0.5)*5     # 5x cuboid
-# dims2 = (np.random.rand(N,3,2)-0.5)*5     # 5x cylinder
-# posos = (np.random.rand(N,23,3)-0.5)*10 #readout at 333 positions
-
-# angs =  (np.random.rand(N,18)-0.5)*2*10 # each step rote by max 10 deg
-# axs =   (np.random.rand(N,18,3)-0.5)
-# anchs = (np.random.rand(N,18,3)-0.5)*5.5
-# movs =  (np.random.rand(N,18,3)-0.5)*0.5
-# rvs = (np.random.rand(N,3)-.5)*0.1
-
-# B = []
-# for mag,dim2,dim3,ang,ax,anch,mov,poso,rv in zip(mags,dims2,dims3,angs,axs,anchs,movs,posos,rvs):
-#     rot = R.from_rotvec(rv)
-#     pm1b = magpy.magnet.Cuboid(mag[0],dim3[0])
-#     pm2b = magpy.magnet.Cuboid(mag[1],dim3[1])
-#     pm3b = magpy.magnet.Cuboid(mag[2],dim3[2])
-#     pm4b = magpy.magnet.Cylinder(mag[3],dim2[0])
-#     pm5b = magpy.magnet.Cylinder(mag[4],dim2[1])
-#     pm6b = magpy.magnet.Cylinder(mag[5],dim2[2])
-
-#     # 18 subsequent operations
-#     for a,aa,aaa,mv in zip(ang,ax,anch,mov):
-#         for pm in [pm1b,pm2b,pm3b,pm4b,pm5b,pm6b]:
-#             pm.move(mv).rotate_from_angax(a,aa,aaa).rotate(rot,aaa)
-#     B += [magpy.getB([pm1b,pm2b,pm3b,pm4b,pm5b,pm6b], poso, sumup=True, niter=100)]
-# B = np.array(B)
-# inp = [mags,dims2,dims3,posos,angs,axs,anchs,movs,rvs,B]
-# pickle.dump(inp,open('testdata_Collection.p', 'wb'))
+# # # GENERATE TESTDATA
+# # N = 5
+# # mags = (np.random.rand(N,6,3)-0.5)*1000
+# # dims3 = np.random.rand(N,3,3)*5     # 5x cuboid
+# # dims2 = np.random.rand(N,3,2)*5     # 5x cylinder
+# # posos = (np.random.rand(N,23,3)-0.5)*10 #readout at 333 positions
+
+# # angs =  (np.random.rand(N,18)-0.5)*2*10 # each step rote by max 10 deg
+# # axs =   (np.random.rand(N,18,3)-0.5)
+# # anchs = (np.random.rand(N,18,3)-0.5)*5.5
+# # movs =  (np.random.rand(N,18,3)-0.5)*0.5
+# # rvs = (np.random.rand(N,3)-.5)*0.1
+
+# # B = []
+# # for mag,dim2,dim3,ang,ax,anch,mov,poso,rv in zip(
+# #        mags,dims2,dims3,angs,axs,anchs,movs,posos,rvs):
+# #     rot = R.from_rotvec(rv)
+# #     pm1b = magpy.magnet.Cuboid(mag[0],dim3[0])
+# #     pm2b = magpy.magnet.Cuboid(mag[1],dim3[1])
+# #     pm3b = magpy.magnet.Cuboid(mag[2],dim3[2])
+# #     pm4b = magpy.magnet.Cylinder(mag[3],dim2[0])
+# #     pm5b = magpy.magnet.Cylinder(mag[4],dim2[1])
+# #     pm6b = magpy.magnet.Cylinder(mag[5],dim2[2])
+
+# #     # 18 subsequent operations
+# #     for a,aa,aaa,mv in zip(ang,ax,anch,mov):
+# #         for pm in [pm1b,pm2b,pm3b,pm4b,pm5b,pm6b]:
+# #             pm.move(mv).rotate_from_angax(a,aa,aaa).rotate(rot,aaa)
+# #     B += [magpy.getB([pm1b,pm2b,pm3b,pm4b,pm5b,pm6b], poso, sumup=True)]
+# # B = np.array(B)
+# # inp = [mags,dims2,dims3,posos,angs,axs,anchs,movs,rvs,B]
+# # pickle.dump(inp,open('testdata_Collection.p', 'wb'))
 
 
 def test_Collection_basics():
diff --git a/tests/test_obj_Cuboid.py b/tests/test_obj_Cuboid.py
index 1638e51df..15d204198 100644
--- a/tests/test_obj_Cuboid.py
+++ b/tests/test_obj_Cuboid.py
@@ -5,45 +5,45 @@
 from magpylib.magnet import Cuboid
 from magpylib._src.obj_classes.class_Sensor import Sensor
 
-# """data generation for test_Cuboid()"""
 
-# N = 100
+# # # """data generation for test_Cuboid()"""
 
-# mags = (np.random.rand(N,3)-0.5)*1000
-# dims = (np.random.rand(N,3)-0.5)*5
-# posos = (np.random.rand(N,333,3)-0.5)*10 #readout at 333 positions
+# # N = 100
 
-# angs =  (np.random.rand(N,18)-0.5)*2*10 # each step rote by max 10 deg
-# axs =   (np.random.rand(N,18,3)-0.5)
-# anchs = (np.random.rand(N,18,3)-0.5)*5.5
-# movs =  (np.random.rand(N,18,3)-0.5)*0.5
+# # mags  = (np.random.rand(N,3)-0.5)*1000
+# # dims  =  np.random.rand(N,3)*5
+# # posos = (np.random.rand(N,333,3)-0.5)*10 #readout at 333 positions
 
-# B = []
-# for mag,dim,ang,ax,anch,mov,poso in zip(mags,dims,angs,axs,anchs,movs,posos):
-#     pm = Cuboid(mag,dim)
+# # angs  =  (np.random.rand(N,18)-0.5)*2*10 # each step rote by max 10 deg
+# # axs   =   (np.random.rand(N,18,3)-0.5)
+# # anchs = (np.random.rand(N,18,3)-0.5)*5.5
+# # movs  =  (np.random.rand(N,18,3)-0.5)*0.5
 
-#     # 18 subsequent operations
-#     for a,aa,aaa,mv in zip(ang,ax,anch,mov):
-#         pm.move(mv).rotate_from_angax(a,aa,aaa)
+# # B = []
+# # for mag,dim,ang,ax,anch,mov,poso in zip(mags,dims,angs,axs,anchs,movs,posos):
+# #     pm = magpy.magnet.Cuboid(mag,dim)
 
-#     B += [pm.getB(poso)]
-# B = np.array(B)
+# #     # 18 subsequent operations
+# #     for a,aa,aaa,mv in zip(ang,ax,anch,mov):
+# #         pm.move(mv).rotate_from_angax(a,aa,aaa)
 
-# inp = [mags,dims,posos,angs,axs,anchs,movs,B]
+# #     B += [pm.getB(poso)]
+# # B = np.array(B)
 
-# pickle.dump(inp,open(os.path.abspath('./../'testdata_Cuboid.p', 'wb'))
+# # inp = [mags,dims,posos,angs,axs,anchs,movs,B]
 
+# # pickle.dump(inp, open(os.path.abspath('testdata_Cuboid.p'), 'wb'))
 
 def test_Cuboid_basics():
-    """ test Cuboid fundamentals, test against magpylib2 fields
+    """ test Cuboid fundamentals
     """
-    # data generated below
+    # data generated in comment above
     data = pickle.load(open(os.path.abspath('./tests/testdata/testdata_Cuboid.p'), 'rb'))
     mags,dims,posos,angs,axs,anchs,movs,B = data
 
     btest = []
     for mag,dim,ang,ax,anch,mov,poso in zip(mags,dims,angs,axs,anchs,movs,posos):
-        pm = Cuboid(mag,dim)
+        pm = magpy.magnet.Cuboid(mag,np.abs(dim))
 
         # 18 subsequent operations
         for a,aa,aaa,mv in zip(ang,ax,anch,mov):
@@ -52,7 +52,7 @@ def test_Cuboid_basics():
         btest += [pm.getB(poso)]
     btest = np.array(btest)
 
-    assert np.allclose(B, btest), "test_Cuboid failed big time"
+    np.testing.assert_allclose(B, btest)
 
 
 def test_Cuboid_add():
@@ -87,3 +87,23 @@ def test_repr_cuboid():
     pm1.style.name = 'cuboid_01'
     assert pm1.__repr__()[:6] == 'Cuboid', 'Cuboid repr failed'
     assert "name='cuboid_01'" in pm1.__repr__(), 'Cuboid repr failed'
+
+
+def test_cuboid_object_vs_lib():
+    """
+    includes a test of the input copy problem
+    """
+
+    a = 1
+    mag = np.array([(10,20,30)])
+    dim = np.array([(  a,  a,  a)])
+    pos = np.array([(2*a,2*a,2*a)])
+    B0 = magpy.core.magnet_cuboid_field(mag, dim, pos)
+    H0 = magpy.core.magnet_cuboid_field(mag, dim, pos, field='H')
+
+    src = magpy.magnet.Cuboid(mag[0], dim[0])
+    B1 = src.getB(pos)
+    H1 = src.getH(pos)
+
+    np.testing.assert_allclose(B0[0], B1)
+    np.testing.assert_allclose(H0[0], H1)
diff --git a/tests/test_obj_Cylinder.py b/tests/test_obj_Cylinder.py
index 9decf7646..84f99e395 100644
--- a/tests/test_obj_Cylinder.py
+++ b/tests/test_obj_Cylinder.py
@@ -1,7 +1,6 @@
 import numpy as np
 import magpylib as magpy
 from magpylib.magnet import Cylinder
-from magpylib._src.fields.field_wrap_BH_level2_dict import getB_dict, getH_dict
 
 
 def test_Cylinder_add():
@@ -79,27 +78,27 @@ def test_Cylinder_getBH():
         B1 = src2.getB(poso)
         H1 = src2.getH(poso)
 
-        B2 = getB_dict(
-            source_type='Cylinder',
+        B2 = magpy.getB(
+            'Cylinder',
+            poso,
             magnetization=mag,
-            dimension=d2,
-            observer=poso)
-        H2 = getH_dict(
-            source_type='Cylinder',
+            dimension=d2,)
+        H2 = magpy.getH(
+            'Cylinder',
+            poso,
             magnetization=mag,
-            dimension=d2,
-            observer=poso)
+            dimension=d2,)
 
-        B3 = getB_dict(
-            source_type='CylinderSegment',
+        B3 = magpy.getB(
+            'CylinderSegment',
+            poso,
             magnetization=mag,
-            dimension=d5,
-            observer=poso)
-        H3 = getH_dict(
-            source_type='CylinderSegment',
+            dimension=d5,)
+        H3 = magpy.getH(
+            'CylinderSegment',
+            poso,
             magnetization=mag,
-            dimension=d5,
-            observer=poso)
+            dimension=d5,)
 
         assert np.allclose(B1, B2)
         assert np.allclose(B1, B3)
diff --git a/tests/test_obj_Loop.py b/tests/test_obj_Loop.py
index 023b73620..7fdfb4811 100644
--- a/tests/test_obj_Loop.py
+++ b/tests/test_obj_Loop.py
@@ -47,16 +47,16 @@ def test_Loop_basic_H():
     assert np.allclose(H, Htest)
 
 
-def test_Cicular_problem_positions():
-    """ Loop on z and on loop
-    """
-    src = magpy.current.Loop(current=1, diameter=2)
-    sens = magpy.Sensor()
-    sens.move([[0,1,0],[1,0,0]], start=1)
-
-    B = src.getB(sens)
-    Btest = np.array([[0,0,0.6283185307179586], [0,0,0], [0,0,0]])
-    assert np.allclose(B, Btest)
+# def test_Cicular_problem_positions():
+#     """ Loop on z and on loop
+#     """
+#     src = magpy.current.Loop(current=1, diameter=2)
+#     sens = magpy.Sensor()
+#     sens.move([[0,1,0],[1,0,0]], start=1)
+
+#     B = src.getB(sens)
+#     Btest = np.array([[0,0,0.6283185307179586], [0,0,0], [0,0,0]])
+#     assert np.allclose(B, Btest)
 
 
 def test_repr():
diff --git a/tests/test_obj_Sphere.py b/tests/test_obj_Sphere.py
index e87b911d1..82dd6f5fe 100644
--- a/tests/test_obj_Sphere.py
+++ b/tests/test_obj_Sphere.py
@@ -5,33 +5,33 @@
 from magpylib import Collection
 import magpylib as magpy
 
-# """data generation for test_Sphere()"""
+# # """data generation for test_Sphere()"""
 
-# N = 100
+# # N = 100
 
-# mags = (np.random.rand(N,3)-0.5)*1000
-# dims = (np.random.rand(N)-0.5)*5
-# posos = (np.random.rand(N,333,3)-0.5)*10 #readout at 333 positions
+# # mags = (np.random.rand(N,3)-0.5)*1000
+# # dims = np.random.rand(N)*5
+# # posos = (np.random.rand(N,333,3)-0.5)*10 #readout at 333 positions
 
-# angs =  (np.random.rand(N,18)-0.5)*2*10 # each step rote by max 10 deg
-# axs =   (np.random.rand(N,18,3)-0.5)
-# anchs = (np.random.rand(N,18,3)-0.5)*5.5
-# movs =  (np.random.rand(N,18,3)-0.5)*0.5
+# # angs =  (np.random.rand(N,18)-0.5)*2*10 # each step rote by max 10 deg
+# # axs =   (np.random.rand(N,18,3)-0.5)
+# # anchs = (np.random.rand(N,18,3)-0.5)*5.5
+# # movs =  (np.random.rand(N,18,3)-0.5)*0.5
 
-# B = []
-# for mag,dim,ang,ax,anch,mov,poso in zip(mags,dims,angs,axs,anchs,movs,posos):
-#     pm = Sphere(mag,dim)
+# # B = []
+# # for mag,dim,ang,ax,anch,mov,poso in zip(mags,dims,angs,axs,anchs,movs,posos):
+# #     pm = magpy.magnet.Sphere(mag,dim)
 
-#     # 18 subsequent operations
-#     for a,aa,aaa,mv in zip(ang,ax,anch,mov):
-#         pm.move(mv).rotate_from_angax(a,aa,aaa)
+# #     # 18 subsequent operations
+# #     for a,aa,aaa,mv in zip(ang,ax,anch,mov):
+# #         pm.move(mv).rotate_from_angax(a,aa,aaa)
 
-#     B += [pm.getB(poso)]
-# B = np.array(B)
+# #     B += [pm.getB(poso)]
+# # B = np.array(B)
 
-# inp = [mags,dims,posos,angs,axs,anchs,movs,B]
+# # inp = [mags,dims,posos,angs,axs,anchs,movs,B]
 
-# pickle.dump(inp,open('testdata_Sphere.p', 'wb'))
+# # pickle.dump(inp, open('testdata_Sphere.p', 'wb'))
 
 def test_Sphere_basics():
     """ test Cuboid fundamentals, test against magpylib2 fields
@@ -84,3 +84,19 @@ def test_repr():
     """
     pm3 = magpy.magnet.Sphere((1,2,3),3)
     assert pm3.__repr__()[:6] == 'Sphere', 'Sphere repr failed'
+
+
+def test_sphere_object_vs_lib():
+    """
+    tests object vs lib computation
+    this also checks if np.int (from array slice) is allowed as input
+    """
+    mag = np.array([(10,20,30)])
+    dia = np.array([1])
+    pos = np.array([(2,2,2)])
+    B1 = magpy.core.magnet_sphere_field(mag, dia, pos)[0]
+
+    src = magpy.magnet.Sphere(mag[0], dia[0])
+    B2 = src.getB(pos)
+
+    np.testing.assert_allclose(B1, B2)
diff --git a/tests/test_physics_consistency.py b/tests/test_physics_consistency.py
index 8ff6210b2..432f3c0bd 100644
--- a/tests/test_physics_consistency.py
+++ b/tests/test_physics_consistency.py
@@ -1,6 +1,6 @@
 import numpy as np
 import magpylib as magpy
-from magpylib._src.fields.field_wrap_BH_level2_dict import getB_dict
+
 
 def test_dipole_approximation():
     """ test if all source fields converge towards the correct dipole field at distance
@@ -191,7 +191,7 @@ def test_Line_vs_Loop():
     # field from line currents
     Bls = []
     for p in po:
-        Bl = getB_dict(source_type='Line', observer=p, current=1,
+        Bl = magpy.getB('Line', p, current=1,
             segment_start=ps, segment_end=pe)
         Bls += [np.sum(Bl, axis=0)]
     Bls = np.array(Bls)
@@ -222,7 +222,7 @@ def Binf(i0, pos):
     pe = (0,0,1000000)
     Bls, Binfs = [], []
     for p in pos_obs:
-        Bls += [getB_dict(source_type='Line', observer=p, current=1,
+        Bls += [magpy.getB('Line', p, current=1,
             segment_start=ps, segment_end=pe)]
         Binfs += [Binf(1,p)]
     Bls = np.array(Bls)
diff --git a/tests/testdata/testdata_Collection.p b/tests/testdata/testdata_Collection.p
index 96aaa670e..3392d8910 100644
Binary files a/tests/testdata/testdata_Collection.p and b/tests/testdata/testdata_Collection.p differ
diff --git a/tests/testdata/testdata_Cuboid.p b/tests/testdata/testdata_Cuboid.p
index 5ab9da5c5..e0be44efc 100644
Binary files a/tests/testdata/testdata_Cuboid.p and b/tests/testdata/testdata_Cuboid.p differ
diff --git a/tests/testdata/testdata_Sphere.p b/tests/testdata/testdata_Sphere.p
index 0417604f1..73cd7212c 100644
Binary files a/tests/testdata/testdata_Sphere.p and b/tests/testdata/testdata_Sphere.p differ