Repository: magpylib/magpylib (language: python)

Issue:
Add `getM()` and `getJ()` - separate inside/outside checks from core functions
  Core functions are trimmed to performance, the mask generation for inside-outside-edge checks take a lot of time. In most cases users know whether they are inside or outside. One more use of it is to refactor the code to enable the calculation of `getM` and `getJ`.

- [x] can inside-outside checks be separated from core computation ?
- [x] can they be turned on and off in a practical manner ?
- [x] Add `getM()`
- [x] Add `getJ()`

> [!NOTE]  
> `in_out` is currently an unused option in `magnet_trimesh_field`

### to clarify:

- #703

- `Circle`, `Polyline`, `Triangle`, `Dipole` -> what are `getM` and `getJ` outputs ? 
- `Sphere` -> cannot use a generic `convert_BHMJ` because of `polarization/3` factor. What about this factor when calculating  `getM` and `getJ`
Move to SI Units  + Redefine magnetization and add codependency with polarization
# Related Issues
- #698 
- #694 

# Notes
This PR includes two changes:
- Changing all inputs and outputs to SI units!
- Adding the `polarization` attribute to magnets which is co-dependent with the `magnetization` attribute

**BREAKING CHANGES** - this will be part of Magpylib v5
Display option for multi-pixel sensors
![image](https://user-images.githubusercontent.com/47384012/112679222-33b45300-8e6c-11eb-8da0-9c23676d7a30.png)



Discussion / hints:
- moving to #647
Reopening this as not addressed. Can still be tracked in  #647



Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9c43c6237..8f4e665e4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,12 +4,15 @@ All notable changes to magpylib are documented here.
 # Changelog
 
 ## [5.0.0dev]
+- The Magpylib inputs and outputs are now in **SI Units**. The `magnetization` term has also been redefined and is now codependent with the `polarization` term ([#712](https://github.com/magpylib/magpylib/issues/712))
+- Added `getM` (magnetization) and `getJ` (polarization) top level functions. Also an `in_out` (inside/ouside) parameter is added to specify the location of the observers relative to the magnet body in order to increase performance ([#717](https://github.com/magpylib/magpylib/issues/717))
+-
 
 ## [4.5.1] - 2023-12-28
-- Fixed a field computatio issue where H-field resulting from axial magnetization is computed incorrectly inside of Cylinders ([#703](https://github.com/magpylib/magpylib/issues/703))
+- Fixed a field computation issue where H-field resulting from axial magnetization is computed incorrectly inside of Cylinders ([#703](https://github.com/magpylib/magpylib/issues/703))
 
 ## [4.5.0] - 2023-12-13
-- Add optional handedness parameter for Sensors ([#687](https://github.com/magpylib/magpylib/pull/687))
+- Added optional handedness parameter for Sensors ([#687](https://github.com/magpylib/magpylib/pull/687))
 - Renaming classes: `Line`→`Polyline`, `Loop`→`Circle`. Old names are still valid but will issue a `DeprecationWarning` and will eventually be removed in the next major version ([#690](https://github.com/magpylib/magpylib/pull/690))
 - Rework CI/CD workflows ([#686](https://github.com/magpylib/magpylib/pull/686))
 
@@ -29,7 +32,7 @@ All notable changes to magpylib are documented here.
 - Many minor graphic improvements ([#663](https://github.com/magpylib/magpylib/pull/663), [#649](https://github.com/magpylib/magpylib/issues/649), [#653](https://github.com/magpylib/magpylib/issues/653))
 - `legend` style option ([#650](https://github.com/magpylib/magpylib/issues/650))
 - Changed unit naming in text to comply with DIN Norm 641 ([#614](https://github.com/magpylib/magpylib/issues/614))
-- Improving the documentation now boasting a contribution guide, a news-blog, an example and tutorial gallery, a getting started section and many other improvements ([#621](https://github.com/magpylib/magpylib/issues/621), [#596](https://github.com/magpylib/magpylib/issues/596), [#580](https://github.com/magpylib/magpylib/issues/580))
+- Improved documentation now boasting a contribution guide, a news-blog, an example and tutorial gallery, a getting started section and many other improvements ([#621](https://github.com/magpylib/magpylib/issues/621), [#596](https://github.com/magpylib/magpylib/issues/596), [#580](https://github.com/magpylib/magpylib/issues/580))
 - Improved numerical stability of `CylinderSegement`, ([#648](https://github.com/magpylib/magpylib/issues/648), [#651](https://github.com/magpylib/magpylib/issues/651))
 
 
diff --git a/docs/_pages/docu/docu_magpylib_api.md b/docs/_pages/docu/docu_magpylib_api.md
index a7dd66ba1..753a0fdf6 100644
--- a/docs/_pages/docu/docu_magpylib_api.md
+++ b/docs/_pages/docu/docu_magpylib_api.md
@@ -15,6 +15,8 @@ Magpylib requires no special input format. All scalar types (`int`, `float`, ...
 (docu-units)=
 ## Units
 
+The important quantity <span style="color: orange">**mu_0**</span> lies at the top level. It's value is not $4 \pi 10^{-7}$ since [the redefinition of the SI base units](https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units), but a value close to it.
+
 For historical reasons Magpylib used non-SI units until Version 4. Starting with version 5 all inputs and outputs are SI-based.
 
 ::::{grid} 3
@@ -26,30 +28,30 @@ For historical reasons Magpylib used non-SI units until Version 4. Starting with
 :columns: 10
 | PHYSICAL QUANTITY | MAGPYLIB PARAMETER | UNITS from v5| UNITS until v4|
 |:---:|:---:|:---:|:---:|
-| Magnetic Polarization $\vec{J}$  | `polarization`                | **T**      | -        |
-| Magnetization $\vec{M}$          | `magnetization`               | **A/m**    | mT       |
+| Magnetic Polarization $\vec{J}$  | `polarization`, `getJ()`      | **T**      | -        |
+| Magnetization $\vec{M}$          | `magnetization`, `getM()`     | **A/m**    | mT       |
 | Electric Current $i_0$           | `current`                     | **A**      | A        |
 | Magnetic Dipole Moment $\vec{m}$ | `moment`                      | **A·m²**   | mT·mm³   |
 | B-field $\vec{B}$                | `getB()`                      | **T**      | mT       |
 | H-field $\vec{H}$                | `getH()`                      | **A/m**    | kA/m     |
-| Length-inputs                    | `position`, `dimension`, ...  | **m**      | mm       |
+| Length-inputs                    | `position`, `dimension`, `vertices`, ...  | **m**      | mm       |
 | Angle-inputs                     | `angle`, `dimension`, ...     | **°**      | °        |
 :::
 
 ::::
 
 ```{warning}
-Up to version 4, Magpylib was unfortunately contributing to the naming confusion in magnetism that is explained well [here](https://www.e-magnetica.pl/doku.php/confusion_between_b_and_h). The input `magnetization` in Magpylib was refering to the magnetic polarization (and not the magnetization), the difference being only in the physical unit. From version 5 onwards this was fixed.
+Up to version 4, Magpylib was unfortunately contributing to the naming confusion in magnetism that is explained well [here](https://www.e-magnetica.pl/doku.php/confusion_between_b_and_h). The input `magnetization` in Magpylib < v5 was refering to the magnetic polarization (and not the magnetization), the difference being only in the physical unit. From version 5 onwards this is fixed.
 ```
 
 ```{note}
-All input and output units in Magpylib (version 5 and higher) are SI-based, see table above. However, for advanced use one should be aware that the analytical solutions are scale invariant - _"a magnet with 1 mm sides creates the same field at 1 mm distance as a magnet with 1 m sides at 1 m distance"_. The choice of length input unit is therefore not relevant, but it is critical to abide by the same length unit across one computation.
+All input and output units in Magpylib (version 5 and higher) are SI-based, see table above. However, for advanced use one should be aware that the analytical solutions are scale invariant - _"a magnet with 1 mm sides creates the same field at 1 mm distance as a magnet with 1 m sides at 1 m distance"_. The choice of length input unit is therefore not relevant, but it is critical to keep the same length unit for all inputs in one computation.
 
-In addition, `getB` returns the same unit as given by the `polarization` input. With polarization input in mT, getB will return mT as well. At the same time when the `magnetization` input is kA/m, then `getH` returns kA/m as well. The B/H-field outputs are related to a M/J-inputs via a factor of $µ_0=4\pi\cdot 10^{-7}$.
+In addition, `getB` returns the same unit as given by the `polarization` input. With polarization input in mT, getB will return mT as well. At the same time when the `magnetization` input is kA/m, then `getH` returns kA/m as well. The B/H-field outputs are related to a M/J-inputs via a factor of $µ_0$.
 ```
 
 ```{note}
-The connection between the magnetic polarization J, the magnetization M and the material parameters of a permanent magnet are shown in {ref}`gallery-tutorial-magnetmodel`.
+The connection between the magnetic polarization J, the magnetization M and the material parameters of a real permanent magnet are shown in {ref}`gallery-tutorial-magnetmodel`.
 ```
 
 <!-- ################################################################## -->
@@ -88,7 +90,7 @@ See {ref}`docu-position` for more information on these features.
 
 See {ref}`docu-graphics` for more information on graphic output, default styles and customization possibilities.
 
-* The <span style="color: orange">**getB()**</span> and <span style="color: orange">**getH()**</span> methods give quick access to field computation.
+* The <span style="color: orange">**getB()**</span>, <span style="color: orange">**getH()**</span>, <span style="color: orange">**getJ()**</span> and <span style="color: orange">**getM()**</span> methods give quick access to field computation.
 
 See {ref}`docu-field-computation` for more information.
 
@@ -603,7 +605,7 @@ The tutorial {ref}`gallery-tutorial-paths` shows intuitive good practice example
 
 <hr style="border:3px solid gray">
 
-Magnetic field computation in Magpylib is done via two top-level functions <span style="color: orange">**getB**</span> and <span style="color: orange">**getH**</span>.
+Magnetic field computation in Magpylib is done via four top-level functions <span style="color: orange">**getB**</span>, <span style="color: orange">**getH**</span>, <span style="color: orange">**getJ**</span> and, <span style="color: orange">**getM**</span>.
 
 
 ```python
@@ -614,7 +616,15 @@ magpylib.getB(sources, observers, squeeze=True, pixel_agg=None, output="ndarray"
 magpylib.getH(sources, observers, squeeze=True, pixel_agg=None, output="ndarray")
 ```
 
-that compute the magnetic field generated by `sources` as seen by the `observers` in their local coordinates. `sources` can be any Magpylib source object (e.g. magnets) or a flat list thereof. `observers` can be an array of position vectors with shape `(n1,n2,n3,...,3)`, any Magpylib observer object (e.g. sensors), or a flat list thereof. The following code shows a minimal example for Magplyib field computation.
+```python
+magpylib.getJ(sources, observers, squeeze=True, pixel_agg=None, output="ndarray")
+```
+
+```python
+magpylib.getM(sources, observers, squeeze=True, pixel_agg=None, output="ndarray")
+```
+
+that compute the respective field generated by `sources` as seen by the `observers` in their local coordinates. `sources` can be any Magpylib source object (e.g. magnets) or a flat list thereof. `observers` can be an array of position vectors with shape `(n1,n2,n3,...,3)`, any Magpylib observer object (e.g. sensors), or a flat list thereof. The following code shows a minimal example for Magplyib field computation.
 
 ::::{grid}
 :gutter: 5
@@ -642,10 +652,10 @@ print(B)
 :::
 ::::
 
-By default, `getB` returns the B-field in units of T, `getH` the H-field in units of A/m.
+By default, `getB` returns the B-field in units of T, `getH` the H-field in units of A/m, `getJ` the magnetic polarization in T and, `getM` the magnetization in A/m, assuming that all inputs are given in SI units. 
 
 ```{note}
-In reality, `getB` returns the same unit as given by the `polarization` input. For example, with polarization input in mT, getB will return mT as well. At the same time when the `magnetization` input is kA/m, then `getH` returns kA/m as well. The B/H-field outputs are related to a M/J-inputs via a factor of $µ_0=4\pi\cdot 10^{-7}$.
+In reality, `getB` returns the same unit as given by the `polarization` input. For example, with polarization input in mT, getB will return mT as well. At the same time when the `magnetization` input is kA/m, then `getH` returns kA/m as well. The B/H-field outputs are related to a M/J-inputs via a factor of $µ_0$.
 ```
 
 The output of a field computation `magpy.getB(sources, observers)` is by default a Numpy array of shape `(l, m, k, n1, n2, n3, ..., 3)` where `l` is the number of input sources, `m` the (maximal) object path length, `k` the number of observers, `n1,n2,n3,...` the sensor pixel shape or the shape of the observer position array input and `3` the three magnetic field components $(B_x, B_y, B_z)$.
@@ -653,14 +663,14 @@ The output of a field computation `magpy.getB(sources, observers)` is by default
 
 * `squeeze`: If True (default) all axes of length 1 in the output (e.g. only a single source) are squeezed.
 
-* `pixel_agg`: Select a compatible numpy aggregator function (e.g. `"min"`, `"mean"`) that is applied to the output. For example, with `pixel_agg="mean"` the mean field of all observer points is returned. Only with this option it is possible to supply `getB` and `getH` with multiple observers that have different pixel shapes.
+* `pixel_agg`: Select a compatible numpy aggregator function (e.g. `"min"`, `"mean"`) that is applied to the output. For example, with `pixel_agg="mean"` the mean field of all observer points is returned. With this option it is possible to supply `getBHJM` with multiple observers that have different pixel shapes.
 
 * `output`: Change the output format. Options are `"ndarray"` (default, returns a Numpy ndarray) and `"dataframe"` (returns a 2D-table Pandas DataFrame).
 
 ```{note}
 Magpylib collects all inputs (object parameters), and vectorizes them for the computation which reduces the computation time dramatically for large inputs.
 
-Try to make all field computations with as few calls to `getB` and `getH` as possible. Avoid Python loops at all costs!
+Try to make all field computations with as few calls to `getBHJM` as possible. Avoid Python loops at all costs!
 ```
 
 The tutorial {ref}`gallery-tutorial-field-computation` shows good practices with Magpylib field computation.
@@ -669,7 +679,7 @@ The tutorial {ref}`gallery-tutorial-field-computation` shows good practices with
 (docu-functional-interface)=
 ## Direct interface
 
-Users can bypass the object oriented functionality of Magpylib and instead compute the field for n given parameter sets. This is done by providing the following inputs to the top level functions `getB` and `getH`,
+Users can bypass the object oriented functionality of Magpylib and instead compute the field for n given parameter sets. This is done by providing the following inputs to the top level functions `getB`, `getH`, `getJ` and, `getM`.
 
 1. `sources`: a string denoting the source type. Allowed values are the Magpylib source class names, see {ref}`docu-classes`.
 2. `observers`: array-like of shape (3,) or (n,3) giving the observer positions.
@@ -723,35 +733,44 @@ At the heart of Magpylib lies a set of core functions that are our implementatio
 :gutter: 1
 
 :::{grid-item}
-<span style="color: orange">**current_polyline_field(**</span> `field`, `observers`, `current`, `segment_start`, `segment_end`<span style="color: orange">**)**</span>
+<span style="color: orange">**magnet_cuboid_Bfield(**</span> `observers`, `dimensions`, `polarizations`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**current_circle_field(**</span> `field`, `observers`, `current`, `diameter`<span style="color: orange">**)**</span>
+<span style="color: orange">**magnet_cylinder_axial_Bfield(**</span> `z0`, `r`, `z`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**magnet_cuboid_field(**</span> `field`, `observers`, `polarization`, `dimension`<span style="color: orange">**)**</span>
+<span style="color: orange">**magnet_cylinder_diametral_Hfield(**</span> `z0`, `r`, `z`, `phi`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**magnet_cylinder_field(**</span> `field`, `observers`, `polarization`, `dimension`<span style="color: orange">**)**</span>
+<span style="color: orange">**magnet_cylinder_segmet_Hfield(**</span> `observers`, `dimensions`, `magnetizations`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**magnet_cylinder_segment_field(**</span> `field`, `observers`, `polarization`, `dimension`<span style="color: orange">**)**</span>
+<span style="color: orange">**magnet_sphere_Bfield(**</span>`observers`, `diameters`, `polarizations`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**magnet_sphere_field(**</span> `field`, `observers`, `polarization`, `diameter`<span style="color: orange">**)**</span>
+<span style="color: orange">**current_circle_Hfield(**</span>`r0`, `r`, `z`, `i0`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**magnet_tetrahedron_field(**</span> `field`, `observers`, `polarization`, `vertices`<span style="color: orange">**)**</span>
+<span style="color: orange">**current_polyline_Hfield(**</span>`observers`, `segments_start`, `segments_end`, `currents`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**dipole_field(**</span> `field`, `observers`, `moment`<span style="color: orange">**)**</span>
+<span style="color: orange">**dipole_Hfield(**</span>`observers`, `moments`<span style="color: orange">**)**</span>
 :::
+
 :::{grid-item}
-<span style="color: orange">**triangle_field(**</span> `field`, `observers`, `polarization`, `vertices`<span style="color: orange">**)**</span>
+<span style="color: orange">**triangle_Bfield(**</span>`observers`, `vertices`, `polarizations`<span style="color: orange">**)**</span>
 :::
+
 ::::
 
-The input `field` is either `"B"` or `"H`. All other inputs must be Numpy ndarrays of shape (n,x). Details can be found in the respective function docstrings. The following example demonstrates the core interface.
+All inputs must be Numpy ndarrays of shape (n,x). Details can be found in the respective function docstrings. The following example demonstrates the core interface.
 
 ::::{grid}
 :gutter: 5
@@ -768,17 +787,16 @@ import magpylib as magpy
 # All inputs and outputs in SI units
 
 # prepare input
-pol = np.array([(1,0,0)]*3)
-dim = np.array([(1,2)]*3)
-obs = np.array([(0,0,0)]*3)
+z0 = np.array([1,1])
+r  = np.array([1,1])
+z  = np.array([2,2])
 
 # compute field with core functions
-B = magpy.core.magnet_cylinder_field(field='B', observers=obs, polarization=pol, dimension=dim)
+B = magpy.core.magnet_cylinder_axial_Bfield(z0=z0, r=r, z=z).T
 
 print(B)
-#  --> [[0.5527864 0.        0.       ]
-#       [0.5527864 0.        0.       ]
-#       [0.5527864 0.        0.       ]]
+#  --> [[0.05561469 0.         0.06690167]
+#       [0.05561469 0.         0.06690167]]
 ```
 :::
 ::::
diff --git a/docs/_pages/gallery/gallery_tutorial_custom.md b/docs/_pages/gallery/gallery_tutorial_custom.md
index 66511aa70..33681cef1 100644
--- a/docs/_pages/gallery/gallery_tutorial_custom.md
+++ b/docs/_pages/gallery/gallery_tutorial_custom.md
@@ -60,8 +60,7 @@ def mono_field(field, observers):
     if field == "B":
         return B  # unit T
     elif field == "H":
-        mu0 = 4 * np.pi * 1e-7
-        H = B / mu0  # unit A/m
+        H = B / magpy.mu_0  # unit A/m
         return H
     else:
         raise ValueError("Field Value must be either B or H")
@@ -158,8 +157,7 @@ class Monopole(magpy.misc.CustomSource):
             if field == "B":
                 return B  # unit T
             elif field == "H":
-                mu0 = 4 * np.pi * 1e-7
-                H = B / mu0  # unit A/m
+                H = B / magpy.mu_0  # unit A/m
                 return H
             else:
                 raise ValueError("Field Value must be either B or H")
diff --git a/magpylib/__init__.py b/magpylib/__init__.py
index fa518fc42..531c90da2 100644
--- a/magpylib/__init__.py
+++ b/magpylib/__init__.py
@@ -37,6 +37,8 @@
     "misc",
     "getB",
     "getH",
+    "getM",
+    "getJ",
     "Sensor",
     "Collection",
     "show",
@@ -47,13 +49,15 @@
     "__credits__",
     "core",
     "graphics",
+    "mu_0",
 ]
 
 # create interface to outside of package
+from scipy.constants import mu_0
 from magpylib._src.defaults.defaults_utility import SUPPORTED_PLOTTING_BACKENDS
 from magpylib import magnet, current, misc, core, graphics
 from magpylib._src.defaults.defaults_classes import default_settings as defaults
-from magpylib._src.fields import getB, getH
+from magpylib._src.fields import getB, getH, getM, getJ
 from magpylib._src.obj_classes.class_Sensor import Sensor
 from magpylib._src.obj_classes.class_Collection import Collection
 from magpylib._src.display.display import show, show_context
diff --git a/magpylib/_src/defaults/defaults_utility.py b/magpylib/_src/defaults/defaults_utility.py
index 77e33dabf..f9017dd7c 100644
--- a/magpylib/_src/defaults/defaults_utility.py
+++ b/magpylib/_src/defaults/defaults_utility.py
@@ -154,7 +154,7 @@ def magic_to_dict(kwargs, separator="_") -> dict:
                 new_kwargs[keys[0]] = val
     for k, v in new_kwargs.items():
         if isinstance(v, dict):
-            new_kwargs[k] = magic_to_dict(v)
+            new_kwargs[k] = magic_to_dict(v, separator=separator)
     return new_kwargs
 
 
diff --git a/magpylib/_src/display/display.py b/magpylib/_src/display/display.py
index 1a6b37f72..235af3d37 100644
--- a/magpylib/_src/display/display.py
+++ b/magpylib/_src/display/display.py
@@ -140,7 +140,7 @@ def get_show_func(backend):
     )
 
 
-ROW_COL_SPECIFIC_NAMES = ("row", "col", "output", "sumup", "pixel_agg")
+ROW_COL_SPECIFIC_NAMES = ("row", "col", "output", "sumup", "pixel_agg", "in_out")
 
 
 def infer_backend(canvas):
@@ -305,10 +305,10 @@ def show(
     output: tuple or string, default="model3d"
         Can be a string or a tuple of strings specifying the plot output type. By default
         `output='model3d'` displays the 3D representations of the objects. If output is a tuple of
-        strings it must be a combination of 'B' or 'H' and 'x', 'y' and/or 'z'. When having multiple
-        coordinates, the field value is the combined vector length (e.g. `('Bx', 'Hxy', 'Byz')`)
-        'Bxy' is equivalent to sqrt(|Bx|^2 + |By|^2). A 2D line plot is then represented
-        accordingly if the objects contain at least one source and one sensor.
+        strings it must be a combination of 'B', 'H', 'M' or 'J' and 'x', 'y' and/or 'z'. When
+        having multiple coordinates, the field value is the combined vector length
+        (e.g. `('Bx', 'Hxy', 'Byz')`) 'Bxy' is equivalent to sqrt(|Bx|^2 + |By|^2). A 2D line plot
+        is then represented accordingly if the objects contain at least one source and one sensor.
 
     sumup: bool, default=True
         If True, sums the field values of the sources. Applies only if `output` is not `'model3d'`.
diff --git a/magpylib/_src/display/traces_generic.py b/magpylib/_src/display/traces_generic.py
index b267f4ff6..d6543cd32 100644
--- a/magpylib/_src/display/traces_generic.py
+++ b/magpylib/_src/display/traces_generic.py
@@ -253,6 +253,7 @@ def get_generic_traces_2D(
     col=None,
     sumup=True,
     pixel_agg=None,
+    in_out="auto",
     style_path_frames=None,
 ):
     """draws and animates sensor values over a path in a subplot"""
@@ -279,9 +280,9 @@ def get_generic_traces_2D(
         field_str, *coords_str = out
         if not coords_str:
             coords_str = list("xyz")
-        if field_str not in ("B", "H") and set(coords_str).difference(set("xyz")):
+        if field_str not in "BHMJ" and set(coords_str).difference(set("xyz")):
             raise ValueError(
-                "The `output` parameter must start with 'B' or 'H' "
+                "The `output` parameter must start with 'B', 'H', 'M', 'J' "
                 "and be followed by a combination of 'x', 'y', 'z' (e.g. 'Bxy' or ('Bxy', 'Hz') )"
                 f"\nreceived {out!r} instead"
             )
@@ -298,6 +299,7 @@ def get_generic_traces_2D(
         field=field_str,
         pixel_agg=pixel_agg,
         output="ndarray",
+        in_out=in_out,
     )
     BH_array = BH_array.swapaxes(1, 2)  # swap axes to have sensors first, path second
 
@@ -317,8 +319,7 @@ def get_obj_list_str(objs):
 
     def get_label_and_color(obj):
         style = obj.style
-        label = getattr(style, "label", None)
-        label = repr(obj) if not label else label
+        label = get_legend_label(obj, style=style)
         color = getattr(style, "color", None)
         return label, color
 
diff --git a/magpylib/_src/display/traces_utility.py b/magpylib/_src/display/traces_utility.py
index a77139ed8..d21f4edae 100644
--- a/magpylib/_src/display/traces_utility.py
+++ b/magpylib/_src/display/traces_utility.py
@@ -26,7 +26,8 @@ def get_legend_label(obj, style=None, suffix=True):
         desc = style.description.text
         if not desc:
             desc = getattr(obj, "_default_style_description", "")
-        suff = f" ({desc})"
+        if desc:
+            suff = f" ({desc})"
     return f"{name}{suff}"
 
 
@@ -581,6 +582,7 @@ def process_show_input_objs(objs, **kwargs):
         "output": "model3d",
         "sumup": True,
         "pixel_agg": "mean",
+        "in_out": "auto",
     }
     max_rows = max_cols = 1
     flat_objs = []
diff --git a/magpylib/_src/fields/__init__.py b/magpylib/_src/fields/__init__.py
index 91b9819b5..48f1c6e52 100644
--- a/magpylib/_src/fields/__init__.py
+++ b/magpylib/_src/fields/__init__.py
@@ -1,6 +1,6 @@
 """_src.fields"""
 
-__all__ = ["getB", "getH"]
+__all__ = ["getB", "getH", "getM", "getJ"]
 
 # create interface to outside of package
-from magpylib._src.fields.field_wrap_BH import getB, getH
+from magpylib._src.fields.field_wrap_BH import getB, getH, getM, getJ
diff --git a/magpylib/_src/fields/field_BH_circle.py b/magpylib/_src/fields/field_BH_circle.py
index e53b6b2d1..fcd3163d0 100644
--- a/magpylib/_src/fields/field_BH_circle.py
+++ b/magpylib/_src/fields/field_BH_circle.py
@@ -2,88 +2,104 @@
 Implementations of analytical expressions for the magnetic field of
 a circular current loop. Computation details in function docstrings.
 """
-import warnings
-
 import numpy as np
+from scipy.constants import mu_0 as MU0
 
-from magpylib._src.exceptions import MagpylibDeprecationWarning
 from magpylib._src.fields.special_cel import cel_iter
 from magpylib._src.input_checks import check_field_input
 from magpylib._src.utility import cart_to_cyl_coordinates
 from magpylib._src.utility import cyl_field_to_cart
-from magpylib._src.utility import MU0
 
 
 # CORE
-def current_circle_field(
-    *,
-    field: str,
-    observers: np.ndarray,
-    diameter: np.ndarray,
-    current: np.ndarray,
+def current_circle_Hfield(
+    r0: np.ndarray,
+    r: np.ndarray,
+    z: np.ndarray,
+    i0: np.ndarray,
 ) -> np.ndarray:
-    """Magnetic field of a circular (line) current loops.
-
-    The loop lies in the z=0 plane with the coordinate origin at its center.
+    """B-field of a circular line-current loops in Cylindrical Coordinates.
 
-    SI units are used for all inputs and outputs.
+    The loops lies in the z=0 plane with the coordinate origin at their centers.
+    The output is proportional to the electrical currents i0, and independent of
+    the length units chosen for observers and loop radii.
 
     Parameters
     ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
+    r0: ndarray, shape (n,)
+        Radii of loops.
 
-    observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
+    r: ndarray, shape (n,)
+        Radial positions of observers.
 
-    diameter: ndarray, shape (n,)
-        Diameter of loop in units of m.
+    z: ndarray, shape (n,)
+        Axial positions of observers.
 
-    current: ndarray, shape (n,)
-        Electrical current in units of A.
+    i0: ndarray, shape (n,)
+        Electrical currents.
 
     Returns
     -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the field of three different circular loops at three different positions.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> H = magpy.core.current_circle_field(
-    ...     field='H',
-    ...     observers=np.array([(0,0,0), (1,1,1), (2,2,2)]),
-    ...     diameter=np.array([1,2,3]),
-    ...     current=np.array([1,1,2])
-    ... )
-    >>> print(H)
-    [[0.         0.         1.        ]
-     [0.0496243  0.0496243  0.02124542]
-     [0.02833835 0.02833835 0.00654999]]
+    B-field: ndarray, shape (n,3)
+        B-field generated by Loops at observer positions.
 
     Notes
     -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
+    Implementation based on "Numerically stable and computationally
+    efficient expression for the magnetic field of a current loop.", M.Ortner et al,
+    Magnetism 2023, 3(1), 11-31.
+    """
+    n5 = len(r)
+
+    # express through ratios (make dimensionless, avoid large/small input values, stupid)
+    r = r / r0
+    z = z / r0
+
+    # field computation from paper
+    z2 = z**2
+    x0 = z2 + (r + 1) ** 2
+    k2 = 4 * r / x0
+    q2 = (z2 + (r - 1) ** 2) / x0
+
+    k = np.sqrt(k2)
+    q = np.sqrt(q2)
+    p = 1 + q
+    pf = k / np.sqrt(r) / q2 / 20 / r0 * 1e-6 * i0
+
+    # cel* part
+    cc = k2 * k2
+    ss = 2 * cc * q / p
+    Hr = pf * z / r * cel_iter(q, p, np.ones(n5), cc, ss, p, q)
+
+    # cel** part
+    cc = k2 * (k2 - (q2 + 1) / r)
+    ss = 2 * k2 * q * (k2 / p - p / r)
+    Hz = -pf * cel_iter(q, p, np.ones(n5), cc, ss, p, q)
+
+    # input is I -> output must be H-field
+    return np.row_stack((Hr, np.zeros(n5), Hz)) * 795774.7154594767  # *1e7/4/np.pi
 
-    Implementation based on "Numerically stable and computationally efficient expression for
-    the magnetic field of a current loop.", M.Ortner et al, Submitted to MDPI Magnetism, 2022
+
+def BHJM_circle(
+    field: str,
+    observers: np.ndarray,
+    diameter: np.ndarray,
+    current: np.ndarray,
+) -> np.ndarray:
+    """
+    - translate circle core to BHJM
+    - treat special cases
     """
 
-    bh = check_field_input(field, "current_circle_field()")
+    # allocate
+    BHJM = np.zeros_like(observers, dtype=float)
+
+    check_field_input(field)
+    if field in "MJ":
+        return BHJM
 
     r, phi, z = cart_to_cyl_coordinates(observers)
     r0 = np.abs(diameter / 2)
-    n = len(r0)
-
-    # allocate result
-    Br_tot, Bz_tot = np.zeros((2, n))
 
     # Special cases:
     # case1: loop radius is 0 -> return (0,0,0)
@@ -94,67 +110,30 @@ def current_circle_field(
     mask3 = r == 0
     if np.any(mask3):
         mask4 = mask3 * ~mask1  # only relevant if not also case1
-        Bz_tot[mask4] = (
-            0.6283185307179587
-            * r0[mask4] ** 2
-            / (z[mask4] ** 2 + r0[mask4] ** 2) ** (3 / 2)
+        BHJM[mask4, 2] = (
+            (r0[mask4] ** 2 / (z[mask4] ** 2 + r0[mask4] ** 2) ** (3 / 2))
+            * current[mask4]
+            * 0.5
         )
 
     # general case
     mask5 = ~np.logical_or(np.logical_or(mask1, mask2), mask3)
     if np.any(mask5):
-        r0 = r0[mask5]
-        r = r[mask5]
-        z = z[mask5]
-        n5 = len(r0)
-
-        # express through ratios (make dimensionless, avoid large/small input values, stupid)
-        r = r / r0
-        z = z / r0
-
-        # field computation from paper
-        z2 = z**2
-        x0 = z2 + (r + 1) ** 2
-        k2 = 4 * r / x0
-        q2 = (z2 + (r - 1) ** 2) / x0
-
-        k = np.sqrt(k2)
-        q = np.sqrt(q2)
-        p = 1 + q
-        pf = k / np.sqrt(r) / q2 / 20 / r0
-
-        # cel* part
-        cc = k2 * k2
-        ss = 2 * cc * q / p
-        Br_tot[mask5] = pf * z / r * cel_iter(q, p, np.ones(n5), cc, ss, p, q)
-
-        # cel** part
-        cc = k2 * (k2 - (q2 + 1) / r)
-        ss = 2 * k2 * q * (k2 / p - p / r)
-        Bz_tot[mask5] = -pf * cel_iter(q, p, np.ones(n5), cc, ss, p, q)
-
-    # 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 * 1e-6
-    ).T  # ugly but fast
-
-    # B or H field
-    if bh:
-        return B_cart
-
-    return B_cart / MU0
-
-
-def current_loop_field(*args, **kwargs):
-    """current_loop_field is deprecated, see current_circle_field"""
-
-    warnings.warn(
-        (
-            "current_loop_field is deprecated and will be removed in a future version, "
-            "use current_circle_field instead."
-        ),
-        MagpylibDeprecationWarning,
-        stacklevel=2,
+        BHJM[mask5] = current_circle_Hfield(
+            r0=r0[mask5],
+            r=r[mask5],
+            z=z[mask5],
+            i0=current[mask5],
+        ).T
+
+    BHJM[:, 0], BHJM[:, 1] = cyl_field_to_cart(phi, BHJM[:, 0])
+
+    if field == "H":
+        return BHJM
+
+    if field == "B":
+        return BHJM * MU0
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
     )
-    return current_circle_field(*args, **kwargs)
diff --git a/magpylib/_src/fields/field_BH_cuboid.py b/magpylib/_src/fields/field_BH_cuboid.py
index deec0e068..128a81dfe 100644
--- a/magpylib/_src/fields/field_BH_cuboid.py
+++ b/magpylib/_src/fields/field_BH_cuboid.py
@@ -3,70 +3,43 @@
 magnetized Cuboids. Computation details in function docstrings.
 """
 import numpy as np
+from scipy.constants import mu_0 as MU0
 
 from magpylib._src.input_checks import check_field_input
-from magpylib._src.utility import MU0
+
+# pylint: disable=too-many-statements
 
 
 # CORE
-def magnet_cuboid_field(
-    *,
-    field: str,
+def magnet_cuboid_Bfield(
     observers: np.ndarray,
-    dimension: np.ndarray,
-    polarization: np.ndarray,
-) -> np.ndarray:
-    """Magnetic field of homogeneously magnetized cuboids.
-
-    The cuboid sides are parallel to the coordinate axes. The geometric center of the
-    cuboid lies in the origin.
+    dimensions: np.ndarray,
+    polarizations: np.ndarray,
+):
+    """B-field of homogeneously magnetized cuboids in Cartesian Coordinates.
 
-    SI units are used for all inputs and outputs.
+    The cuboids sides are parallel to the coordinate axes. The geometric centers of the
+    cuboids lie in the origin. The output is proportional to the polarization magnitude
+    and independent of the length units chosen for observers and dimensions.
 
     Parameters
     ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
     observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
+        Observer positions (x,y,z) in Cartesian coordinates.
 
-    dimension: ndarray, shape (n,3)
-        Length of Cuboid sides in units of m.
+    dimensions: ndarray, shape (n,3)
+        Length of cuboids sides.
 
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T.
+    polarizations: ndarray, shape (n,3)
+        Magnetic polarization vectors.
 
     Returns
     -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the field of three different instances.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.core.magnet_cuboid_field(
-    ...     field='B',
-    ...     observers=np.array([(1,2,0), (2,3,4), (0,0,0)]),
-    ...     dimension=np.array([(2,2,2), (3,3,3), (4,4,4)]),
-    ...     polarization=np.array([(0,0,1), (1,0,0), (0,0,1)]),
-    ... )
-    >>> print(B)
-    [[ 0.          0.         -0.05227894]
-     [-0.00820941  0.00849123  0.011429  ]
-     [ 0.          0.          0.66666667]]
+    B-Field: ndarray, shape (n,3)
+        B-field generated by Cuboids at observer positions.
 
     Notes
     -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
     Field computations via magnetic surface charge density. Published
     several times with similar expressions:
 
@@ -89,181 +62,202 @@ def magnet_cuboid_field(
 
     Cichon: IEEE Sensors Journal, vol. 19, no. 7, April 1, 2019, p.2509
     """
-    # pylint: disable=too-many-statements
+    pol_x, pol_y, pol_z = polarizations.T
+    a, b, c = dimensions.T / 2
+    x, y, z = np.copy(observers).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(pol_x), 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 polarization-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-polarization
+    #    the 'missing' third sign is hidden in ff1x
+    bx_pol_x = pol_x * ff1x * qsigns[:, 0, 0]
+    by_pol_x = pol_x * ff2z * qsigns[:, 0, 1]
+    bz_pol_x = pol_x * ff2y * qsigns[:, 0, 2]
+    # contributions from y-polarization
+    bx_pol_y = pol_y * ff2z * qsigns[:, 1, 0]
+    by_pol_y = pol_y * ff1y * qsigns[:, 1, 1]
+    bz_pol_y = -pol_y * ff2x * qsigns[:, 1, 2]
+    # contributions from z-polarization
+    bx_pol_z = pol_z * ff2y * qsigns[:, 2, 0]
+    by_pol_z = -pol_z * ff2x * qsigns[:, 2, 1]
+    bz_pol_z = pol_z * ff1z * qsigns[:, 2, 2]
+
+    # summing all contributions
+    bx_tot = bx_pol_x + bx_pol_y + bx_pol_z
+    by_tot = by_pol_x + by_pol_y + by_pol_z
+    bz_tot = bz_pol_x + bz_pol_y + bz_pol_z
+
+    # 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
+
+    B /= 4 * np.pi
+    return B
+
+
+def BHJM_magnet_cuboid(
+    field: str,
+    observers: np.ndarray,
+    dimension: np.ndarray,
+    polarization: np.ndarray,
+) -> np.ndarray:
+    """
+    - translate cuboid core field to BHJM
+    - treat special cases
+    """
 
-    bh = check_field_input(field, "magnet_cuboid_field()")
+    RTOL_SURFACE = 1e-15  # relative distance tolerance to be considered on surface
 
-    magx, magy, magz = polarization.T
+    check_field_input(field)
+
+    pol_x, pol_y, pol_z = polarization.T
     a, b, c = np.abs(dimension.T) / 2
     x, y, z = observers.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))
+    # allocate for output
+    BHJM = polarization.astype(float)
 
-    # SPECIAL CASE 1: mag = (0,0,0)
-    mask1 = (magx == 0) * (magy == 0) * (magz == 0)  # 2x faster than np.all()
+    # SPECIAL CASE 1: polarization = (0,0,0)
+    mask_pol_not_null = ~(
+        (pol_x == 0) * (pol_y == 0) * (pol_z == 0)
+    )  # 2x faster than np.all()
 
     # SPECIAL CASE 2: 0 in dimension
-    mask2 = (a * b * c).astype(bool)
+    mask_dim_not_null = (a * b * c).astype(bool)
 
     # SPECIAL CASE 3: observer lies on-edge/corner
-    # -> 1e-15 to account for numerical imprecision 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 = polarization[mask_gen].T
-        a, b, c = dimension[mask_gen].T / 2
-        x, y, z = np.copy(observers[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-polarization
-        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-polarization
-        bx_magy = magy * ff2z * qsigns[:, 1, 0]
-        by_magy = magy * ff1y * qsigns[:, 1, 1]
-        bz_magy = -magy * ff2x * qsigns[:, 1, 2]
-        # contributions from z-polarization
-        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 polarization vector
-    mask_inside = mx2 & my2 & mz2
-    B[mask_inside] -= polarization[mask_inside]
-    H = B / MU0  # T -> A/m
-    return H
+    #   EPSILON to account for numerical imprecision when e.g. rotating
+    #   /a /b /c to account for the missing scaling (EPSILON is large when
+    #    a is e.g. EPSILON itself)
+
+    # on-surface is not a special case
+    mask_surf_x = abs(x_dist := abs(x) - a) < RTOL_SURFACE * a  # on surface
+    mask_surf_y = abs(y_dist := abs(y) - b) < RTOL_SURFACE * b  # on surface
+    mask_surf_z = abs(z_dist := abs(z) - c) < RTOL_SURFACE * c  # on surface
+
+    # inside-outside
+    mask_inside_x = x_dist < RTOL_SURFACE * a
+    mask_inside_y = y_dist < RTOL_SURFACE * b
+    mask_inside_z = z_dist < RTOL_SURFACE * c
+    mask_inside = mask_inside_x & mask_inside_y & mask_inside_z
+
+    # on edge (requires on-surface and inside-outside)
+    mask_xedge = mask_surf_y & mask_surf_z & mask_inside_x
+    mask_yedge = mask_surf_x & mask_surf_z & mask_inside_y
+    mask_zedge = mask_surf_x & mask_surf_y & mask_inside_z
+    mask_not_edge = ~(mask_xedge | mask_yedge | mask_zedge)
+
+    mask_gen = mask_pol_not_null & mask_dim_not_null & mask_not_edge
+
+    if field == "J":
+        BHJM[~mask_inside] = 0
+        return BHJM
+
+    if field == "M":
+        BHJM[~mask_inside] = 0
+        return BHJM / MU0
+
+    BHJM *= 0  # return (0,0,0) for all special cases
+    BHJM[mask_gen] = magnet_cuboid_Bfield(
+        observers=observers[mask_gen],
+        dimensions=dimension[mask_gen],
+        polarizations=polarization[mask_gen],
+    )
+    if field == "B":
+        return BHJM
+
+    if field == "H":
+        BHJM[mask_inside] -= polarization[mask_inside]
+        return BHJM / MU0
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
diff --git a/magpylib/_src/fields/field_BH_cylinder.py b/magpylib/_src/fields/field_BH_cylinder.py
index 0ab22a438..214d3f1d6 100644
--- a/magpylib/_src/fields/field_BH_cylinder.py
+++ b/magpylib/_src/fields/field_BH_cylinder.py
@@ -4,6 +4,7 @@
 """
 # pylint: disable = no-name-in-module
 import numpy as np
+from scipy.constants import mu_0 as MU0
 from scipy.special import ellipe
 from scipy.special import ellipk
 
@@ -11,28 +12,35 @@
 from magpylib._src.input_checks import check_field_input
 from magpylib._src.utility import cart_to_cyl_coordinates
 from magpylib._src.utility import cyl_field_to_cart
-from magpylib._src.utility import MU0
 
 
-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
-    CS. The geometric center of the Cylinder is in the origin.
+# CORE
+def magnet_cylinder_axial_Bfield(z0: np.ndarray, r: np.ndarray, z: np.ndarray) -> list:
+    """B-field of axially magnetized cylinders in Cylinder Coordinates.
 
-    Implementation from [Derby2009].
+    The cylinder axes coincide with the z-axis of the Cylindrical CS and the
+    geometric center of the cylinder lies in the origin. Length inputs are
+    made dimensionless by division over the cylinder radii. Unit polarization
+    is assumed.
 
     Parameters
     ----------
-    dim: ndarray, shape (n,2)
-        dimension of cylinder (d, h), diameter and height, in units of m
-    pos_obs: ndarray, shape (n,2)
-        position of observer (r,z) in cylindrical coordinates in units of m
+    z0: ndarray, shape (n,)
+        Ratios of half cylinder heights over cylinder radii.
+
+    r: ndarray, shape (n,)
+        Ratios of radial observer positions over cylinder radii.
+
+    z: Ratios of axial observer positions over cylinder radii.
 
     Returns
     -------
-    B-field: ndarray
-        B-field array of shape (n,2) in cylindrical coordinates (Br,Bz) in units of T.
+    B-field: ndarray, (Br, Bphi, Bz)
+        B-field generated by Cylinders at observer positions.
+
+    Notes
+    -----
+    Implementation based on Derby, American Journal of Physics 78.3 (2010): 229-235.
     """
     n = len(z0)
 
@@ -48,10 +56,10 @@ def fieldB_cylinder_axial(z0: np.ndarray, r: np.ndarray, z: np.ndarray) -> list:
     gamma = dmr / dpr
     one = np.ones(n)
 
-    # radial field (unit magnetization)
+    # radial field (unit polarization)
     Br = (cel(k1, one, one, -one) / sq1 - cel(k0, one, one, -one) / sq0) / np.pi
 
-    # axial field (unit magnetization)
+    # axial field (unit polarization)
     Bz = (
         1
         / dpr
@@ -62,40 +70,46 @@ def fieldB_cylinder_axial(z0: np.ndarray, r: np.ndarray, z: np.ndarray) -> list:
         / np.pi
     )
 
-    return Br, Bz
+    return np.row_stack((Br, np.zeros(n), Bz))
 
 
-def fieldH_cylinder_diametral(
+# CORE
+def magnet_cylinder_diametral_Hfield(
     z0: np.ndarray,
     r: np.ndarray,
-    phi: np.ndarray,
     z: np.ndarray,
-) -> np.ndarray:
-    """
-    H-field in Cylindrical CS of Cylinder magnet with homogenous
-    diametral unit magnetization. The Cylinder axis coincides with the z-axis of the
-    CS. The geometric center of the Cylinder is in the origin.
-
-    Implementation from [Rauber2021].
-
-    H-Field computed analytically via the magnetic scalar potential. Final integration
-    reduced to complete elliptic integrals.
+    phi: np.ndarray,
+) -> list:
+    """B-field of diametrally magnetized cylinders in Cylinder Coordinates.
 
-    Numerical Instabilities: See discussion on GitHub.
+    The cylinder axes coincide with the z-axis of the Cylindrical CS and the
+    geometric center of the cylinder lies in the origin. Length inputs are
+    made dimensionless by division over the cylinder radii. Unit magnetization
+    is assumed.
 
     Parameters
     ----------
-    dim: ndarray, shape (n,2)
-        dimension of cylinder (d, h), diameter and height, in units of m
-    tetta: ndarray, shape (n,)
-        angle between magnetization vector and x-axis in [rad]. M = (cos(tetta), sin(tetta), 0)
-    obs_pos: ndarray, shape (n,3)
-        position of observer (r,phi,z) in cylindrical coordinates in units of m and rad
+    z0: ndarray, shape (n)
+        Ratios of cylinder heights over cylinder radii.
+
+    r: ndarray, shape (n)
+        Ratios of radial observer positions over cylinder radii.
+
+    z: ndarray, shape (n)
+        Ratios of axial observer positions over cylinder radii.
+
+    phi: Azimuth angles between observers and magnetization directions.
 
     Returns
     -------
-    H-field: ndarray
-        H-field array of shape (n,3) in cylindrical coordinates (Hr, Hphi, Hz) in units of A/m.
+    H-Field: np.ndarray, (Hr, Hphi, Hz)
+        H-field generated by Cylinders at observer positions.
+
+    Notes
+    -----
+    Implementation partially based on Caciagli: Journal of Magnetism and
+    Magnetic Materials 456 (2018): 423-432, and [Ortner, Leitner, Rauber]
+    (unpublished).
     """
     # pylint: disable=too-many-statements
 
@@ -245,149 +259,112 @@ def fieldH_cylinder_diametral(
             )
         )
 
-    return Hr, Hphi, Hz
+    return np.row_stack((Hr, Hphi, Hz))
 
 
-# CORE
-def magnet_cylinder_field(
-    *,
+def BHJM_magnet_cylinder(
     field: str,
     observers: np.ndarray,
     dimension: np.ndarray,
     polarization: np.ndarray,
 ) -> np.ndarray:
-    """Magnetic field of homogeneously magnetized cylinders.
-
-    The cylinder axis coincides with the z-axis and the geometric center of the
-    cylinder lies in the origin.
-
-    SI units are used for all inputs and outputs.
-
-    Parameters
-    ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
-    observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
-
-    dimension: ndarray, shape (n,2)
-        Cylinder dimension (d,h) with diameter d and height h in units of m.
-
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T.
-
-    Returns
-    -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the B-field of two different cylinder magnets at position (1,0,0).
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.core.magnet_cylinder_field(
-    ...     field='B',
-    ...     observers=np.array([(1,0,0), (1,0,0)]),
-    ...     dimension=np.array([(1,1), (1,3)]),
-    ...     polarization=np.array([(0,0,1), (.5,0,.5)]),
-    ... )
-    >>> print(B)
-    [[ 0.          0.         -0.05185272]
-     [ 0.06821654  0.         -0.01576545]]
-
-    Notes
-    -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
-    Axial implementation based on
-
-    Derby: American Journal of Physics 78.3 (2010): 229-235.
-
-    Diametral implementation based on
-
-    Caciagli: Journal of Magnetism and Magnetic Materials 456 (2018): 423-432.
-
-    Leitner/Rauber/Orter: WIP
+    """
+    - Translate cylinder core fields to BHJM
+    - special cases
     """
 
-    bh = check_field_input(field, "magnet_cylinder_field()")
+    check_field_input(field)
 
     # transform to Cy CS --------------------------------------------
     r, phi, z = cart_to_cyl_coordinates(observers)
     r0, z0 = dimension.T / 2
 
     # scale invariance (make dimensionless)
-    r = np.copy(r / r0)
-    z = np.copy(z / r0)
-    z0 = np.copy(z0 / r0)
+    r = r / r0
+    z = z / r0
+    z0 = z0 / r0
 
-    # allocate field vectors ----------------------------------------
-    Br, Bphi, Bz = np.zeros((3, len(r)))
+    # allocate for output
+    BHJM = polarization.astype(float)
 
-    # 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
+    # inside/outside
+    mask_between_bases = np.abs(z) <= z0  # in-between top and bottom plane
+    mask_inside_hull = r <= 1  # inside Cylinder hull plane
+    mask_inside = mask_between_bases & mask_inside_hull
 
-    # special case: mag = 0
-    mask0 = np.linalg.norm(polarization, axis=1) == 0
+    if field == "J":
+        BHJM[~mask_inside] = 0
+        return BHJM
 
-    # special case: on Cylinder edge
-    mask_edge = m0 & m1
+    if field == "M":
+        BHJM[~mask_inside] = 0
+        return BHJM / MU0
 
-    # general case
-    mask_gen = ~mask0 & ~mask_edge
+    # SPECIAL CASE 1: on Cylinder edge
+    mask_on_hull = np.isclose(r, 1, rtol=1e-15, atol=0)  # on Cylinder hull plane
+    mask_on_bases = np.isclose(abs(z), z0, rtol=1e-15, atol=0)  # on top or bottom plane
+    mask_not_on_edge = ~(mask_on_hull & mask_on_bases)
 
     # axial/transv polarization cases
-    magx, magy, magz = polarization.T
-    mask_tv = (magx != 0) | (magy != 0)
-    mask_ax = magz != 0
+    pol_x, pol_y, pol_z = polarization.T
+    mask_pol_tv = (pol_x != 0) | (pol_y != 0)
+    mask_pol_ax = pol_z != 0
 
-    # inside/outside
-    mask_inside = m2 & m3
+    # SPECIAL CASE 2: pol = 0
+    mask_pol_not_null = ~((pol_x == 0) * (pol_y == 0) * (pol_z == 0))
+
+    # general case
+    mask_gen = mask_pol_not_null & mask_not_on_edge
 
     # general case masks
-    mask_tv = mask_tv & mask_gen
-    mask_ax = mask_ax & mask_gen
+    mask_pol_tv = mask_pol_tv & mask_gen
+    mask_pol_ax = mask_pol_ax & mask_gen
     mask_inside = mask_inside & mask_gen
 
-    # transversal polarization 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]
-        )
+    BHJM *= 0
 
-        # 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 polarization 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_ax * mask_inside] -= magz[mask_ax * mask_inside]
-    return np.concatenate(((Bx,), (By,), (Bz,)), axis=0).T / MU0
+    # transversal polarization contributions -----------------------
+    if any(mask_pol_tv):
+        pol_xy = np.sqrt(pol_x**2 + pol_y**2)[mask_pol_tv]
+        tetta = np.arctan2(pol_y[mask_pol_tv], pol_x[mask_pol_tv])
+
+        BHJM[mask_pol_tv] = (
+            magnet_cylinder_diametral_Hfield(
+                z0=z0[mask_pol_tv],
+                r=r[mask_pol_tv],
+                z=z[mask_pol_tv],
+                phi=phi[mask_pol_tv] - tetta,
+            )
+            * pol_xy
+        ).T
+
+    # axial polarization contributions ----------------------------
+    if any(mask_pol_ax):
+        BHJM[mask_pol_ax] += (
+            magnet_cylinder_axial_Bfield(
+                z0=z0[mask_pol_ax],
+                r=r[mask_pol_ax],
+                z=z[mask_pol_ax],
+            )
+            * pol_z[mask_pol_ax]
+        ).T
+
+    BHJM[:, 0], BHJM[:, 1] = cyl_field_to_cart(phi, BHJM[:, 0], BHJM[:, 1])
+
+    # add/subtract Mag when inside for B/H
+    if field == "B":
+        mask_tv_inside = mask_pol_tv * mask_inside
+        if any(mask_tv_inside):  # tv computes H-field
+            BHJM[mask_tv_inside, 0] += pol_x[mask_tv_inside]
+            BHJM[mask_tv_inside, 1] += pol_y[mask_tv_inside]
+        return BHJM
+
+    if field == "H":
+        mask_ax_inside = mask_pol_ax * mask_inside
+        if any(mask_ax_inside):  # ax computes B-field
+            BHJM[mask_ax_inside, 2] -= pol_z[mask_ax_inside]
+        return BHJM / MU0
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
diff --git a/magpylib/_src/fields/field_BH_cylinder_segment.py b/magpylib/_src/fields/field_BH_cylinder_segment.py
index f3f8cf9ce..c7460d770 100644
--- a/magpylib/_src/fields/field_BH_cylinder_segment.py
+++ b/magpylib/_src/fields/field_BH_cylinder_segment.py
@@ -4,13 +4,13 @@
 # pylint: disable=no-name-in-module
 # pylint: disable=too-many-statements
 import numpy as np
+from scipy.constants import mu_0 as MU0
 from scipy.special import ellipeinc
 from scipy.special import ellipkinc
 
-from magpylib._src.fields.field_BH_cylinder import magnet_cylinder_field
+from magpylib._src.fields.field_BH_cylinder import BHJM_magnet_cylinder
 from magpylib._src.fields.special_el3 import el3_angle
 from magpylib._src.input_checks import check_field_input
-from magpylib._src.utility import MU0
 
 
 def arctan_k_tan_2(k, phi):
@@ -2132,38 +2132,53 @@ def case235(r, r_i, r_bar_i, phi_bar_j, phi_bar_M, phi_bar_Mj, theta_M, z_bar_k)
     return results
 
 
-def magnet_cylinder_segment_core(
-    mag: np.ndarray, dim: np.ndarray, obs_pos: np.ndarray
+# CORE
+def magnet_cylinder_segment_Hfield(
+    observers: np.ndarray,
+    dimensions: np.ndarray,
+    magnetizations: np.ndarray,
 ) -> np.ndarray:
-    """
-    Put all solutions properly together (see paper Florian).
+    """Magnetic field of homogeneously magnetized cylinder ring segments
+    in Cartesian Coordinates.
+
+    The cylinder axes coincide with the z-axis of the Cylindrical CS and the
+    geometric center of the cylinder lies in the origin. The result is
+    proportional to the magnetization magnitude, and independent of the
+    length units used for dimensions and observers.
 
     Parameters
     ----------
-    mag: ndarray, shape (n,3)
-        magnetization vector (M, phi, th) in spherical CS, units: T, rad
-    obs_pos : ndarray, shape (n,3)
-        observer positions (r,phi,z) in cy CS, units: m, rad
-    dim: ndarray, shape (n,6)
-        segment dimension (r1,r2,phi1,phi2,z1,z2) in cy CS , units: m, rad
+    observers : ndarray, shape (n,3)
+        Observer positions (r,phi,z) in Cylinder coordinates, where phi is
+        given in rad.
+
+    dimensions: ndarray, shape (n,6)
+        Segment dimension [(r1,r2,phi1,phi2,z1,z2), ...]. r1 and r2 are
+        inner and outer radii. phi1 and phi2 are azimuth section angles
+        in rad. z1 and z2 are z-values of bottom and top.
+
+    magnetizations: ndarray, shape (n,3)
+        Magnetization vectors [(M, phi, th), ...] in spherical CS. M is the
+        magnitude of magnetization, phi and th are azimuth and polar angles
+        in rad.
 
     Returns
     -------
-    H-field: ndarray
-        H-field in cylindrical coordinates (Hr, Hphi, Hz), shape (n,3) in units of A/m.
+    H-field: ndarray, shape (n,3)
+        H-field generated by Cylinder Segments at observer positions.
 
     Notes
     -----
-    Field computed in the surface charge picture. Final integrals reduced to incomplete
-    elliptic integrals.
+    Implementation based on F.Slanovc, Journal of Magnetism and Magnetic
+    Materials, Volume 559, 1 October 2022, 169482
     """
 
     # tile inputs into 8-stacks (boundary cases)
-    r, phi, z = np.repeat(obs_pos, 8, axis=0).T
-    r_i = np.repeat(dim[:, :2], 4)
-    phi_j = np.repeat(np.tile(dim[:, 2:4], 2), 2)
-    z_k = np.ravel(np.tile(dim[:, 4:6], 4))
-    _, phi_M, theta_M = np.repeat(mag, 8, axis=0).T
+    r, phi, z = np.repeat(observers, 8, axis=0).T
+    r_i = np.repeat(dimensions[:, :2], 4)
+    phi_j = np.repeat(np.tile(dimensions[:, 2:4], 2), 2)
+    z_k = np.ravel(np.tile(dimensions[:, 4:6], 4))
+    _, phi_M, theta_M = np.repeat(magnetizations, 8, axis=0).T
 
     # initialize results array with nan
     result = np.empty((len(r), 3, 3))
@@ -2292,33 +2307,32 @@ def magnet_cylinder_segment_core(
     result = np.sum(result[:, (1, 2, 4, 7)] - result[:, (0, 3, 5, 6)], axis=(1, 3))
 
     # multiply with magnetization amplitude
-    result = result.T * mag[:, 0] / (4 * np.pi)
+    result = result.T * magnetizations[:, 0] * 1e-7
 
     return result.T
 
 
-def magnet_cylinder_segment_field_internal(
+def BHJM_cylinder_segment_internal(
     field: str,
     observers: np.ndarray,
     polarization: np.ndarray,
     dimension: np.ndarray,
 ) -> np.ndarray:
     """
-    internal version of magnet_cylinder_segment_field used for object oriented interface.
+    internal version of BHJM_cylinder_segment used for object oriented interface.
 
     Falls back to magnet_cylinder_field whenever the section angles describe the full
     360° cylinder.
     """
-    n = len(polarization)
 
-    BHfinal = np.zeros((n, 3))
+    BHfinal = np.zeros_like(observers, dtype=float)
 
     r1, r2, h, phi1, phi2 = dimension.T
 
     # case1: segment
     mask1 = (phi2 - phi1) < 360
 
-    BHfinal[mask1] = magnet_cylinder_segment_field(
+    BHfinal[mask1] = BHJM_cylinder_segment(
         field=field,
         observers=observers[mask1],
         polarization=polarization[mask1],
@@ -2327,7 +2341,7 @@ def magnet_cylinder_segment_field_internal(
 
     # case2: full cylinder
     mask1x = ~mask1
-    BHfinal[mask1x] = magnet_cylinder_field(
+    BHfinal[mask1x] = BHJM_magnet_cylinder(
         field=field,
         observers=observers[mask1x],
         polarization=polarization[mask1x],
@@ -2336,7 +2350,7 @@ def magnet_cylinder_segment_field_internal(
 
     # case2a: hollow cylinder <- should be vectorized together with above
     mask2 = (r1 != 0) & mask1x
-    BHfinal[mask2] -= magnet_cylinder_field(
+    BHfinal[mask2] -= BHJM_magnet_cylinder(
         field=field,
         observers=observers[mask2],
         polarization=polarization[mask2],
@@ -2346,72 +2360,19 @@ def magnet_cylinder_segment_field_internal(
     return BHfinal
 
 
-# CORE
-def magnet_cylinder_segment_field(
-    *,
+def BHJM_cylinder_segment(
     field: str,
     observers: np.ndarray,
     dimension: np.ndarray,
     polarization: np.ndarray,
 ) -> np.ndarray:
-    """Magnetic field of homogeneously magnetized cylinder segments.
-
-    The cylinder axis coincides with the z-axis of the global coordinate
-    system. The geometric center of the cylinder lies in the origin.
-
-    SI units are used for all inputs and outputs.
-
-    Parameters
-    ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
-    observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
-
-    dimension: ndarray, shape (n,5)
-        Dimensions of CylinderSegments (r1,r2,h,phi1,phi2) with inner radius
-        r1, outer radius r2, height h in units of m and the two segment
-        angles phi1 < phi2 in units of deg.
-
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T.
-
-    Returns
-    -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the field of two different cylinder segment magnets at position (1,1,1).
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.core.magnet_cylinder_segment_field(
-    ...     field='B',
-    ...     observers=np.array([(1,1,1), (1,1,1)]),
-    ...     dimension=np.array([(0,1,2,-90,90), (1,2,4,35,125)]),
-    ...     polarization=np.array([(0,0,1), (.5,.5,0)]),
-    ... )
-    >>> print(B)
-    [[ 0.07046526  0.08373724 -0.0198113 ]
-     [ 0.29846023  0.20757316  0.00349617]]
-
-    Notes
-    -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
-    Implementation based on F.Slanovc, Journal of Magnetism and Magnetic
-    Materials, Volume 559, 1 October 2022, 169482
     """
-    bh = check_field_input(field, "magnet_cylinder_segment_field()")
+    - translate cylinder segment field to BHJM
+    - special cases catching
+    """
+    check_field_input(field)
 
-    BHfinal = np.zeros((len(polarization), 3))
+    BHJM = polarization.astype(float)
 
     r1, r2, h, phi1, phi2 = dimension.T
     r1 = abs(r1)
@@ -2431,6 +2392,8 @@ def magnet_cylinder_segment_field(
 
     # determine when points lie inside and on surface of magnet -------------
 
+    # mask_inside = None
+    # if in_out == "auto":
     # phip1 in [-2pi,0], phio2 in [0,2pi]
     phio1 = phi
     phio2 = phi - np.sign(phi) * 2 * np.pi
@@ -2452,42 +2415,68 @@ def magnet_cylinder_segment_field(
     )  # top / bottom
     mask_surf_r = (close(r, r1) | close(r, r2)) & mask_phi_in & mask_z_in  # in / out
     mask_surf_phi = (mask_phi1 | mask_phi2) & mask_r_in & mask_z_in  # in / out
-    mask_on_surface = mask_surf_z | mask_surf_r | mask_surf_phi
-    mask_not_on_surf = ~mask_on_surface
+    mask_not_on_surf = ~(mask_surf_z | mask_surf_r | mask_surf_phi)
 
     # inside
     mask_inside = mask_r_in & mask_phi_in & mask_z_in
+    # else:
+    #     mask_inside = np.full(len(observers), in_out == "inside")
+    #     mask_not_on_surf = np.full(len(observers), True)
+    # ACHTUNG @alex
+    #   1. inside and not_on_surface are not the same! Cant just put to true.
 
-    # return 0 when all points are on surface --------------------------------
-    if np.all(mask_on_surface):
-        return BHfinal
+    # return 0 when all points are on surface
+    if not np.any(mask_not_on_surf):
+        return BHJM * 0
+
+    if field == "J":
+        BHJM[~mask_inside] = 0
+        return BHJM
+
+    if field == "M":
+        BHJM[~mask_inside] = 0
+        return BHJM / MU0
+
+    BHJM *= 0
 
     # redefine input if there are some surface-points -------------------------
-    magg = polarization[mask_not_on_surf]
+    pol = polarization[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]
 
     # transform mag to spherical CS -----------------------------------------
-    m = np.sqrt(magg[:, 0] ** 2 + magg[:, 1] ** 2 + magg[:, 2] ** 2)
-    phi_m = np.arctan2(magg[:, 1], magg[:, 0])
-    th_m = np.arctan2(np.sqrt(magg[:, 0] ** 2 + magg[:, 1] ** 2), magg[:, 2])
+    m = np.sqrt(pol[:, 0] ** 2 + pol[:, 1] ** 2 + pol[:, 2] ** 2) / MU0  # J -> M
+    phi_m = np.arctan2(pol[:, 1], pol[:, 0])
+    th_m = np.arctan2(np.sqrt(pol[:, 0] ** 2 + pol[:, 1] ** 2), pol[:, 2])
     mag_sph = np.concatenate(((m,), (phi_m,), (th_m,)), axis=0).T
 
     # compute H and transform to cart CS -------------------------------------
-    H_cy = magnet_cylinder_segment_core(mag_sph, dim, pos_obs_cy)
+    H_cy = magnet_cylinder_segment_Hfield(
+        magnetizations=mag_sph, dimensions=dim, observers=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)
-    H = np.concatenate(((Hx,), (Hy,), (Hz,)), axis=0).T
+    BHJM[mask_not_on_surf] = np.concatenate(((Hx,), (Hy,), (Hz,)), axis=0).T / MU0
 
-    # return B or H --------------------------------------------------------
-    if not bh:
-        BHfinal[mask_not_on_surf] = H / MU0
-        return BHfinal
+    if field == "H":
+        return BHJM
 
-    B = H
-    BHfinal[mask_not_on_surf] = B
-    maskX = mask_inside * mask_not_on_surf
-    BHfinal[maskX] += polarization[maskX]
-    return BHfinal
+    if field == "B":
+        BHJM *= MU0
+        BHJM[mask_inside] += polarization[mask_inside]
+        BHJM[~mask_not_on_surf] *= 0
+        return BHJM
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
+
+    # return convert_HBMJ(
+    #     output_field_type=field,
+    #     polarization=polarization,
+    #     input_field_type="H",
+    #     field_values=H_all,
+    #     mask_inside=mask_inside & mask_not_on_surf,
+    # )
diff --git a/magpylib/_src/fields/field_BH_dipole.py b/magpylib/_src/fields/field_BH_dipole.py
index 70cdd48d9..f5e72427c 100644
--- a/magpylib/_src/fields/field_BH_dipole.py
+++ b/magpylib/_src/fields/field_BH_dipole.py
@@ -2,86 +2,90 @@
 Core implementation of dipole field
 """
 import numpy as np
+from scipy.constants import mu_0 as MU0
 
 from magpylib._src.input_checks import check_field_input
-from magpylib._src.utility import MU0
 
 
 # CORE
-def dipole_field(
-    *,
-    field: str,
+def dipole_Hfield(
     observers: np.ndarray,
-    moment: np.ndarray,
+    moments: np.ndarray,
 ) -> np.ndarray:
     """Magnetic field of a dipole moments.
 
     The dipole moment lies in the origin of the coordinate system.
-
-    SI units are used for all inputs and outputs.
+    The output is proportional to the moment input, and is independent
+    of length units used for observers (and moment) input considering
+    that the moment is proportional to [L]**2.
+    Returns np.inf for all non-zero moment components in the origin.
 
     Parameters
     ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
     observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
+        Observer positions (x,y,z) in Cartesian coordinates.
 
-    moment: ndarray, shape (n,3)
-        Dipole moment vector in units of A·m².
+    moments: ndarray, shape (n,3)
+        Dipole moment vector.
 
     Returns
     -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the B-field of two different dipole-observer instances.
-
-    >>> import magpylib as magpy
-    >>> import numpy as np
-    >>> B = magpy.core.dipole_field(
-    ...     field="B",
-    ...     observers=np.array([(1,2,3), (-1,-2,-3)]),
-    ...     moment=np.array([(0,0,1e6), (1e5,0,1e7)])
-    ... )
-    >>> print(B)
-    [[0.00122722 0.00245444 0.00177265]
-     [0.01212221 0.02462621 0.01784923]]
+    H-field: ndarray, shape (n,3)
+        H-field of Dipole in Cartesian coordinates.
 
     Notes
     -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
     The moment of a magnet is given by its volume*magnetization.
     """
-    bh = check_field_input(field, "dipole_field()")
 
     x, y, z = observers.T
     r = np.sqrt(x**2 + y**2 + z**2)  # faster than np.linalg.norm
     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 * observers, axis=1) * observers.T / r**5
-            - moment.T / r**3
-        ).T * 1e-7
-
-    # when r=0 return np.inf in all non-zero moment directions
+        H = (
+            (
+                3 * np.sum(moments * observers, axis=1) * observers.T / r**5
+                - moments.T / r**3
+            ).T
+            / 4
+            / np.pi
+        )
+
+    # when r=0 return np.inf in all non-zero moments directions
     mask1 = r == 0
     if np.any(mask1):
         with np.errstate(divide="ignore", invalid="ignore"):
-            B[mask1] = moment[mask1] / 0.0
-            np.nan_to_num(B, copy=False, posinf=np.inf, neginf=np.NINF)
-
-    if bh:
-        return B
+            H[mask1] = moments[mask1] / 0.0
+            np.nan_to_num(H, copy=False, posinf=np.inf, neginf=np.NINF)
 
-    H = B / MU0
     return H
+
+
+def BHJM_dipole(
+    field: str,
+    observers: np.ndarray,
+    moment: np.ndarray,
+) -> np.ndarray:
+    """
+    - translate dipole field to BHJM
+    """
+    check_field_input(field)
+
+    if field in "MJ":
+        return np.zeros_like(observers, dtype=float)
+
+    BHJM = dipole_Hfield(
+        observers=observers,
+        moments=moment,
+    )
+
+    if field == "H":
+        return BHJM
+
+    if field == "B":
+        return BHJM * MU0
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
diff --git a/magpylib/_src/fields/field_BH_polyline.py b/magpylib/_src/fields/field_BH_polyline.py
index 2ca9bb7ab..afe029b6b 100644
--- a/magpylib/_src/fields/field_BH_polyline.py
+++ b/magpylib/_src/fields/field_BH_polyline.py
@@ -1,14 +1,11 @@
 """
 Implementations of analytical expressions of line current segments
 """
-import warnings
-
 import numpy as np
 from numpy.linalg import norm
+from scipy.constants import mu_0 as MU0
 
-from magpylib._src.exceptions import MagpylibDeprecationWarning
 from magpylib._src.input_checks import check_field_input
-from magpylib._src.utility import MU0
 
 
 def current_vertices_field(
@@ -34,7 +31,7 @@ def current_vertices_field(
     - B-field (ndarray nx3): B-field vectors at pos_obs in units of T
     """
     if vertices is None:
-        return current_polyline_field(
+        return BHJM_current_polyline(
             field=field,
             observers=observers,
             current=current,
@@ -45,7 +42,7 @@ def current_vertices_field(
     nvs = np.array([f.shape[0] for f in vertices])  # lengths of vertices sets
     if all(v == nvs[0] for v in nvs):  # if all vertices sets have the same lenghts
         n0, n1, *_ = vertices.shape
-        BH = current_polyline_field(
+        BH = BHJM_current_polyline(
             field=field,
             observers=np.repeat(observers, n1 - 1, axis=0),
             current=np.repeat(current, n1 - 1, axis=0),
@@ -56,7 +53,7 @@ def current_vertices_field(
         BH = np.sum(BH, axis=1)
     else:
         split_indices = np.cumsum(nvs - 1)[:-1]  # remove last to avoid empty split
-        BH = current_polyline_field(
+        BH = BHJM_current_polyline(
             field=field,
             observers=np.repeat(observers, nvs - 1, axis=0),
             current=np.repeat(current, nvs - 1, axis=0),
@@ -69,99 +66,44 @@ def current_vertices_field(
 
 
 # CORE
-def current_polyline_field(
-    *,
-    field: str,
+def current_polyline_Hfield(
     observers: np.ndarray,
-    segment_start: np.ndarray,
-    segment_end: np.ndarray,
-    current: np.ndarray,
+    segments_start: np.ndarray,
+    segments_end: np.ndarray,
+    currents: np.ndarray,
 ) -> np.ndarray:
-    """Magnetic field of line current segments.
+    """B-field of line current segments.
 
     The current flows from start to end positions. The field is set to (0,0,0) on a
-    line segment.
-
-    SI units are used for all inputs and outputs.
+    line segment. The output is proportional to the current and independent of the
+    length units chosen for observers and dimensions.
 
     Parameters
     ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
     observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
+        Observer positions (x,y,z) in Cartesian coordinates.
 
-    segment_start: ndarray, shape (n,3)
-        Polyline start positions (x,y,z) in Cartesian coordinates in units of m.
+    segments_start: ndarray, shape (n,3)
+        Polyline start positions (x,y,z) in Cartesian coordinates.
 
-    segment_end: ndarray, shape (n,3)
-        Polyline end positions (x,y,z) in Cartesian coordinates in units of m.
+    segments_end: ndarray, shape (n,3)
+        Polyline end positions (x,y,z) in Cartesian coordinates.
 
-    current: ndarray, shape (n,)
-        Electrical current in units of A.
+    currents: ndarray, shape (n,)
+        Electrical currents in segements.
 
     Returns
     -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/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
-    >>> H = magpy.core.current_polyline_field(
-    ...     field='H',
-    ...     observers=np.array([( 0,0,1)]*3),
-    ...     segment_start=np.array([(-.5,0,0), (.5,0,0), (1.5,0,0)]),
-    ...     segment_end=np.array([(-1.5,0,0), (-.5,0,0), (.5,0,0)]),
-    ...     current=np.array([1,1,1]),
-    ... )
-    >>> print(H)
-    [[ 0.          0.03062433 -0.        ]
-     [ 0.          0.07117625 -0.        ]
-     [ 0.          0.03062433 -0.        ]]
+    B-Field: ndarray, shape (n,3)
+        B-field generated by current segements at observer positions.
 
     Notes
     -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
     Field computation via law of Biot Savart. See also countless online resources.
     eg. http://www.phys.uri.edu/gerhard/PHY204/tsl216.pdf
     """
-    # pylint: disable=too-many-statements
-    bh = check_field_input(field, "current_polyline_field()")
-
-    # allocate for special case treatment
-    ntot = len(current)
-    field_all = np.zeros((ntot, 3))
-
-    # Check for zero-length segments (or discontinuous)
-    mask_nan_start = np.isnan(segment_start).all(axis=1)
-    mask_nan_end = np.isnan(segment_end).all(axis=1)
-    mask_equal = np.all(segment_start == segment_end, axis=1)
-    mask0 = mask_equal | mask_nan_start | mask_nan_end
-
-    if np.all(mask0):
-        return field_all
-
-    # continue only with non-zero segments
-    if np.any(mask0):
-        not_mask0 = ~mask0  # avoid multiple computation of ~mask
-        current = current[not_mask0]
-        segment_start = segment_start[not_mask0]
-        segment_end = segment_end[not_mask0]
-        observers = observers[not_mask0]
-
     # rename
-    p1, p2, po = segment_start, segment_end, observers
+    p1, p2, po = segments_start, segments_end, observers
 
     # make dimensionless (avoid all large/small input problems) by introducing
     # the segment length as characteristic length scale.
@@ -179,8 +121,6 @@ def current_polyline_field(
 
     # separate on-line cases (-> B=0)
     mask1 = norm_o4 < 1e-15  # account for numerical issues
-    if np.all(mask1):
-        return field_all
 
     # continue only with general off-line cases
     if np.any(mask1):
@@ -191,7 +131,7 @@ def current_polyline_field(
         p4 = p4[not_mask1]
         norm_12 = norm_12[not_mask1]
         norm_o4 = norm_o4[not_mask1]
-        current = current[not_mask1]
+        currents = currents[not_mask1]
 
     # determine field direction
     cros = np.cross(p2 - p1, po - p4)
@@ -220,29 +160,66 @@ def current_polyline_field(
     mask4 = ~mask2 * ~mask3
     deltaSin[mask4] = abs(sinTh1[mask4] + sinTh2[mask4])
 
-    field = (deltaSin / norm_o4 * eB.T / norm_12 * current * 1e-7).T
+    # B = (deltaSin / norm_o4 * eB.T / norm_12 * current * 1e-7).T
+
+    # avoid array creation if possible
+    if np.any(mask1):
+        H = np.zeros_like(observers, dtype=float)
+        H[~mask1] = (deltaSin / norm_o4 * eB.T / norm_12 * currents / (4 * np.pi)).T
+        return H
+    return (deltaSin / norm_o4 * eB.T / norm_12 * currents / (4 * np.pi)).T
+
+
+def BHJM_current_polyline(
+    field: str,
+    observers: np.ndarray,
+    segment_start: np.ndarray,
+    segment_end: np.ndarray,
+    current: np.ndarray,
+) -> np.ndarray:
+    """
+    - translate Polyline field to BHJM
+    - treat some special cases
+    """
+    # pylint: disable=too-many-statements
+
+    check_field_input(field)
+
+    BHJM = np.zeros_like(observers, dtype=float)
+
+    if field in "MJ":
+        return BHJM
 
-    # broadcast general case results into allocated vector
-    mask0[~mask0] = mask1
-    field_all[~mask0] = field
+    # Check for zero-length segments (or discontinuous)
+    mask_nan_start = np.isnan(segment_start).all(axis=1)
+    mask_nan_end = np.isnan(segment_end).all(axis=1)
+    mask_equal = np.all(segment_start == segment_end, axis=1)
+    mask0 = mask_equal | mask_nan_start | mask_nan_end
+    not_mask0 = ~mask0  # avoid multiple computation of ~mask
+
+    if np.all(mask0):
+        return BHJM
 
-    # return B
-    if bh:
-        return field_all
+    # continue only with non-zero segments
+    if np.any(mask0):
+        current = current[not_mask0]
+        segment_start = segment_start[not_mask0]
+        segment_end = segment_end[not_mask0]
+        observers = observers[not_mask0]
 
-    # return H
-    return field_all / MU0
+    BHJM[not_mask0] = current_polyline_Hfield(
+        observers=observers,
+        segments_start=segment_start,
+        segments_end=segment_end,
+        currents=current,
+    )
 
+    if field == "H":
+        return BHJM
 
-def current_line_field(*args, **kwargs):
-    """current_line_field is deprecated, see current_polyline_field"""
+    if field == "B":
+        return BHJM * MU0
 
-    warnings.warn(
-        (
-            "current_line_field is deprecated and will be removed in a future version, "
-            "use current_polyline_field instead."
-        ),
-        MagpylibDeprecationWarning,
-        stacklevel=2,
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
     )
-    return current_polyline_field(*args, **kwargs)
diff --git a/magpylib/_src/fields/field_BH_sphere.py b/magpylib/_src/fields/field_BH_sphere.py
index 39d261cf0..af7090a3a 100644
--- a/magpylib/_src/fields/field_BH_sphere.py
+++ b/magpylib/_src/fields/field_BH_sphere.py
@@ -3,103 +3,101 @@
 magnetized Spheres. Computation details in function docstrings.
 """
 import numpy as np
+from scipy.constants import mu_0 as MU0
 
 from magpylib._src.input_checks import check_field_input
-from magpylib._src.utility import MU0
 
 
 # CORE
-def magnet_sphere_field(
-    *,
-    field: str,
+def magnet_sphere_Bfield(
     observers: np.ndarray,
-    diameter: np.ndarray,
-    polarization: np.ndarray,
+    diameters: np.ndarray,
+    polarizations: np.ndarray,
 ) -> np.ndarray:
-    """Magnetic field of homogeneously magnetized spheres.
+    """Magnetic field of homogeneously magnetized spheres in Cartesian Coordinates.
 
-    The center of the sphere lies in the origin of the coordinate system.
-
-    SI units are used for all inputs and outputs.
+    The center of the spheres lie in the origin of the coordinate system. The output
+    is proportional to the polarization magnitudes, and independent of the length
+    units chosen for observers and diameters.
 
     Parameters
     ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
     observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
+        Observer positions (x,y,z) in Cartesian coordinates.
 
-    diameter: ndarray, shape (n,)
-        Sphere diameters in units of m.
+    diameters: ndarray, shape (n,)
+        Sphere diameters.
 
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T.
+    polarizations: ndarray, shape (n,3)
+        Magnetic polarization vectors.
 
     Returns
     -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the field of two different spherical magnets at position (1,1,1),
-    inside and outside.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.core.magnet_sphere_field(
-    ...     field='B',
-    ...     observers=np.array([(1,1,1), (1,1,1)]),
-    ...     diameter=np.array([1,5]),
-    ...     polarization=np.array([(1,2,3), (0,0,3)]),
-    ... )
-    >>> print(B)
-    [[0.04009377 0.03207501 0.02405626]
-     [0.         0.         2.        ]]
+    B-field: ndarray, shape (n,3)
+        B-field generated by Spheres at observer positions.
 
     Notes
     -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
     The field corresponds to a dipole field on the outside and is 2/3*mag
     in the inside (see e.g. "Theoretical Physics, Bertelmann").
     """
+    return BHJM_magnet_sphere(
+        field="B",
+        observers=observers,
+        diameter=diameters,
+        polarization=polarizations,
+    )
 
-    bh = check_field_input(field, "magnet_sphere_field()")
 
-    # all special cases r0=0 and mag=0 automatically covered
+def BHJM_magnet_sphere(
+    field: str,
+    observers: np.ndarray,
+    diameter: np.ndarray,
+    polarization: np.ndarray,
+) -> np.ndarray:
+    """
+    - compute sphere field and translate to BHJM
+    - magnet sphere field, cannot be moved to a core function, because
+    core computation requires inside-outside check, but BHJM translation also.
+    Would require 2 checks, or forwarding the masks ... both not ideal
+    """
+    check_field_input(field)
 
     x, y, z = np.copy(observers.T)
     r = np.sqrt(x**2 + y**2 + z**2)  # faster than np.linalg.norm
-    r0 = abs(diameter) / 2
+    r_sphere = abs(diameter) / 2
 
     # inside field & allocate
-    B = polarization * 2 / 3
+    BHJM = polarization.astype(float)
+    out = r > r_sphere
+
+    if field == "J":
+        BHJM[out] = 0.0
+        return BHJM
+
+    if field == "M":
+        BHJM[out] = 0.0
+        return BHJM / MU0
+
+    BHJM *= 2 / 3
+
+    BHJM[out] = (
+        (
+            3 * np.sum(polarization[out] * observers[out], axis=1) * observers[out].T
+            - polarization[out].T * r[out] ** 2
+        )
+        / r[out] ** 5
+        * r_sphere[out] ** 3
+        / 3
+    ).T
 
-    # overwrite outside field entries
-    mask_out = r >= r0
+    if field == "B":
+        return BHJM
 
-    mag1 = polarization[mask_out]
-    obs1 = observers[mask_out]
-    r1 = r[mask_out]
-    r01 = r0[mask_out]
+    if field == "H":
+        BHJM[~out] -= polarization[~out]
+        return BHJM / MU0
 
-    field_out = (
-        (3 * (np.sum(mag1 * obs1, axis=1) * obs1.T) / r1**5 - mag1.T / r1**3)
-        * r01**3
-        / 3
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
     )
-    B[mask_out] = field_out.T
-
-    if bh:
-        return B
-
-    # adjust and return H
-    B[~mask_out] = -polarization[~mask_out] / 3
-    H = B / MU0
-    return H
diff --git a/magpylib/_src/fields/field_BH_tetrahedron.py b/magpylib/_src/fields/field_BH_tetrahedron.py
index 6a0f68458..f6d776162 100644
--- a/magpylib/_src/fields/field_BH_tetrahedron.py
+++ b/magpylib/_src/fields/field_BH_tetrahedron.py
@@ -3,8 +3,9 @@
 magnetized tetrahedra. Computation details in function docstrings.
 """
 import numpy as np
+from scipy.constants import mu_0 as MU0
 
-from magpylib._src.fields.field_BH_triangle import triangle_field
+from magpylib._src.fields.field_BH_triangle import BHJM_triangle
 from magpylib._src.input_checks import check_field_input
 
 
@@ -38,11 +39,17 @@ def check_chirality(points: np.ndarray) -> np.ndarray:
     return points
 
 
-def point_inside(points: np.ndarray, vertices: np.ndarray) -> np.ndarray:
+def point_inside(points: np.ndarray, vertices: np.ndarray, in_out: str) -> np.ndarray:
     """
     Takes points, as well as the vertices of a tetrahedra.
     Returns boolean array indicating whether the points are inside the tetrahedra.
     """
+    if in_out == "inside":
+        return np.array([True] * len(points))
+
+    if in_out == "outside":
+        return np.array([False] * len(points))
+
     mat = vertices[:, 1:].swapaxes(0, 1) - vertices[:, 0]
     mat = np.transpose(mat.swapaxes(0, 1), (0, 2, 1))
 
@@ -57,73 +64,36 @@ def point_inside(points: np.ndarray, vertices: np.ndarray) -> np.ndarray:
     return inside
 
 
-# CORE
-def magnet_tetrahedron_field(
-    *,
+def BHJM_magnet_tetrahedron(
     field: str,
     observers: np.ndarray,
     vertices: np.ndarray,
     polarization: np.ndarray,
+    in_out="auto",
 ) -> np.ndarray:
     """
-    Magnetic field generated by a homogeneously magnetized tetrahedra.
-
-    SI units are used for all inputs and outputs.
-
-    Parameters
-    ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
-    observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
-
-    vertices: ndarray, shape (n,4,3)
-        Vertices of the individual tetrahedrons [(pos1a, pos1b, pos1c, pos1d),
-        (pos2a, pos2b, pos2c, pos2d), ...] given in units of m.
+    - compute tetrahedron field from Triangle field
+    - translate to BHJM
+    - treat special cases
+    """
 
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T.
+    check_field_input(field)
 
-    Returns
-    -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the B-field of two different tetrahedron-observer instances
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.core.magnet_tetrahedron_field(
-    ...     field='B',
-    ...     observers=np.array([(1,2,3), (2,3,4)]),
-    ...     vertices=np.array([((-1,0,0), (1,-1,0), (1,1,0), (0,0,1))]*2),
-    ...     polarization=np.array([(222,333,444), (111,112,113)]),
-    ... )
-    >>> print(B)
-    [[0.19075398 0.8240532  1.18170862]
-     [0.03125701 0.08445416 0.1178967 ]]
-
-    Notes
-    -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
-    The tetrahedron is built up via 4 faces applying the Triangle class, making sure that
-    all normal vectors point outwards, and providing inside-outside evaluation to
-    distinguish between B- and H-field.
-    """
+    # allocate - try not to generate more arrays
+    BHJM = polarization.astype(float)
 
-    bh = check_field_input(field, "magnet_tetrahedron_field()")
+    if field == "J":
+        mask_inside = point_inside(observers, vertices, in_out)
+        BHJM[~mask_inside] = 0
+        return BHJM
 
-    n = len(observers)
+    if field == "M":
+        mask_inside = point_inside(observers, vertices, in_out)
+        BHJM[~mask_inside] = 0
+        return BHJM / MU0
 
     vertices = check_chirality(vertices)
+
     tri_vertices = np.concatenate(
         (
             vertices[:, (0, 2, 1), :],
@@ -133,23 +103,28 @@ def magnet_tetrahedron_field(
         ),
         axis=0,
     )
-    tri_fields = triangle_field(
+    tri_field = BHJM_triangle(
         field=field,
         observers=np.tile(observers, (4, 1)),
         vertices=tri_vertices,
         polarization=np.tile(polarization, (4, 1)),
     )
-    tetra_field = (  # slightly faster than reshape + sum
-        tri_fields[:n]
-        + tri_fields[n : 2 * n]
-        + tri_fields[2 * n : 3 * n]
-        + tri_fields[3 * n :]
+    n = len(observers)
+    BHJM = (  # slightly faster than reshape + sum
+        tri_field[:n]
+        + tri_field[n : 2 * n]
+        + tri_field[2 * n : 3 * n]
+        + tri_field[3 * n :]
     )
 
-    if not bh:
-        return tetra_field
+    if field == "H":
+        return BHJM
 
-    # if B, and inside magnet add polarization vector
-    mask_inside = point_inside(observers, vertices)
-    tetra_field[mask_inside] += polarization[mask_inside]
-    return tetra_field
+    if field == "B":
+        mask_inside = point_inside(observers, vertices, in_out)
+        BHJM[mask_inside] += polarization[mask_inside]
+        return BHJM
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
diff --git a/magpylib/_src/fields/field_BH_triangle.py b/magpylib/_src/fields/field_BH_triangle.py
index 0691043b9..289b938fa 100644
--- a/magpylib/_src/fields/field_BH_triangle.py
+++ b/magpylib/_src/fields/field_BH_triangle.py
@@ -4,9 +4,9 @@
 """
 # pylance: disable=Code is unreachable
 import numpy as np
+from scipy.constants import mu_0 as MU0
 
 from magpylib._src.input_checks import check_field_input
-from magpylib._src.utility import MU0
 
 
 def vcross3(a: np.ndarray, b: np.ndarray) -> np.ndarray:
@@ -72,77 +72,49 @@ def solid_angle(R: np.ndarray, r: np.ndarray) -> np.ndarray:
     return np.where(abs(result) > 6.2831853, 0, result)
 
 
-# CORE
-def triangle_field(
-    *,
-    field: str,
+def triangle_Bfield(
     observers: np.ndarray,
     vertices: np.ndarray,
-    polarization: np.ndarray,
+    polarizations: np.ndarray,
 ) -> np.ndarray:
     """Magnetic field generated by homogeneously magnetically charged triangular surfaces.
 
-    SI units are used for all inputs and outputs.
+    The charge is proportional to the projection of the polarization vectors onto the
+    triangle surfaces. The order of the triangle vertices defines the sign of the
+    surface normal vector (right-hand-rule). The output is proportional to the
+    polarization magnitude, and independent of the length units chosen for observers
+    and vertices.
 
-    Can be used to compute the field of a homogeneously magnetized bodies
-    with triangular surface mesh. In this case each Triangle must be defined so that the
-    normal vector points outwards.
+    Can be used to compute the field of a homogeneously magnetized bodies with triangular
+    surface mesh. In this case each Triangle must be defined so that the surface normal
+    vector points outwards.
 
     Parameters
     ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
     observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
+        Observer positions (x,y,z) in Cartesian coordinates.
 
-    dimensions: ndarray, shape (n,3)
-        Length of Cuboid sides in units of m.
+    vertices: ndarray, shape (n,3,3)
+        Triangle vertex positions ((P11,P12,P13), (P21, P22, P23), ...) in Cartesian
+        coordinates.
 
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T. The triangle surface charge is the
-        projection of the polarization vector onto the surface normal vector. The order
-        of the vertices defines the sign of the normal vector (right-hand-rule).
+    polarizations: ndarray, shape (n,3)
+        Magnetic polarization vectors.
 
     Returns
     -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Examples
-    --------
-    Compute the B-field of two different triangle-observer instances.
-
-    >>> import numpy as np
-    >>> import magpylib as magpy
-    >>> B = magpy.core.triangle_field(
-    ...     field='B',
-    ...     observers=np.array([(-.1,.2,.1), (.1,.2,.1)]),
-    ...     vertices=np.array([((-1,0,0), (1,-1,0), (1,1,0))]*2),
-    ...     polarization=np.array([(.22,.33,.44), (.33,.44,.55)]),
-    ... )
-    >>> print(B)
-    [[-0.0548087   0.05350955  0.17683832]
-     [-0.04252323  0.05292106  0.23092368]]
+    B-field: ndarray, shape (n,3)
+        B-field generated by Triangles at observer positions.
 
     Notes
     -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
     Field computations implemented from Guptasarma, Geophysics, 1999, 64:1, 70-74.
     Corners give (nan, nan, nan). Edges and in-plane perp components are set to 0.
     Loss of precision when approaching a triangle as (x-edge)**2 :(
     Loss of precision with distance from the triangle as distance**3 :(
     """
-    # pylint: disable=too-many-statements
-    bh = check_field_input(field, "triangle_field()")
-
     n = norm_vector(vertices)
-    sigma = np.einsum("ij, ij->i", n, polarization)  # vectorized inner product
+    sigma = np.einsum("ij, ij->i", n, polarizations)  # vectorized inner product
 
     # vertex <-> observer
     R = np.swapaxes(vertices, 0, 1) - observers
@@ -186,10 +158,46 @@ def triangle_field(
         )
     PQR = np.einsum("ij, ijk -> jk", I, L)
     B = sigma * (n.T * solid_angle(R, r) - vcross3(n, PQR).T)
+    B = B / np.pi / 4.0
+
+    return B.T
 
-    # return B or compute and return H -------------
-    if bh:
-        return B.T / np.pi / 4.0
 
-    H = B.T / 4 / np.pi / MU0
-    return H
+def BHJM_triangle(
+    field: str,
+    observers: np.ndarray,
+    vertices: np.ndarray,
+    polarization: np.ndarray,
+) -> np.ndarray:
+    """
+    - translate triangle core field to BHJM
+    """
+    check_field_input(field)
+
+    BHJM = polarization.astype(float) * 0.0
+
+    if field == "M":
+        return BHJM
+
+    if field == "J":
+        return BHJM
+
+    BHJM = triangle_Bfield(
+        observers=observers,
+        vertices=vertices,
+        polarizations=polarization,
+    )
+
+    # new MU0 problem:
+    #   input is polarization -> output has MU0 on it and must be B
+    #   H will then be connected via MU0
+
+    if field == "B":
+        return BHJM
+
+    if field == "H":
+        return BHJM / MU0
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
diff --git a/magpylib/_src/fields/field_BH_triangularmesh.py b/magpylib/_src/fields/field_BH_triangularmesh.py
index bd0de6caf..9f2da9ce1 100644
--- a/magpylib/_src/fields/field_BH_triangularmesh.py
+++ b/magpylib/_src/fields/field_BH_triangularmesh.py
@@ -3,12 +3,13 @@
 Computation details in function docstrings.
 """
 # pylint: disable=too-many-nested-blocks
+# pylint: disable=too-many-branches
 # pylance: disable=Code is unreachable
 import numpy as np
 import scipy.spatial
+from scipy.constants import mu_0 as MU0
 
-from magpylib._src.fields.field_BH_triangle import triangle_field
-from magpylib._src.utility import MU0
+from magpylib._src.fields.field_BH_triangle import BHJM_triangle
 
 
 def calculate_centroid(vertices, faces):
@@ -494,9 +495,7 @@ def mask_inside_trimesh(points: np.ndarray, faces: np.ndarray) -> np.ndarray:
     return mask_inside
 
 
-# CORE LIKE - but is not a core function!
-def magnet_trimesh_field(
-    *,
+def BHJM_magnet_trimesh(
     field: str,
     observers: np.ndarray,
     mesh: np.ndarray,
@@ -504,103 +503,75 @@ def magnet_trimesh_field(
     in_out="auto",
 ) -> np.ndarray:
     """
-    Core-like function that computes the field of triangular meshes using the triangle_field
-
-    !!!Closed meshes are assumed (input comes only from TriangularMesh class)!!!
-    This is the reasons that this is not a core function
-
-    SI units are used for all inputs and outputs.
-
-    Parameters
-    ----------
-    field: str, default=`'B'`
-        If `field='B'` return B-field in units of T, if `field='H'` return H-field
-        in units of A/m.
-
-    observers: ndarray, shape (n,3)
-        Observer positions (x,y,z) in Cartesian coordinates in units of m.
-
-    mesh: ndarray, shape (n,n1,3,3) or ragged sequence
-        Triangular mesh of shape [(x1,y1,z1), (x2,y2,z2), (x3,y3,z3)].
-        `mesh` can be a ragged sequence of mesh-children with different lengths.
-
-    polarization: ndarray, shape (n,3)
-        Magnetic polarization vectors in units of T.
-
-    in_out: {'auto', 'inside', 'outside'}
-        Tells if the points are inside or outside the enclosing mesh for the correct B/H-field
-        calculation. By default `in_out='auto'` and the inside/outside mask is automatically
-        generated using a ray tracing algorithm to determine which observers are inside and which
-        are outside the closed body. For performance reasons, one can define `in_out='outside'`
-        or `in_out='inside'` if it is known in advance that all observers satisfy the same
-        condition.
-
-    Returns
-    -------
-    B-field or H-field: ndarray, shape (n,3)
-        B- or H-field of source in Cartesian coordinates in units of T or A/m.
-
-    Notes
-    -----
-    Advanced unit use: The input unit of magnetization and polarization
-    gives the output unit of H and B. All results are independent of the
-    length input units. One must be careful, however, to use consistently
-    the same length unit throughout a script.
-
-    Field computations via publication:
-    Guptasarma: GEOPHYSICS 1999 64:1, 70-74
+    - Compute triangular mesh field from triangle fields.
+    - Closed meshes are assumed (input comes only from TriangularMesh class)
+    - Field computations via publication: Guptasarma: GEOPHYSICS 1999 64:1, 70-74
     """
-    if mesh.ndim != 1:  # all vertices objects have same number of children
-        n0, n1, *_ = mesh.shape
-        vertices_tiled = mesh.reshape(-1, 3, 3)
-        observers_tiled = np.repeat(observers, n1, axis=0)
-        polarization_tiled = np.repeat(polarization, n1, axis=0)
-        B = triangle_field(
-            field="B",
-            observers=observers_tiled,
-            vertices=vertices_tiled,
-            polarization=polarization_tiled,
-        )
-        B = B.reshape((n0, n1, 3))
-        B = np.sum(B, axis=1)
+    if field in "BH":
+        if mesh.ndim != 1:  # all vertices objects have same number of children
+            n0, n1, *_ = mesh.shape
+            vertices_tiled = mesh.reshape(-1, 3, 3)
+            observers_tiled = np.repeat(observers, n1, axis=0)
+            polarization_tiled = np.repeat(polarization, n1, axis=0)
+            BHJM = BHJM_triangle(
+                field="B",
+                observers=observers_tiled,
+                vertices=vertices_tiled,
+                polarization=polarization_tiled,
+            )
+            BHJM = BHJM.reshape((n0, n1, 3))
+            BHJM = np.sum(BHJM, axis=1)
+        else:
+            nvs = [f.shape[0] for f in mesh]  # length of vertex set
+            split_indices = np.cumsum(nvs)[:-1]  # remove last to avoid empty split
+            vertices_tiled = np.concatenate([f.reshape((-1, 3, 3)) for f in mesh])
+            observers_tiled = np.repeat(observers, nvs, axis=0)
+            polarization_tiled = np.repeat(polarization, nvs, axis=0)
+            BHJM = BHJM_triangle(
+                field="B",
+                observers=observers_tiled,
+                vertices=vertices_tiled,
+                polarization=polarization_tiled,
+            )
+            b_split = np.split(BHJM, split_indices)
+            BHJM = np.array([np.sum(bh, axis=0) for bh in b_split])
     else:
-        nvs = [f.shape[0] for f in mesh]  # length of vertex set
-        split_indices = np.cumsum(nvs)[:-1]  # remove last to avoid empty split
-        vertices_tiled = np.concatenate([f.reshape((-1, 3, 3)) for f in mesh])
-        observers_tiled = np.repeat(observers, nvs, axis=0)
-        polarization_tiled = np.repeat(polarization, nvs, axis=0)
-        B = triangle_field(
-            field="B",
-            observers=observers_tiled,
-            vertices=vertices_tiled,
-            polarization=polarization_tiled,
-        )
-        b_split = np.split(B, split_indices)
-        B = np.array([np.sum(bh, axis=0) for bh in b_split])
+        BHJM = np.zeros_like(observers, dtype=float)
+
+    if field == "H":
+        return BHJM / MU0
+
+    if in_out == "auto":
+        prev_ind = 0
+        # group similar meshes for inside-outside evaluation and adding B
+        for new_ind, _ in enumerate(BHJM):
+            if (
+                new_ind == len(BHJM) - 1
+                or mesh[new_ind].shape != mesh[prev_ind].shape
+                or not np.all(mesh[new_ind] == mesh[prev_ind])
+            ):
+                if new_ind == len(BHJM) - 1:
+                    new_ind = len(BHJM)
+                mask_inside = mask_inside_trimesh(
+                    observers[prev_ind:new_ind], mesh[prev_ind]
+                )
+                # if inside magnet add polarization vector
+                BHJM[prev_ind:new_ind][mask_inside] += polarization[prev_ind:new_ind][
+                    mask_inside
+                ]
+                prev_ind = new_ind
+    elif in_out == "inside":
+        BHJM += polarization
 
     if field == "B":
-        if in_out == "auto":
-            prev_ind = 0
-            # group similar meshes for inside-outside evaluation and adding B
-            for new_ind, _ in enumerate(B):
-                if (
-                    new_ind == len(B) - 1
-                    or mesh[new_ind].shape != mesh[prev_ind].shape
-                    or not np.all(mesh[new_ind] == mesh[prev_ind])
-                ):
-                    if new_ind == len(B) - 1:
-                        new_ind = len(B)
-                    inside_mask = mask_inside_trimesh(
-                        observers[prev_ind:new_ind], mesh[prev_ind]
-                    )
-                    # if inside magnet add polarization vector
-                    B[prev_ind:new_ind][inside_mask] += polarization[prev_ind:new_ind][
-                        inside_mask
-                    ]
-                    prev_ind = new_ind
-        elif in_out == "inside":
-            B += polarization
-        return B
-
-    H = B / MU0
-    return H
+        return BHJM
+
+    if field == "J":
+        return BHJM
+
+    if field == "M":
+        return BHJM / MU0
+
+    raise ValueError(  # pragma: no cover
+        "`output_field_type` must be one of ('B', 'H', 'M', 'J'), " f"got {field!r}"
+    )
diff --git a/magpylib/_src/fields/field_wrap_BH.py b/magpylib/_src/fields/field_wrap_BH.py
index a266eefcd..79c94358c 100644
--- a/magpylib/_src/fields/field_wrap_BH.py
+++ b/magpylib/_src/fields/field_wrap_BH.py
@@ -1,9 +1,16 @@
+# pylint: disable=cyclic-import
+# pylint: disable=too-many-lines
 """Field computation structure:
 
-level0:(field_BH_XXX.py files)
+level_core:(field_BH_XXX.py files)
     - pure vectorized field computations from literature
     - all computations in source CS
-    - distinguish B/H
+
+level0a:(BHJM_XX)
+    - distinguish between B, H, J and M
+
+level0b:(BHJM_internal_XX)
+    - connect BHJM-level to level1
 
 level1(getBH_level1):
     - apply transformation to global CS
@@ -40,7 +47,6 @@
 
 level5(sens.getB, sens.getH): <--- USER INTERFACE
 """
-# pylint: disable=cyclic-import
 import numbers
 import warnings
 from itertools import product
@@ -60,6 +66,7 @@
 from magpylib._src.utility import format_obj_input
 from magpylib._src.utility import format_src_inputs
 from magpylib._src.utility import get_registered_sources
+from magpylib._src.utility import has_parameter
 
 
 def tile_group_property(group: list, n_pp: int, prop_name: str):
@@ -138,6 +145,10 @@ def getBH_level1(
     # transform obs_pos into source CS
     pos_rel_rot = orientation.apply(observers - position, inverse=True)
 
+    # filter arguments
+    if not has_parameter(field_func, "in_out"):  # in_out passed only to magnets
+        kwargs.pop("in_out", None)
+
     # compute field
     BH = field_func(field=field, observers=pos_rel_rot, **kwargs)
 
@@ -149,38 +160,9 @@ def getBH_level1(
 
 
 def getBH_level2(
-    sources, observers, *, field, sumup, squeeze, pixel_agg, output, **kwargs
+    sources, observers, *, field, sumup, squeeze, pixel_agg, output, in_out, **kwargs
 ) -> np.ndarray:
     """Compute field for given sources and observers.
-
-    Parameters
-    ----------
-    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) sfor all sources.
-    squeeze : bool, default=True:
-        If True output is squeezed (axes of length 1 are eliminated)
-    pixel_agg : str
-        A compatible numpy aggregator string (e.g. `'min'`, `'max'`, `'mean'`)
-        which applies on pixel output values.
-    field : {'B', 'H'}
-        'B' computes B field, 'H' computes H-field
-    output: str, default='ndarray'
-        Output type, which must be one of `('ndarray', 'dataframe')`. By default a multi-
-        dimensional array ('ndarray') is returned. If 'dataframe' is chosen, the function
-        returns a 2D-table as a `pandas.DataFrame` object (the Pandas library must be
-        installed).
-
-    Returns
-    -------
-    field: ndarray, shape squeeze((L,M,K,N1,N2,...,3)), field of L sources, M path
-    positions, K sensors and N1xN2x.. observer positions and 3 field components.
-
     Info:
     -----
     - generates a 1D list of sources (collections flattened) and a 1D list of sensors from input
@@ -205,6 +187,7 @@ def getBH_level2(
             observers=observers,
             field=field,
             squeeze=squeeze,
+            in_out=in_out,
             **kwargs,
         )
 
@@ -225,6 +208,19 @@ def getBH_level2(
     check_dimensions(src_list)
     check_excitations(src_list)
 
+    # make sure that given in_out there is a Tetrahedron class or a TriangularMesh
+    #   class in sources. Else throw a warning
+    if in_out != "auto":
+        from magpylib._src.obj_classes.class_magnet_Tetrahedron import Tetrahedron
+
+        if not any(isinstance(src, (Tetrahedron, TriangularMesh)) for src in src_list):
+            warnings.warn(
+                "Argument `in_out` for field computation was set, but is ignored"
+                " in the computation. `in_out` has an effect only for magnet classes"
+                " Tetrahedron and TriangularMesh.",
+                UserWarning,
+            )
+
     # make sure that TriangularMesh sources have a closed mesh when getB is called - warn if not
     if field == "B":
         for src in src_list:
@@ -338,9 +334,10 @@ def getBH_level2(
         lg = len(group["sources"])
         gr = group["sources"]
         src_dict = get_src_dict(gr, n_pix, n_pp, poso)  # compute array dict for level1
+        # compute field
         B_group = getBH_level1(
-            field_func=field_func, field=field, **src_dict
-        )  # compute field
+            field_func=field_func, field=field, in_out=in_out, **src_dict
+        )
         if B_group is None:
             raise MagpylibMissingInput(
                 f"Cannot compute {field}-field because "
@@ -450,6 +447,7 @@ def getBH_dict_level2(
     position=(0, 0, 0),
     orientation=R.identity(),
     squeeze=True,
+    in_out="auto",
     **kwargs: dict,
 ) -> np.ndarray:
     """Functional interface access to vectorized computation
@@ -546,7 +544,7 @@ def getBH_dict_level2(
     kwargs["orientation"] = R.from_quat(kwargs["orientation"])
 
     # compute and return B
-    B = getBH_level1(field=field, field_func=field_func, **kwargs)
+    B = getBH_level1(field=field, field_func=field_func, in_out=in_out, **kwargs)
 
     if B is not None and squeeze:
         return np.squeeze(B)
@@ -560,6 +558,7 @@ def getB(
     squeeze=True,
     pixel_agg=None,
     output="ndarray",
+    in_out="auto",
     **kwargs,
 ):
     """Compute B-field in units of T for given sources and observers.
@@ -604,6 +603,19 @@ def getB(
         `numpy.ndarray` object is returned. If `'dataframe'` is chosen, a `pandas.DataFrame`
         object is returned (the Pandas library must be installed).
 
+    in_out: {'auto', 'inside', 'outside'}
+        This parameter only applies for magnet bodies. It specifies the location of the
+        observers relative to the magnet body, affecting the calculation of the magnetic field.
+        The options are:
+        - 'auto': The location (inside or outside the cuboid) is determined automatically for
+        each observer.
+        - 'inside': All observers are considered to be inside the cuboid; use this for
+            performance optimization if applicable.
+        - 'outside': All observers are considered to be outside the cuboid; use this for
+            performance optimization if applicable.
+        Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+        locations is unknown.
+
     See Also
     --------
     *Functional interface
@@ -681,8 +693,8 @@ def getB(
     >>> src2 = magpy.magnet.Sphere(polarization=(0,0,.1), diameter=.001)
     >>> B = magpy.getB([src1, src2], (.01,.01,.01))
     >>> print(B)
-    [[6.05434592e-06 6.05434592e-06 2.35680448e-08]
-     [8.01875374e-07 8.01875374e-07 1.51582450e-22]]
+    [[ 6.05434592e-06  6.05434592e-06  2.35680448e-08]
+     [ 8.01875374e-07  8.01875374e-07 -9.05619815e-23]]
 
     We can also use sensor objects as observers input:
 
@@ -693,8 +705,8 @@ def getB(
     [[[ 6.05434592e-06  6.05434592e-06  2.35680448e-08]
       [-6.05434592e-06 -6.05434592e-06  2.35680448e-08]]
     <BLANKLINE>
-     [[ 8.01875374e-07  8.01875374e-07  1.51582450e-22]
-      [-8.01875374e-07 -8.01875374e-07  1.51582450e-22]]]
+     [[ 8.01875374e-07  8.01875374e-07 -9.05619815e-23]
+      [-8.01875374e-07 -8.01875374e-07 -9.05619815e-23]]]
 
     Through the functional interface we can compute the same fields for the loop as:
 
@@ -715,18 +727,19 @@ def getB(
     ... )
     >>> print(B)
     [[ 1.66322588e-07  1.66322588e-07  1.61742625e-10]
-     [-4.69451597e-07 -4.69451597e-07  4.70690813e-07]
+     [-4.69451598e-07 -4.69451598e-07  4.70690813e-07]
      [ 7.96993186e-07  1.59398637e-06 -7.91258466e-07]
      [-1.37369334e-06 -1.37369334e-06 -1.36554287e-06]]
     """
     return getBH_level2(
         sources,
         observers,
+        field="B",
         sumup=sumup,
         squeeze=squeeze,
         pixel_agg=pixel_agg,
         output=output,
-        field="B",
+        in_out=in_out,
         **kwargs,
     )
 
@@ -738,6 +751,7 @@ def getH(
     squeeze=True,
     pixel_agg=None,
     output="ndarray",
+    in_out="auto",
     **kwargs,
 ):
     """Compute H-field in units of A/m for given sources and observers.
@@ -782,6 +796,19 @@ def getH(
         `numpy.ndarray` object is returned. If `'dataframe'` is chosen, a `pandas.DataFrame`
         object is returned (the Pandas library must be installed).
 
+    in_out: {'auto', 'inside', 'outside'}
+        This parameter only applies for magnet bodies. It specifies the location of the
+        observers relative to the magnet body, affecting the calculation of the magnetic field.
+        The options are:
+        - 'auto': The location (inside or outside the cuboid) is determined automatically for
+        each observer.
+        - 'inside': All observers are considered to be inside the cuboid; use this for
+            performance optimization if applicable.
+        - 'outside': All observers are considered to be outside the cuboid; use this for
+            performance optimization if applicable.
+        Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+        locations is unknown.
+
     See Also
     --------
     *Functional interface
@@ -852,8 +879,8 @@ def getH(
     >>> src2 = magpy.magnet.Sphere(polarization=(0,0,.1), diameter=.001)
     >>> H = magpy.getH([src1, src2], (.01,.01,.01))
     >>> print(H)
-    [[4.81789540e+00 4.81789540e+00 1.87548541e-02]
-     [6.38112147e-01 6.38112147e-01 1.20625481e-16]]
+    [[ 4.81789540e+00  4.81789540e+00  1.87548541e-02]
+     [ 6.38112147e-01  6.38112147e-01 -7.20669350e-17]]
 
     We can also use sensor objects as observers input:
 
@@ -864,8 +891,8 @@ def getH(
     [[[ 4.81789540e+00  4.81789540e+00  1.87548541e-02]
       [-4.81789540e+00 -4.81789540e+00  1.87548541e-02]]
     <BLANKLINE>
-     [[ 6.38112147e-01  6.38112147e-01  1.20625481e-16]
-      [-6.38112147e-01 -6.38112147e-01  1.20625481e-16]]]
+     [[ 6.38112147e-01  6.38112147e-01 -7.20669350e-17]
+      [-6.38112147e-01 -6.38112147e-01 -7.20669350e-17]]]
 
     Through the functional interface we can compute the same fields for the loop as:
 
@@ -893,10 +920,288 @@ def getH(
     return getBH_level2(
         sources,
         observers,
+        field="H",
         sumup=sumup,
         squeeze=squeeze,
         pixel_agg=pixel_agg,
         output=output,
-        field="H",
+        in_out=in_out,
+        **kwargs,
+    )
+
+
+def getM(
+    sources=None,
+    observers=None,
+    sumup=False,
+    squeeze=True,
+    pixel_agg=None,
+    output="ndarray",
+    in_out="auto",
+    **kwargs,
+):
+    """Compute M-field in units of A/m for given sources and observers.
+
+    Field implementations can be directly accessed (avoiding the object oriented
+    Magpylib interface) by providing a string input `sources=source_type`, array_like
+    positions as `observers` input, and all other necessary input parameters (see below)
+    as kwargs.
+
+    Parameters
+    ----------
+    sources: source and collection objects or 1D list thereof
+        Sources that generate the magnetic field. Can be a single source (or collection)
+        or a 1D list of l sources and/or collection objects.
+
+        Functional interface: input must be one of (`Cuboid`, `Cylinder`, `CylinderSegment`,
+        `Sphere`, `Dipole`, `Circle` or `Polyline`).
+
+    observers: array_like or (list of) `Sensor` objects
+        Can be array_like positions of shape (n1, n2, ..., 3) where the field
+        should be evaluated, a `Sensor` object with pixel shape (n1, n2, ..., 3) or a list
+        of such sensor objects (must all have similar pixel shapes). All positions
+        are given in units of m.
+
+        Functional interface: Input must be array_like with shape (3,) or (n,3) corresponding
+        positions to observer positions in units of m.
+
+    sumup: bool, default=`False`
+        If `True`, the fields of all sources are summed up.
+
+    squeeze: bool, default=`True`
+        If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
+        a single sensor or only a single source) are eliminated.
+
+    pixel_agg: str, default=`None`
+        Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+        which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+        With this option, observers input with different (pixel) shapes is allowed.
+
+    output: str, default=`'ndarray'`
+        Output type, which must be one of (`'ndarray'`, `'dataframe'`). By default a
+        `numpy.ndarray` object is returned. If `'dataframe'` is chosen, a `pandas.DataFrame`
+        object is returned (the Pandas library must be installed).
+
+    in_out: {'auto', 'inside', 'outside'}
+        This parameter only applies for magnet bodies. It specifies the location of the
+        observers relative to the magnet body, affecting the calculation of the magnetic field.
+        The options are:
+        - 'auto': The location (inside or outside the cuboid) is determined automatically for
+        each observer.
+        - 'inside': All observers are considered to be inside the cuboid; use this for
+            performance optimization if applicable.
+        - 'outside': All observers are considered to be outside the cuboid; use this for
+            performance optimization if applicable.
+        Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+        locations is unknown.
+
+    See Also
+    --------
+    *Functional interface
+
+    position: array_like, shape (3,) or (n,3), default=`(0,0,0)`
+        Source position(s) in the global coordinates in units of m.
+
+    orientation: scipy `Rotation` object with length 1 or n, default=`None`
+        Object orientation(s) in the global coordinates. `None` corresponds to
+        a unit-rotation.
+
+    magnetization: array_like, shape (3,) or (n,3)
+        Only source_type in (`Cuboid`, `Cylinder`, `CylinderSegment`, `Sphere`)!
+        Magnetization vector(s) (mu0*M, remanence field) in units of A/m given in
+        the local object coordinates (rotates with object).
+
+    moment: array_like, shape (3) or (n,3), unit A·m²
+        Only source_type == `Dipole`!
+        Magnetic dipole moment(s) in units of A·m² given in the local object coordinates
+        (rotates with object). For homogeneous magnets the relation moment=magnetization*volume
+        holds.
+
+    current: array_like, shape (n,)
+        Only source_type == `Circle` or `Polyline`!
+        Electrical current in units of A.
+
+    dimension: array_like, shape (x,) or (n,x)
+        Only source_type in (`Cuboid`, `Cylinder`, `CylinderSegment`)!
+        Magnet dimension input in units of m and deg. Dimension format x of sources is similar
+        as in object oriented interface.
+
+    diameter: array_like, shape (n,)
+        Only source_type == `Sphere` or `Circle`!
+        Diameter of source in units of m.
+
+    segment_start: array_like, shape (n,3)
+        Only source_type == `Polyline`!
+        Start positions of line current segments in units of m.
+
+    segment_end: array_like, shape (n,3)
+        Only source_type == `Polyline`!
+        End positions of line current segments in units of m.
+
+    Returns
+    -------
+    M-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
+        M-field at each path position (index m) for each sensor (index k) and each sensor pixel
+        position (indices n1, n2, ...) in units of A/m. Sensor pixel positions are equivalent
+        to simple observer positions. Paths of objects that are shorter than index m are
+        considered as static beyond their end.
+
+    Functional interface: ndarray, shape (n,3)
+        M-field for every parameter set in units of A/m.
+
+    Notes
+    -----
+    This function automatically joins all sensor and position inputs together and groups
+    similar sources for optimal vectorization of the computation. For maximal performance
+    call this function as little as possible and avoid using it in loops.
+    """
+    return getBH_level2(
+        sources,
+        observers,
+        field="M",
+        sumup=sumup,
+        squeeze=squeeze,
+        pixel_agg=pixel_agg,
+        output=output,
+        in_out=in_out,
+        **kwargs,
+    )
+
+
+def getJ(
+    sources=None,
+    observers=None,
+    sumup=False,
+    squeeze=True,
+    pixel_agg=None,
+    output="ndarray",
+    in_out="auto",
+    **kwargs,
+):
+    """Compute J-field in units of T for given sources and observers.
+
+    Field implementations can be directly accessed (avoiding the object oriented
+    Magpylib interface) by providing a string input `sources=source_type`, array_like
+    positions as `observers` input, and all other necessary input parameters (see below)
+    as kwargs.
+
+    Parameters
+    ----------
+    sources: source and collection objects or 1D list thereof
+        Sources that generate the magnetic field. Can be a single source (or collection)
+        or a 1D list of l sources and/or collection objects.
+
+        Functional interface: input must be one of (`Cuboid`, `Cylinder`, `CylinderSegment`,
+        `Sphere`, `Dipole`, `Circle` or `Polyline`).
+
+    observers: array_like or (list of) `Sensor` objects
+        Can be array_like positions of shape (n1, n2, ..., 3) where the field
+        should be evaluated, a `Sensor` object with pixel shape (n1, n2, ..., 3) or a list
+        of such sensor objects (must all have similar pixel shapes). All positions
+        are given in units of m.
+
+        Functional interface: Input must be array_like with shape (3,) or (n,3) corresponding
+        positions to observer positions in units of m.
+
+    sumup: bool, default=`False`
+        If `True`, the fields of all sources are summed up.
+
+    squeeze: bool, default=`True`
+        If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
+        a single sensor or only a single source) are eliminated.
+
+    pixel_agg: str, default=`None`
+        Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+        which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+        With this option, observers input with different (pixel) shapes is allowed.
+
+    output: str, default=`'ndarray'`
+        Output type, which must be one of (`'ndarray'`, `'dataframe'`). By default a
+        `numpy.ndarray` object is returned. If `'dataframe'` is chosen, a `pandas.DataFrame`
+        object is returned (the Pandas library must be installed).
+
+    in_out: {'auto', 'inside', 'outside'}
+        This parameter only applies for magnet bodies. It specifies the location of the
+        observers relative to the magnet body, affecting the calculation of the magnetic field.
+        The options are:
+        - 'auto': The location (inside or outside the cuboid) is determined automatically for
+        each observer.
+        - 'inside': All observers are considered to be inside the cuboid; use this for
+            performance optimization if applicable.
+        - 'outside': All observers are considered to be outside the cuboid; use this for
+            performance optimization if applicable.
+        Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+        locations is unknown.
+
+    See Also
+    --------
+    *Functional interface
+
+    position: array_like, shape (3,) or (n,3), default=`(0,0,0)`
+        Source position(s) in the global coordinates in units of m.
+
+    orientation: scipy `Rotation` object with length 1 or n, default=`None`
+        Object orientation(s) in the global coordinates. `None` corresponds to
+        a unit-rotation.
+
+    magnetization: array_like, shape (3,) or (n,3)
+        Only source_type in (`Cuboid`, `Cylinder`, `CylinderSegment`, `Sphere`)!
+        Magnetization vector(s) (mu0*M, remanence field) in units of A/m given in
+        the local object coordinates (rotates with object).
+
+    moment: array_like, shape (3) or (n,3), unit A·m²
+        Only source_type == `Dipole`!
+        Magnetic dipole moment(s) in units of A·m² given in the local object coordinates
+        (rotates with object). For homogeneous magnets the relation moment=magnetization*volume
+        holds.
+
+    current: array_like, shape (n,)
+        Only source_type == `Circle` or `Polyline`!
+        Electrical current in units of A.
+
+    dimension: array_like, shape (x,) or (n,x)
+        Only source_type in (`Cuboid`, `Cylinder`, `CylinderSegment`)!
+        Magnet dimension input in units of m and deg. Dimension format x of sources is similar
+        as in object oriented interface.
+
+    diameter: array_like, shape (n,)
+        Only source_type == `Sphere` or `Circle`!
+        Diameter of source in units of m.
+
+    segment_start: array_like, shape (n,3)
+        Only source_type == `Polyline`!
+        Start positions of line current segments in units of m.
+
+    segment_end: array_like, shape (n,3)
+        Only source_type == `Polyline`!
+        End positions of line current segments in units of m.
+
+    Returns
+    -------
+    J-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
+        J-field at each path position (index m) for each sensor (index k) and each sensor pixel
+        position (indices n1, n2, ...) in units of T. Sensor pixel positions are equivalent
+        to simple observer positions. Paths of objects that are shorter than index m are
+        considered as static beyond their end.
+
+    Functional interface: ndarray, shape (n,3)
+        J-field for every parameter set in units of T.
+
+    Notes
+    -----
+    This function automatically joins all sensor and position inputs together and groups
+    similar sources for optimal vectorization of the computation. For maximal performance
+    call this function as little as possible and avoid using it in loops.
+
+    """
+    return getBH_level2(
+        sources,
+        observers,
+        field="J",
+        sumup=sumup,
+        squeeze=squeeze,
+        pixel_agg=pixel_agg,
+        output=output,
+        in_out=in_out,
         **kwargs,
     )
diff --git a/magpylib/_src/input_checks.py b/magpylib/_src/input_checks.py
index 06ab08183..229752585 100644
--- a/magpylib/_src/input_checks.py
+++ b/magpylib/_src/input_checks.py
@@ -118,17 +118,14 @@ def check_degree_type(inp):
         )
 
 
-def check_field_input(inp, origin):
+def check_field_input(inp):
     """check field input"""
-    if isinstance(inp, str):
-        if inp == "B":
-            return True
-        if inp == "H":
-            return False
-    raise MagpylibBadUserInput(
-        f"{origin} input can only be `field='B'` or `field='H'`.\n"
-        f"Instead received {repr(inp)}."
-    )
+    allowed = tuple("BHMJ")
+    if not (isinstance(inp, str) and inp in allowed):
+        raise MagpylibBadUserInput(
+            f"`field` input can only be one of {allowed}.\n"
+            f"Instead received {repr(inp)}."
+        )
 
 
 def validate_field_func(val):
diff --git a/magpylib/_src/obj_classes/class_BaseDisplayRepr.py b/magpylib/_src/obj_classes/class_BaseDisplayRepr.py
index b69b93e7b..1b404bd3a 100644
--- a/magpylib/_src/obj_classes/class_BaseDisplayRepr.py
+++ b/magpylib/_src/obj_classes/class_BaseDisplayRepr.py
@@ -48,7 +48,7 @@ def _get_description(self, exclude=None):
         for k in list(dict.fromkeys(list(UNITS) + list(params))):
             if not k.startswith("_") and k in params and k not in exclude:
                 unit = UNITS.get(k, None)
-                unit_str = f"{unit}" if unit else ""
+                unit_str = f" {unit}" if unit else ""
                 if k == "position":
                     val = getattr(self, "_position")
                     if val.shape[0] != 1:
@@ -79,7 +79,7 @@ def _get_description(self, exclude=None):
                         val = f"shape{val.shape}"
                 else:
                     val = getattr(self, k)
-                lines.append(f"  • {k}: {val} {unit_str}")
+                lines.append(f"  • {k}: {val}{unit_str}")
         return lines
 
     def describe(self, *, exclude=("style", "field_func"), return_string=False):
diff --git a/magpylib/_src/obj_classes/class_BaseExcitations.py b/magpylib/_src/obj_classes/class_BaseExcitations.py
index 20ce7f91e..18078ebb6 100644
--- a/magpylib/_src/obj_classes/class_BaseExcitations.py
+++ b/magpylib/_src/obj_classes/class_BaseExcitations.py
@@ -51,7 +51,9 @@ def field_func(self, val):
             )
         self._field_func = val
 
-    def getB(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
+    def getB(
+        self, *observers, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
         """Compute the B-field at observers in units of T generated by the source.
 
         SI units are used for all inputs and outputs.
@@ -79,6 +81,19 @@ def getB(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
             returns a 2D-table as a `pandas.DataFrame` object (the Pandas library must be
             installed).
 
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
         Returns
         -------
         B-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
@@ -115,14 +130,17 @@ def getB(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
         return getBH_level2(
             self,
             observers,
+            field="B",
             sumup=False,
             squeeze=squeeze,
             pixel_agg=pixel_agg,
             output=output,
-            field="B",
+            in_out=in_out,
         )
 
-    def getH(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
+    def getH(
+        self, *observers, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
         """Compute the H-field in units of A/m at observers generated by the source.
 
         Parameters
@@ -148,6 +166,19 @@ def getH(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
             returns a 2D-table as a `pandas.DataFrame` object (the Pandas library must be
             installed).
 
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
         Returns
         -------
         H-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
@@ -165,9 +196,9 @@ def getH(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
         >>> src = magpy.magnet.Sphere(polarization=(0,0,1.), diameter=1)
         >>> H = src.getH(((0,0,0), (1,0,0), (2,0,0)))
         >>> print(H)
-        [[      0.               0.         -265258.23848649]
-         [      0.               0.          -33157.27981081]
-         [      0.               0.           -4144.65997635]]
+        [[      0.               0.         -265258.23834209]
+         [      0.               0.          -33157.27979276]
+         [      0.               0.           -4144.6599741 ]]
 
         Compute the H-field at two sensors, each one with two pixels
 
@@ -175,22 +206,147 @@ def getH(self, *observers, squeeze=True, pixel_agg=None, output="ndarray"):
         >>> sens2 = sens1.copy(position=(2,0,0))
         >>> H = src.getH(sens1, sens2)
         >>> print(H)
-        [[[  9702.7918453       0.         -31695.78669464]
-          [ -9702.7918453       0.         -31695.78669464]]
+        [[[  9702.79184001      0.         -31695.78667738]
+          [ -9702.79184001      0.         -31695.78667738]]
         <BLANKLINE>
-         [[   617.83031378      0.          -4098.27441472]
-          [  -617.83031378      0.          -4098.27441472]]]
+         [[   617.83031344      0.          -4098.27441249]
+          [  -617.83031344      0.          -4098.27441249]]]
 
         """
         observers = format_star_input(observers)
         return getBH_level2(
             self,
             observers,
+            field="H",
             sumup=False,
             squeeze=squeeze,
             pixel_agg=pixel_agg,
             output=output,
-            field="H",
+            in_out=in_out,
+        )
+
+    def getM(
+        self, *observers, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
+        """Compute the M-field in units of A/m at observers generated by the source.
+
+        Parameters
+        ----------
+        observers: array_like or (list of) `Sensor` objects
+            Can be array_like positions of shape (n1, n2, ..., 3) where the field
+            should be evaluated, a `Sensor` object with pixel shape (n1, n2, ..., 3) or a list
+            of such sensor objects (must all have similar pixel shapes). All positions
+            are given in units of m.
+
+        squeeze: bool, default=`True`
+            If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
+            only a single source) are eliminated.
+
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observers input with different (pixel) shapes is allowed.
+
+        output: str, default='ndarray'
+            Output type, which must be one of `('ndarray', 'dataframe')`. By default a multi-
+            dimensional array ('ndarray') is returned. If 'dataframe' is chosen, the function
+            returns a 2D-table as a `pandas.DataFrame` object (the Pandas library must be
+            installed).
+
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
+        Returns
+        -------
+        M-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
+            M-field at each path position (index m) for each sensor (index k) and each sensor
+            pixel position (indices n1,n2,...) in units of A/m. Sensor pixel positions are
+            equivalent to simple observer positions. Paths of objects that are shorter than
+            index m will be considered as static beyond their end.
+        """
+        observers = format_star_input(observers)
+        return getBH_level2(
+            self,
+            observers,
+            field="M",
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            output=output,
+            in_out=in_out,
+        )
+
+    def getJ(
+        self, *observers, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
+        """Compute the J-field at observers in units of T generated by the source.
+
+        SI units are used for all inputs and outputs.
+
+        Parameters
+        ----------
+        observers: array_like or (list of) `Sensor` objects
+            Can be array_like positions of shape (n1, n2, ..., 3) where the field
+            should be evaluated, a `Sensor` object with pixel shape (n1, n2, ..., 3) or a list
+            of such sensor objects (must all have similar pixel shapes). All positions are given
+            in units of m.
+
+        squeeze: bool, default=`True`
+            If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
+            only a single source) are eliminated.
+
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observers input with different (pixel) shapes is allowed.
+
+        output: str, default='ndarray'
+            Output type, which must be one of `('ndarray', 'dataframe')`. By default a multi-
+            dimensional array ('ndarray') is returned. If 'dataframe' is chosen, the function
+            returns a 2D-table as a `pandas.DataFrame` object (the Pandas library must be
+            installed).
+
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
+        Returns
+        -------
+        J-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
+            J-field at each path position (index m) for each sensor (index k) and each sensor
+            pixel position (indices n1,n2,...) in units of T. Sensor pixel positions are equivalent
+            to simple observer positions. Paths of objects that are shorter than index m will be
+            considered as static beyond their end.
+        """
+        observers = format_star_input(observers)
+        return getBH_level2(
+            self,
+            observers,
+            field="J",
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            output=output,
+            in_out=in_out,
         )
 
 
diff --git a/magpylib/_src/obj_classes/class_Collection.py b/magpylib/_src/obj_classes/class_Collection.py
index f0a484f23..51a93a4c5 100644
--- a/magpylib/_src/obj_classes/class_Collection.py
+++ b/magpylib/_src/obj_classes/class_Collection.py
@@ -518,7 +518,9 @@ def _validate_getBH_inputs(self, *inputs):
             sources, sensors = self, inputs
         return sources, sensors
 
-    def getB(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
+    def getB(
+        self, *inputs, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
         """Compute B-field for given sources and observers.
 
         SI units are used for all inputs and outputs.
@@ -545,6 +547,19 @@ def getB(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
             `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
             object is returned (the Pandas library must be installed).
 
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
         Returns
         -------
         B-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
@@ -584,14 +599,17 @@ def getB(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
         return getBH_level2(
             sources,
             sensors,
+            field="B",
             sumup=False,
             squeeze=squeeze,
             pixel_agg=pixel_agg,
             output=output,
-            field="B",
+            in_out=in_out,
         )
 
-    def getH(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
+    def getH(
+        self, *inputs, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
         """Compute H-field for given sources and observers.
 
         SI units are used for all inputs and outputs.
@@ -626,6 +644,19 @@ def getH(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
             equivalent to simple observer positions. Paths of objects that are shorter than
             index m are considered as static beyond their end.
 
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
         Examples
         --------
         In this example we create a collection from two sources and two sensors:
@@ -645,11 +676,11 @@ def getH(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
         >>> H = magpy.getH([src1, src2], col)
         >>> H = magpy.getH([src1, src2], [sens1, sens2])
         >>> print(H)
-        [[[    0.             0.         66314.55962162]
-          [    0.             0.         66314.55962162]]
+        [[[    0.             0.         66314.55958552]
+          [    0.             0.         66314.55958552]]
         <BLANKLINE>
-         [[    0.             0.         66314.55962162]
-          [    0.             0.         66314.55962162]]]
+         [[    0.             0.         66314.55958552]
+          [    0.             0.         66314.55958552]]]
         """
 
         sources, sensors = self._validate_getBH_inputs(*inputs)
@@ -657,11 +688,140 @@ def getH(self, *inputs, squeeze=True, pixel_agg=None, output="ndarray"):
         return getBH_level2(
             sources,
             sensors,
+            field="H",
             sumup=False,
             squeeze=squeeze,
             pixel_agg=pixel_agg,
             output=output,
-            field="H",
+            in_out=in_out,
+        )
+
+    def getM(
+        self, *inputs, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
+        """Compute M-field for given sources and observers.
+
+        SI units are used for all inputs and outputs.
+
+        Parameters
+        ----------
+        inputs: source or observer objects
+            Input can only be observers if the collection contains only sources. In this case the
+            collection behaves like a single source.
+            Input can only be sources if the collection contains sensors. In this case the
+            collection behaves like a list of all its sensors.
+
+        squeeze: bool, default=`True`
+            If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
+            only a single source) are eliminated.
+
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observers input with different (pixel) shapes is allowed.
+
+        output: str, default='ndarray'
+            Output type, which must be one of `('ndarray', 'dataframe')`. By default a
+            `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
+            object is returned (the Pandas library must be installed).
+
+        Returns
+        -------
+        M-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
+            M-field at each path position (index m) for each sensor (index k) and each sensor
+            pixel position (indeices n1,n2,...) in units of A/m. Sensor pixel positions are
+            equivalent to simple observer positions. Paths of objects that are shorter than
+            index m are considered as static beyond their end.
+
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+        """
+
+        sources, sensors = self._validate_getBH_inputs(*inputs)
+
+        return getBH_level2(
+            sources,
+            sensors,
+            field="M",
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            output=output,
+            in_out=in_out,
+        )
+
+    def getJ(
+        self, *inputs, squeeze=True, pixel_agg=None, output="ndarray", in_out="auto"
+    ):
+        """Compute J-field for given sources and observers.
+
+        SI units are used for all inputs and outputs.
+
+        Parameters
+        ----------
+        inputs: source or observer objects
+            Input can only be observers if the collection contains only sources. In this case the
+            collection behaves like a single source.
+            Input can only be sources if the collection contains only sensors. In this case the
+            collection behaves like a list of all its sensors.
+
+        squeeze: bool, default=`True`
+            If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
+            only a single source) are eliminated.
+
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observers input with different (pixel) shapes is allowed.
+
+        output: str, default='ndarray'
+            Output type, which must be one of `('ndarray', 'dataframe')`. By default a
+            `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
+            object is returned (the Pandas library must be installed).
+
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
+        Returns
+        -------
+        J-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3) or DataFrame
+            J-field at each path position ( index m) for each sensor (index k) and each
+            sensor pixel position (indices n1,n2,...) in units of T. Sensor pixel positions
+            are equivalent to simple observer positions. Paths of objects that are shorter
+            than index m are considered as static beyond their end.
+        """
+
+        sources, sensors = self._validate_getBH_inputs(*inputs)
+
+        return getBH_level2(
+            sources,
+            sensors,
+            field="J",
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            output=output,
+            in_out=in_out,
         )
 
     @property
diff --git a/magpylib/_src/obj_classes/class_Sensor.py b/magpylib/_src/obj_classes/class_Sensor.py
index 2a4a70c3a..628729547 100644
--- a/magpylib/_src/obj_classes/class_Sensor.py
+++ b/magpylib/_src/obj_classes/class_Sensor.py
@@ -81,7 +81,7 @@ class Sensor(BaseGeo, BaseDisplayRepr):
     >>> print(B)
     [[0.00000000e+00 8.88576588e-05 8.88576588e-05]
      [0.00000000e+00 9.16274003e-05 9.16274003e-05]
-     [0.00000000e+00 1.01415383e-04 1.01415383e-04]]
+     [0.00000000e+00 1.01415384e-04 1.01415384e-04]]
     """
 
     _style_class = SensorStyle
@@ -142,7 +142,13 @@ def handedness(self, val):
         self._handedness = val
 
     def getB(
-        self, *sources, sumup=False, squeeze=True, pixel_agg=None, output="ndarray"
+        self,
+        *sources,
+        sumup=False,
+        squeeze=True,
+        pixel_agg=None,
+        output="ndarray",
+        in_out="auto",
     ):
         """Compute the B-field in units of T as seen by the sensor.
 
@@ -169,6 +175,19 @@ def getB(
             `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
             object is returned (the Pandas library must be installed).
 
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
         Returns
         -------
         B-field: ndarray, shape squeeze(l, m, n1, n2, ..., 3) or DataFrame
@@ -203,21 +222,28 @@ def getB(
         >>> print(B)
         [[0.00000000e+00 8.88576588e-05 8.88576588e-05]
          [0.00000000e+00 9.16274003e-05 9.16274003e-05]
-         [0.00000000e+00 1.01415383e-04 1.01415383e-04]]
+         [0.00000000e+00 1.01415384e-04 1.01415384e-04]]
         """
         sources = format_star_input(sources)
         return getBH_level2(
             sources,
             self,
+            field="B",
             sumup=sumup,
             squeeze=squeeze,
             pixel_agg=pixel_agg,
             output=output,
-            field="B",
+            in_out=in_out,
         )
 
     def getH(
-        self, *sources, sumup=False, squeeze=True, pixel_agg=None, output="ndarray"
+        self,
+        *sources,
+        sumup=False,
+        squeeze=True,
+        pixel_agg=None,
+        output="ndarray",
+        in_out="auto",
     ):
         """Compute the H-field in units of A/m as seen by the sensor.
 
@@ -244,6 +270,19 @@ def getH(
             `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
             object is returned (the Pandas library must be installed).
 
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
         Returns
         -------
         H-field: ndarray, shape squeeze(l, m, n1, n2, ..., 3) or DataFrame
@@ -284,11 +323,144 @@ def getH(
         return getBH_level2(
             sources,
             self,
+            field="H",
             sumup=sumup,
             squeeze=squeeze,
             pixel_agg=pixel_agg,
             output=output,
-            field="H",
+            in_out=in_out,
+        )
+
+    def getM(
+        self,
+        *sources,
+        sumup=False,
+        squeeze=True,
+        pixel_agg=None,
+        output="ndarray",
+        in_out="auto",
+    ):
+        """Compute the M-field in units of A/m as seen by the sensor.
+
+        Parameters
+        ----------
+        sources: source and collection objects or 1D list thereof
+            Sources that generate the magnetic field. Can be a single source (or collection)
+            or a 1D list of l source and/or collection objects.
+
+        sumup: bool, default=`False`
+            If `True`, the fields of all sources are summed up.
+
+        squeeze: bool, default=`True`
+            If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
+            a single sensor or only a single source) are eliminated.
+
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observers input with different (pixel) shapes is allowed.
+
+        output: str, default='ndarray'
+            Output type, which must be one of `('ndarray', 'dataframe')`. By default a
+            `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
+            object is returned (the Pandas library must be installed).
+
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
+        Returns
+        -------
+        M-field: ndarray, shape squeeze(l, m, n1, n2, ..., 3) or DataFrame
+            M-field of each source (index l) at each path position (index m) and each sensor pixel
+            position (indices n1,n2,...) in units of A/m. Paths of objects that are shorter than
+            index m are considered as static beyond their end.
+        """
+        sources = format_star_input(sources)
+        return getBH_level2(
+            sources,
+            self,
+            field="M",
+            sumup=sumup,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            output=output,
+            in_out=in_out,
+        )
+
+    def getJ(
+        self,
+        *sources,
+        sumup=False,
+        squeeze=True,
+        pixel_agg=None,
+        output="ndarray",
+        in_out="auto",
+    ):
+        """Compute the J-field in units of T as seen by the sensor.
+
+        Parameters
+        ----------
+        sources: source and collection objects or 1D list thereof
+            Sources that generate the magnetic field. Can be a single source (or collection)
+            or a 1D list of l source and/or collection objects.
+
+        sumup: bool, default=`False`
+            If `True`, the fields of all sources are summed up.
+
+        squeeze: bool, default=`True`
+            If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
+            a single sensor or only a single source) are eliminated.
+
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observers input with different (pixel) shapes is allowed.
+
+        output: str, default='ndarray'
+            Output type, which must be one of `('ndarray', 'dataframe')`. By default a
+            `numpy.ndarray` object is returned. If 'dataframe' is chosen, a `pandas.DataFrame`
+            object is returned (the Pandas library must be installed).
+
+        in_out: {'auto', 'inside', 'outside'}
+            This parameter only applies for magnet bodies. It specifies the location of the
+            observers relative to the magnet body, affecting the calculation of the magnetic field.
+            The options are:
+            - 'auto': The location (inside or outside the cuboid) is determined automatically for
+            each observer.
+            - 'inside': All observers are considered to be inside the cuboid; use this for
+              performance optimization if applicable.
+            - 'outside': All observers are considered to be outside the cuboid; use this for
+              performance optimization if applicable.
+            Choosing 'auto' is fail-safe but may be computationally intensive if the mix of observer
+            locations is unknown.
+
+        Returns
+        -------
+        J-field: ndarray, shape squeeze(l, m, n1, n2, ..., 3) or DataFrame
+            J-field of each source (index l) at each path position (index m) and each sensor pixel
+            position (indices n1,n2,...) in units of T. Paths of objects that are shorter than
+            index m are considered as static beyond their end.
+        """
+        sources = format_star_input(sources)
+        return getBH_level2(
+            sources,
+            self,
+            field="J",
+            sumup=sumup,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            output=output,
+            in_out=in_out,
         )
 
     @property
diff --git a/magpylib/_src/obj_classes/class_current_Circle.py b/magpylib/_src/obj_classes/class_current_Circle.py
index a4830722f..c83446488 100644
--- a/magpylib/_src/obj_classes/class_current_Circle.py
+++ b/magpylib/_src/obj_classes/class_current_Circle.py
@@ -3,7 +3,7 @@
 
 from magpylib._src.display.traces_core import make_Circle
 from magpylib._src.exceptions import MagpylibDeprecationWarning
-from magpylib._src.fields.field_BH_circle import current_circle_field
+from magpylib._src.fields.field_BH_circle import BHJM_circle
 from magpylib._src.input_checks import check_format_input_scalar
 from magpylib._src.obj_classes.class_BaseExcitations import BaseCurrent
 from magpylib._src.utility import unit_prefix
@@ -68,7 +68,7 @@ class Circle(BaseCurrent):
     >>> print(B)
     [[-1.63585841e-24 -4.44388287e-05  4.44388287e-05]
      [-6.55449367e-24 -4.44688604e-05  4.44688604e-05]
-     [-9.85948764e-24 -4.45190261e-05  4.45190261e-05]]
+     [-9.85948765e-24 -4.45190261e-05  4.45190261e-05]]
 
     The same result is obtained when the rotated source moves along a path away from an
     observer at position (1,1,1). This time we use a `Sensor` object as observer.
@@ -80,10 +80,10 @@ class Circle(BaseCurrent):
     >>> print(B)
     [[-1.63585841e-24 -4.44388287e-05  4.44388287e-05]
      [-6.55449367e-24 -4.44688604e-05  4.44688604e-05]
-     [-9.85948764e-24 -4.45190261e-05  4.45190261e-05]]
+     [-9.85948765e-24 -4.45190261e-05  4.45190261e-05]]
     """
 
-    _field_func = staticmethod(current_circle_field)
+    _field_func = staticmethod(BHJM_circle)
     _field_func_kwargs_ndim = {"current": 1, "diameter": 1}
     get_trace = make_Circle
 
@@ -136,7 +136,7 @@ class Loop(Circle):
     def _field_func(*args, **kwargs):
         """Catch Deprecation warning in getBH_dict"""
         _deprecation_warn()
-        return current_circle_field(*args, **kwargs)
+        return BHJM_circle(*args, **kwargs)
 
     def __init__(self, *args, **kwargs):
         _deprecation_warn()
diff --git a/magpylib/_src/obj_classes/class_current_Polyline.py b/magpylib/_src/obj_classes/class_current_Polyline.py
index bb3b15c93..2d566df86 100644
--- a/magpylib/_src/obj_classes/class_current_Polyline.py
+++ b/magpylib/_src/obj_classes/class_current_Polyline.py
@@ -71,7 +71,7 @@ class Polyline(BaseCurrent):
     >>> B = src.getB([(.01,.01,.01), (.02,.02,.02), (.03,.03,.03)])
     >>> print(B)
     [[-1.04529728e-21  3.50341393e-06 -3.50341393e-06]
-     [-9.28140348e-23  3.62181325e-07 -3.62181325e-07]
+     [-9.28140349e-23  3.62181325e-07 -3.62181325e-07]
      [-1.72744075e-23  1.03643004e-07 -1.03643004e-07]]
 
     The same result is obtained when the rotated source moves along a path away from an
@@ -83,7 +83,7 @@ class Polyline(BaseCurrent):
     >>> B = src.getB(sens)
     >>> print(B)
     [[-1.04529728e-21  3.50341393e-06 -3.50341393e-06]
-     [-9.28140348e-23  3.62181325e-07 -3.62181325e-07]
+     [-9.28140349e-23  3.62181325e-07 -3.62181325e-07]
      [-1.72744075e-23  1.03643004e-07 -1.03643004e-07]]
     """
 
diff --git a/magpylib/_src/obj_classes/class_magnet_Cuboid.py b/magpylib/_src/obj_classes/class_magnet_Cuboid.py
index 2e3e4646d..dd32a7d45 100644
--- a/magpylib/_src/obj_classes/class_magnet_Cuboid.py
+++ b/magpylib/_src/obj_classes/class_magnet_Cuboid.py
@@ -1,6 +1,6 @@
 """Magnet Cuboid class code"""
 from magpylib._src.display.traces_core import make_Cuboid
-from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
+from magpylib._src.fields.field_BH_cuboid import BHJM_magnet_cuboid
 from magpylib._src.input_checks import check_format_input_vector
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
 from magpylib._src.utility import unit_prefix
@@ -61,7 +61,7 @@ class Cuboid(BaseMagnet):
     >>> src = magpy.magnet.Cuboid(polarization=(.5,.6,.7), dimension=(.01,.01,.01))
     >>> H = src.getH((.01,.01,.01))
     >>> print(H)
-    [16149.04136518 14906.80741401 13664.57346284]
+    [16149.04135639 14906.8074059  13664.57345541]
 
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
@@ -74,7 +74,7 @@ class Cuboid(BaseMagnet):
      [-0.03557183  0.00646436  0.14943466]]
     """
 
-    _field_func = staticmethod(magnet_cuboid_field)
+    _field_func = staticmethod(BHJM_magnet_cuboid)
     _field_func_kwargs_ndim = {"polarization": 2, "dimension": 2}
     get_trace = make_Cuboid
 
diff --git a/magpylib/_src/obj_classes/class_magnet_Cylinder.py b/magpylib/_src/obj_classes/class_magnet_Cylinder.py
index 61b289334..c876cc7e7 100644
--- a/magpylib/_src/obj_classes/class_magnet_Cylinder.py
+++ b/magpylib/_src/obj_classes/class_magnet_Cylinder.py
@@ -1,6 +1,6 @@
 """Magnet Cylinder class code"""
 from magpylib._src.display.traces_core import make_Cylinder
-from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_field
+from magpylib._src.fields.field_BH_cylinder import BHJM_magnet_cylinder
 from magpylib._src.input_checks import check_format_input_vector
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
 from magpylib._src.utility import unit_prefix
@@ -60,7 +60,7 @@ class Cylinder(BaseMagnet):
     >>> src = magpy.magnet.Cylinder(polarization=(.1,.2,.3), dimension=(.01,.01))
     >>> H = src.getH((.01,.01,.01))
     >>> print(H)
-    [4849.91343385 3883.17815728 2739.73202386]
+    [4849.91343121 3883.17815517 2739.73202237]
 
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
@@ -85,7 +85,7 @@ class Cylinder(BaseMagnet):
      [1.25715233e-04 2.01445027e-04 1.31238931e-05]]
     """
 
-    _field_func = staticmethod(magnet_cylinder_field)
+    _field_func = staticmethod(BHJM_magnet_cylinder)
     _field_func_kwargs_ndim = {"polarization": 2, "dimension": 2}
     get_trace = make_Cylinder
 
diff --git a/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py b/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py
index e5445c6a4..e3c1ae3b9 100644
--- a/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py
+++ b/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py
@@ -3,7 +3,7 @@
 
 from magpylib._src.display.traces_core import make_CylinderSegment
 from magpylib._src.fields.field_BH_cylinder_segment import (
-    magnet_cylinder_segment_field_internal,
+    BHJM_cylinder_segment_internal,
 )
 from magpylib._src.input_checks import check_format_input_cylinder_segment
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
@@ -75,7 +75,7 @@ class CylinderSegment(BaseMagnet):
     >>> src = magpy.magnet.CylinderSegment(polarization=(.1,.2,.3), dimension=(.01,.02,.01,0,45))
     >>> H = src.getH((.02,.02,.02))
     >>> print(H)
-    [ 807.84692335 1934.22812967 2741.16804712]
+    [ 807.84692247 1934.22812757 2741.16804414]
 
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
@@ -100,7 +100,7 @@ class CylinderSegment(BaseMagnet):
      [ 0.00025439  0.00074332  0.00011683]]
     """
 
-    _field_func = staticmethod(magnet_cylinder_segment_field_internal)
+    _field_func = staticmethod(BHJM_cylinder_segment_internal)
     _field_func_kwargs_ndim = {"polarization": 2, "dimension": 2}
     get_trace = make_CylinderSegment
 
diff --git a/magpylib/_src/obj_classes/class_magnet_Sphere.py b/magpylib/_src/obj_classes/class_magnet_Sphere.py
index b002f0c0f..6f4f4b615 100644
--- a/magpylib/_src/obj_classes/class_magnet_Sphere.py
+++ b/magpylib/_src/obj_classes/class_magnet_Sphere.py
@@ -1,6 +1,6 @@
 """Magnet Sphere class code"""
 from magpylib._src.display.traces_core import make_Sphere
-from magpylib._src.fields.field_BH_sphere import magnet_sphere_field
+from magpylib._src.fields.field_BH_sphere import BHJM_magnet_sphere
 from magpylib._src.input_checks import check_format_input_scalar
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
 from magpylib._src.utility import unit_prefix
@@ -60,7 +60,7 @@ class Sphere(BaseMagnet):
     >>> src = magpy.magnet.Sphere(polarization=(.1,.2,.3), diameter=.01)
     >>> H = src.getH((.01,.01,.01))
     >>> print(H)
-    [3190.56073739 2552.44858992 1914.33644244]
+    [3190.56073566 2552.44858853 1914.33644139]
 
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
@@ -85,7 +85,7 @@ class Sphere(BaseMagnet):
      [8.40017059e-05 1.34701220e-04 8.69866146e-06]]
     """
 
-    _field_func = staticmethod(magnet_sphere_field)
+    _field_func = staticmethod(BHJM_magnet_sphere)
     _field_func_kwargs_ndim = {"polarization": 2, "diameter": 1}
     get_trace = make_Sphere
 
diff --git a/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py b/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py
index 66409c651..8ed542a6c 100644
--- a/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py
+++ b/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 from magpylib._src.display.traces_core import make_Tetrahedron
-from magpylib._src.fields.field_BH_tetrahedron import magnet_tetrahedron_field
+from magpylib._src.fields.field_BH_tetrahedron import BHJM_magnet_tetrahedron
 from magpylib._src.input_checks import check_format_input_vector
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
 
@@ -72,7 +72,7 @@ class Tetrahedron(BaseMagnet):
     >>> src = magpy.magnet.Tetrahedron(polarization=(.1,.2,.3), vertices=verts)
     >>> H = src.getH((.01,.01,.01))
     >>> print(H)
-    [2070.89782733 1656.71826186 1242.5386964 ]
+    [2070.8978262  1656.71826096 1242.53869572]
 
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
@@ -97,7 +97,7 @@ class Tetrahedron(BaseMagnet):
      [ 2.90426931e-05  5.22556994e-05 -1.70596096e-06]]
     """
 
-    _field_func = staticmethod(magnet_tetrahedron_field)
+    _field_func = staticmethod(BHJM_magnet_tetrahedron)
     _field_func_kwargs_ndim = {"polarization": 1, "vertices": 3}
     get_trace = make_Tetrahedron
 
diff --git a/magpylib/_src/obj_classes/class_magnet_TriangularMesh.py b/magpylib/_src/obj_classes/class_magnet_TriangularMesh.py
index 5e987b5f5..eacc13676 100644
--- a/magpylib/_src/obj_classes/class_magnet_TriangularMesh.py
+++ b/magpylib/_src/obj_classes/class_magnet_TriangularMesh.py
@@ -6,6 +6,7 @@
 
 from magpylib._src.display.traces_core import make_TriangularMesh
 from magpylib._src.exceptions import MagpylibMissingInput
+from magpylib._src.fields.field_BH_triangularmesh import BHJM_magnet_trimesh
 from magpylib._src.fields.field_BH_triangularmesh import calculate_centroid
 from magpylib._src.fields.field_BH_triangularmesh import fix_trimesh_orientation
 from magpylib._src.fields.field_BH_triangularmesh import (
@@ -13,7 +14,6 @@
 )
 from magpylib._src.fields.field_BH_triangularmesh import get_intersecting_triangles
 from magpylib._src.fields.field_BH_triangularmesh import get_open_edges
-from magpylib._src.fields.field_BH_triangularmesh import magnet_trimesh_field
 from magpylib._src.input_checks import check_format_input_vector
 from magpylib._src.input_checks import check_format_input_vector2
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
@@ -114,7 +114,7 @@ class TriangularMesh(BaseMagnet):
     [0.00260237 0.00208189 0.00156142]
     """
 
-    _field_func = staticmethod(magnet_trimesh_field)
+    _field_func = staticmethod(BHJM_magnet_trimesh)
     _field_func_kwargs_ndim = {"polarization": 2, "mesh": 3}
     get_trace = make_TriangularMesh
     _style_class = TriangularMeshStyle
diff --git a/magpylib/_src/obj_classes/class_misc_Dipole.py b/magpylib/_src/obj_classes/class_misc_Dipole.py
index ae18537df..192cc8928 100644
--- a/magpylib/_src/obj_classes/class_misc_Dipole.py
+++ b/magpylib/_src/obj_classes/class_misc_Dipole.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 from magpylib._src.display.traces_core import make_Dipole
-from magpylib._src.fields.field_BH_dipole import dipole_field
+from magpylib._src.fields.field_BH_dipole import BHJM_dipole
 from magpylib._src.input_checks import check_format_input_vector
 from magpylib._src.obj_classes.class_BaseExcitations import BaseSource
 from magpylib._src.style import DipoleStyle
@@ -81,7 +81,7 @@ class Dipole(BaseSource):
      [0.0100802  0.01720799 0.00712778]]
     """
 
-    _field_func = staticmethod(dipole_field)
+    _field_func = staticmethod(BHJM_dipole)
     _field_func_kwargs_ndim = {"moment": 2}
     _style_class = DipoleStyle
     get_trace = make_Dipole
diff --git a/magpylib/_src/obj_classes/class_misc_Triangle.py b/magpylib/_src/obj_classes/class_misc_Triangle.py
index 8f14f974c..8e4da601f 100644
--- a/magpylib/_src/obj_classes/class_misc_Triangle.py
+++ b/magpylib/_src/obj_classes/class_misc_Triangle.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 from magpylib._src.display.traces_core import make_Triangle
-from magpylib._src.fields.field_BH_triangle import triangle_field
+from magpylib._src.fields.field_BH_triangle import BHJM_triangle
 from magpylib._src.input_checks import check_format_input_vector
 from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
 from magpylib._src.style import TriangleStyle
@@ -75,7 +75,7 @@ class Triangle(BaseMagnet):
     >>> src = magpy.misc.Triangle(polarization=(.1,.2,.3), vertices=verts)
     >>> H = src.getH((.1,.1,.1))
     >>> print(H)
-    [18.88869831 18.88869831 19.54560637]
+    [18.8886983  18.8886983  19.54560636]
 
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
@@ -100,7 +100,7 @@ class Triangle(BaseMagnet):
      [0.00030049 0.00031044 0.00031044]]
     """
 
-    _field_func = staticmethod(triangle_field)
+    _field_func = staticmethod(BHJM_triangle)
     _field_func_kwargs_ndim = {"polarization": 2, "vertices": 2}
     get_trace = make_Triangle
     _style_class = TriangleStyle
diff --git a/magpylib/_src/utility.py b/magpylib/_src/utility.py
index f4eef2514..3769b33eb 100644
--- a/magpylib/_src/utility.py
+++ b/magpylib/_src/utility.py
@@ -2,7 +2,10 @@
 # pylint: disable=import-outside-toplevel
 # pylint: disable=cyclic-import
 # import numbers
+from functools import lru_cache
+from inspect import signature
 from math import log10
+from typing import Callable
 from typing import Sequence
 
 import numpy as np
@@ -21,8 +24,6 @@ def get_allowed_sources_msg():
 - string {srcs}"""
 
 
-MU0 = 4 * np.pi * 1e-7
-
 ALLOWED_OBSERVER_MSG = """Observers must be either
 - array_like positions of shape (N1, N2, ..., 3)
 - Sensor object
@@ -389,3 +390,10 @@ def open_animation(filepath, embed=True):
         import webbrowser
 
         webbrowser.open(filepath)
+
+
+@lru_cache(maxsize=None)
+def has_parameter(func: Callable, param_name: str) -> bool:
+    """Check if input function has a specific parameter"""
+    sig = signature(func)
+    return param_name in sig.parameters
diff --git a/magpylib/core/__init__.py b/magpylib/core/__init__.py
index 03aed10d3..fc284e7a6 100644
--- a/magpylib/core/__init__.py
+++ b/magpylib/core/__init__.py
@@ -3,27 +3,25 @@
 """
 
 __all__ = [
-    "dipole_field",
-    "current_circle_field",
-    "current_loop_field",
-    "current_line_field",
-    "current_polyline_field",
-    "magnet_sphere_field",
-    "magnet_cuboid_field",
-    "magnet_cylinder_field",
-    "magnet_cylinder_segment_field",
-    "triangle_field",
-    "magnet_tetrahedron_field",
+    "magnet_cuboid_Bfield",
+    "magnet_cylinder_axial_Bfield",
+    "magnet_cylinder_diametral_Hfield",
+    "magnet_cylinder_segment_Hfield",
+    "magnet_sphere_Bfield",
+    "current_circle_Hfield",
+    "current_polyline_Hfield",
+    "dipole_Hfield",
+    "triangle_Bfield",
 ]
 
-from magpylib._src.fields.field_BH_dipole import dipole_field
-from magpylib._src.fields.field_BH_circle import current_circle_field
-from magpylib._src.fields.field_BH_circle import current_loop_field
-from magpylib._src.fields.field_BH_polyline import current_line_field
-from magpylib._src.fields.field_BH_polyline import current_polyline_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
-from magpylib._src.fields.field_BH_triangle import triangle_field
-from magpylib._src.fields.field_BH_tetrahedron import magnet_tetrahedron_field
+from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_Bfield
+from magpylib._src.fields.field_BH_cylinder import magnet_cylinder_axial_Bfield
+from magpylib._src.fields.field_BH_cylinder import magnet_cylinder_diametral_Hfield
+from magpylib._src.fields.field_BH_cylinder_segment import (
+    magnet_cylinder_segment_Hfield,
+)
+from magpylib._src.fields.field_BH_sphere import magnet_sphere_Bfield
+from magpylib._src.fields.field_BH_circle import current_circle_Hfield
+from magpylib._src.fields.field_BH_polyline import current_polyline_Hfield
+from magpylib._src.fields.field_BH_dipole import dipole_Hfield
+from magpylib._src.fields.field_BH_triangle import triangle_Bfield


Test patch:
diff --git a/tests/test_core_misc.py b/tests/test_BHMJ_level.py
similarity index 76%
rename from tests/test_core_misc.py
rename to tests/test_BHMJ_level.py
index 3a9c7526a..b23555dd2 100644
--- a/tests/test_core_misc.py
+++ b/tests/test_BHMJ_level.py
@@ -1,21 +1,18 @@
 import numpy as np
-import pytest
 from numpy.testing import assert_allclose
-
-import magpylib as magpy
-from magpylib._src.exceptions import MagpylibDeprecationWarning
+from scipy.constants import mu_0 as MU0
+
+from magpylib._src.fields.field_BH_circle import BHJM_circle
+from magpylib._src.fields.field_BH_cuboid import BHJM_magnet_cuboid
+from magpylib._src.fields.field_BH_cylinder import BHJM_magnet_cylinder
+from magpylib._src.fields.field_BH_cylinder_segment import BHJM_cylinder_segment
+from magpylib._src.fields.field_BH_dipole import BHJM_dipole
+from magpylib._src.fields.field_BH_polyline import BHJM_current_polyline
 from magpylib._src.fields.field_BH_polyline import current_vertices_field
-from magpylib._src.fields.field_BH_triangularmesh import magnet_trimesh_field
-from magpylib._src.utility import MU0
-from magpylib.core import current_circle_field
-from magpylib.core import current_line_field
-from magpylib.core import current_loop_field
-from magpylib.core import current_polyline_field
-from magpylib.core import magnet_cuboid_field
-from magpylib.core import magnet_cylinder_segment_field
-from magpylib.core import magnet_sphere_field
-from magpylib.core import magnet_tetrahedron_field
-from magpylib.core import triangle_field
+from magpylib._src.fields.field_BH_sphere import BHJM_magnet_sphere
+from magpylib._src.fields.field_BH_tetrahedron import BHJM_magnet_tetrahedron
+from magpylib._src.fields.field_BH_triangle import BHJM_triangle
+from magpylib._src.fields.field_BH_triangularmesh import BHJM_magnet_trimesh
 
 
 #######################################################################################
@@ -25,7 +22,21 @@
 # NEW V5 BASIC FIELD COMPUTATION TESTS
 
 
-def test_magnet_cuboid_field_BH():
+def helper_check_HBMJ_consistency(func, **kw):
+    """
+    helper function to check H,B,M,J field consistencies
+    returns H, B, M, J
+    """
+    B = func(field="B", **kw)
+    H = func(field="H", **kw)
+    M = func(field="M", **kw)
+    J = func(field="J", **kw)
+    np.testing.assert_allclose(M * MU0, J)
+    np.testing.assert_allclose(B, MU0 * H + J)
+    return H, B, M, J
+
+
+def test_BHJM_magnet_cuboid():
     """test cuboid field"""
     pol = np.array(
         [
@@ -54,24 +65,15 @@ def test_magnet_cuboid_field_BH():
             (1, -1, 0),
             (1, -1, 0),
             (1, 2, 3),
-            (0, 0, 0),
+            (0, 0, 0),  # inside
         ]
     )
-
-    B = magnet_cuboid_field(
-        field="B",
-        observers=obs,
-        polarization=pol,
-        dimension=dim,
-    )
-    H = magnet_cuboid_field(
-        field="H",
-        observers=obs,
-        polarization=pol,
-        dimension=dim,
-    )
-    J = np.array([(0, 0, 0)] * 5 + [(1, 2, 3)])
-    np.testing.assert_allclose(B, MU0 * H + J)
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "dimension": dim,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_magnet_cuboid, **kw)
 
     Btest = [
         [0.0, 0.0, 0.0],
@@ -93,8 +95,15 @@ def test_magnet_cuboid_field_BH():
     ]
     np.testing.assert_allclose(H, Htest, rtol=1e-5)
 
+    Jtest = np.array([(0, 0, 0)] * 5 + [(1, 2, 3)])
+    np.testing.assert_allclose(J, Jtest, rtol=1e-5)
 
-def test_magnet_cylinder_field_BH():
+    # H_inout = BHJM_magnet_cuboid(field="H", in_out="outside", **kw)
+    # Htest_inout = Htest + Jtest / MU0
+    # np.testing.assert_allclose(H_inout, Htest_inout, rtol=1e-5)
+
+
+def test_BHJM_magnet_cylinder():
     """test cylinder field computation"""
     pol = np.array(
         [
@@ -117,23 +126,16 @@ def test_magnet_cylinder_field_BH():
             (1, 2, 3),
             (1, -1, 0),
             (1, 1, 1),
-            (0, 0, 0),
+            (0, 0, 0),  # inside
         ]
     )
-    B = magpy.core.magnet_cylinder_field(
-        field="B",
-        observers=obs,
-        polarization=pol,
-        dimension=dim,
-    )
-    H = magpy.core.magnet_cylinder_field(
-        field="H",
-        observers=obs,
-        polarization=pol,
-        dimension=dim,
-    )
-    J = np.array([(0, 0, 0), (0, 0, 0), (0, 0, 0), (1, 1, 1)])
-    np.testing.assert_allclose(B, MU0 * H + J)
+
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "dimension": dim,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_magnet_cylinder, **kw)
 
     Btest = [
         [0.0, 0.0, 0.0],
@@ -151,9 +153,16 @@ def test_magnet_cylinder_field_BH():
     ]
     np.testing.assert_allclose(H, Htest)
 
+    Jtest = np.array([(0, 0, 0), (0, 0, 0), (0, 0, 0), (1, 1, 1)])
+    np.testing.assert_allclose(J, Jtest)
+
+    # H_inout = BHJM_magnet_cylinder(field="H", in_out="outside", **kw)
+    # Htest_inout = Htest - Jtest / MU0
+    # np.testing.assert_allclose(H_inout, Htest_inout, rtol=1e-5)
+
 
-def test_magnet_sphere_field_BH():
-    """test magnet_sphere_field"""
+def test_BHJM_magnet_sphere():
+    """test BHJM_magnet_sphere"""
     pol = np.array(
         [
             (0, 0, 0),
@@ -167,24 +176,17 @@ def test_magnet_sphere_field_BH():
         [
             (1, 2, 3),
             (1, -1, 0),
-            (0, -1, 0),
-            (1, -1, 0.5),
+            (0, -1, 0),  # inside
+            (1, -1, 0.5),  # inside
         ]
     )
-    B = magnet_sphere_field(
-        field="B",
-        observers=obs,
-        diameter=dia,
-        polarization=pol,
-    )
-    H = magnet_sphere_field(
-        field="H",
-        observers=obs,
-        diameter=dia,
-        polarization=pol,
-    )
-    J = np.array([(0, 0, 0), (0, 0, 0), pol[2], pol[3]])
-    np.testing.assert_allclose(B, MU0 * H + J)
+
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "diameter": dia,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_magnet_sphere, **kw)
 
     Btest = [
         [0.0, 0.0, 0.0],
@@ -202,6 +204,9 @@ def test_magnet_sphere_field_BH():
     ]
     np.testing.assert_allclose(H, Htest)
 
+    Jtest = [(0, 0, 0), (0, 0, 0), (2, 3, -1), (2, 3, -1)]
+    np.testing.assert_allclose(J, Jtest)
+
 
 def test_field_cylinder_segment_BH():
     """CylinderSegmetn field test"""
@@ -226,23 +231,16 @@ def test_field_cylinder_segment_BH():
             (1, 2, 3),
             (1, -1, 0),
             (0, -1, 0),
-            (1, -1, 0.5),
+            (1, -1, 0.5),  # inside
         ]
     )
-    B = magnet_cylinder_segment_field(
-        field="B",
-        observers=obs,
-        dimension=dim,
-        polarization=pol,
-    )
-    H = magnet_cylinder_segment_field(
-        field="H",
-        observers=obs,
-        dimension=dim,
-        polarization=pol,
-    )
-    J = np.array([(0, 0, 0)] * 3 + [pol[3]])
-    np.testing.assert_allclose(B, MU0 * H + J)
+
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "dimension": dim,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_cylinder_segment, **kw)
 
     Btest = [
         [0.0, 0.0, 0.0],
@@ -260,8 +258,11 @@ def test_field_cylinder_segment_BH():
     ]
     np.testing.assert_allclose(H, Htest, rtol=1e-6)
 
+    Jtest = [(0, 0, 0), (0, 0, 0), (0, 0, 0), (2, 3, -1)]
+    np.testing.assert_allclose(J, Jtest)
+
 
-def test_triangle_field_BH():
+def test_BHJM_triangle_BH():
     """Test of triangle field core function"""
     pol = np.array(
         [
@@ -287,19 +288,13 @@ def test_triangle_field_BH():
             (2, 3, 1),
         ]
     )
-    B = triangle_field(
-        field="B",
-        observers=obs,
-        vertices=vert,
-        polarization=pol,
-    )
-    H = triangle_field(
-        field="H",
-        observers=obs,
-        vertices=vert,
-        polarization=pol,
-    )
-    np.testing.assert_allclose(B, MU0 * H)
+
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "vertices": vert,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_triangle, **kw)
 
     Btest = [
         [0.0, 0.0, 0.0],
@@ -317,6 +312,9 @@ def test_triangle_field_BH():
     ]
     np.testing.assert_allclose(H, Htest, rtol=1e-06)
 
+    Jtest = [(0, 0, 0), (0, 0, 0), (0, 0, 0), (0, 0, 0)]
+    np.testing.assert_allclose(J, Jtest)
+
 
 def test_magnet_tetrahedron_field_BH():
     """Test of tetrahedron field core function"""
@@ -324,8 +322,9 @@ def test_magnet_tetrahedron_field_BH():
         [
             (0, 0, 0),
             (1, 2, 3),
-            (-1, 0.5, 0.1),
+            (-1, 0.5, 0.1),  # inside
             (2, 2, -1),
+            (3, 2, 1),  # inside
         ]
     )
     vert = np.array(
@@ -334,6 +333,7 @@ def test_magnet_tetrahedron_field_BH():
             [(0, 0, 0), (0, 1, 0), (1, 0, 0), (0, 0, 1)],
             [(-1, 0, -1), (1, 1, -1), (1, -1, -1), (0, 0, 1)],
             [(-1, 0, -1), (1, 1, -1), (1, -1, -1), (0, 0, 1)],
+            [(-10, 0, -10), (10, 10, -10), (10, -10, -10), (0, 0, 10)],
         ]
     )
     obs = np.array(
@@ -342,28 +342,22 @@ def test_magnet_tetrahedron_field_BH():
             (1, 1, 1),
             (0, 0, 0),
             (2, 0, 0),
+            (1, 2, 3),
         ]
     )
-    B = magnet_tetrahedron_field(
-        field="B",
-        observers=obs,
-        vertices=vert,
-        polarization=pol,
-    )
-    H = magnet_tetrahedron_field(
-        field="H",
-        observers=obs,
-        vertices=vert,
-        polarization=pol,
-    )
-    J = np.array([(0, 0, 0)] * 2 + [pol[2]] + [(0, 0, 0)])
-    np.testing.assert_allclose(B, MU0 * H + J)
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "vertices": vert,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_magnet_tetrahedron, **kw)
 
     Btest = [
         [0.0, 0.0, 0.0],
         [0.02602367, 0.02081894, 0.0156142],
         [-0.69704332, 0.20326329, 0.11578416],
         [0.04004769, -0.03186713, 0.03854207],
+        [2.09887014, 1.42758632, 0.8611617],
     ]
     np.testing.assert_allclose(B, Btest, rtol=1e-06)
 
@@ -372,12 +366,16 @@ def test_magnet_tetrahedron_field_BH():
         [20708.97827326, 16567.1826186, 12425.38696395],
         [241085.26350642, -236135.56979233, 12560.63814427],
         [31868.94160192, -25359.05664996, 30670.80436549],
+        [-717096.35551784, -455512.33538799, -110484.00786285],
     ]
     np.testing.assert_allclose(H, Htest, rtol=1e-06)
 
+    Jtest = [(0, 0, 0), (0, 0, 0), (-1, 0.5, 0.1), (0, 0, 0), (3, 2, 1)]
+    np.testing.assert_allclose(J, Jtest, rtol=1e-06)
 
-def test_magnet_trimesh_field_BH():
-    """Test of magnet_trimesh_field core-like function"""
+
+def test_BHJM_magnet_trimesh_BH():
+    """Test of BHJM_magnet_trimesh core-like function"""
 
     mesh1 = [
         [
@@ -401,7 +399,7 @@ def test_magnet_trimesh_field_BH():
             [-0.5703949332237244, -0.25462886691093445, 0.7809056639671326],
         ],
     ]
-    mesh2 = [
+    mesh2 = [  # inside
         [
             [0.9744000434875488, 0.15463787317276, 0.16319207847118378],
             [-0.12062954157590866, -0.8440634608268738, -0.522499144077301],
@@ -423,23 +421,15 @@ def test_magnet_trimesh_field_BH():
             [-0.47620221972465515, -0.0791524201631546, 0.8757661581039429],
         ],
     ]
-    meshes = np.array([mesh1, mesh2])
+    mesh = np.array([mesh1, mesh2])
     pol = np.array([(1, 2, 3), (3, 2, 1)])
     obs = np.array([(1, 2, 3), (0, 0, 0)])
-    B = magnet_trimesh_field(
-        field="B",
-        observers=obs,
-        mesh=meshes,
-        polarization=pol,
-    )
-    H = magnet_trimesh_field(
-        field="H",
-        observers=obs,
-        mesh=meshes,
-        polarization=pol,
-    )
-    J = np.array([(0, 0, 0), (3, 2, 1)])
-    np.testing.assert_allclose(B, MU0 * H + J)
+    kw = {
+        "observers": obs,
+        "polarization": pol,
+        "mesh": mesh,
+    }
+    H, B, _, J = helper_check_HBMJ_consistency(BHJM_magnet_trimesh, **kw)
 
     Btest = [
         [1.54452002e-03, 3.11861149e-03, 4.68477835e-03],
@@ -453,22 +443,18 @@ def test_magnet_trimesh_field_BH():
     ]
     np.testing.assert_allclose(H, Htest)
 
+    Jtest = [(0, 0, 0), (3, 2, 1)]
+    np.testing.assert_allclose(J, Jtest, rtol=1e-06)
 
-def test_current_circle_field_BH():
+
+def test_BHJM_circle():
     """Test of current circle field core function"""
-    B = magpy.core.current_circle_field(
-        field="B",
-        observers=np.array([(1, 1, 1), (2, 2, 2), (3, 3, 3)]),
-        diameter=np.array([2, 4, 6]),
-        current=np.array([1, 1, 2]) * 1e3,
-    )
-    H = magpy.core.current_circle_field(
-        field="H",
-        observers=np.array([(1, 1, 1), (2, 2, 2), (3, 3, 3)]),
-        diameter=np.array([2, 4, 6]),
-        current=np.array([1, 1, 2]) * 1e3,
-    )
-    np.testing.assert_allclose(B, MU0 * H)
+    kw = {
+        "observers": np.array([(1, 1, 1), (2, 2, 2), (3, 3, 3)]),
+        "current": np.array([1, 1, 2]) * 1e3,
+        "diameter": np.array([2, 4, 6]),
+    }
+    H, B, M, _ = helper_check_HBMJ_consistency(BHJM_circle, **kw)
 
     Btest = (
         np.array(
@@ -494,25 +480,21 @@ def test_current_circle_field_BH():
     )
     np.testing.assert_allclose(H, Htest)
 
+    Mtest = [(0, 0, 0)] * 3
+    np.testing.assert_allclose(M, Mtest, rtol=1e-06)
+
 
-def test_current_polyline_field_BH():
+def test_BHJM_current_polyline():
     """Test of current polyline field core function"""
     vert = np.array([(-1.5, 0, 0), (-0.5, 0, 0), (0.5, 0, 0), (1.5, 0, 0)])
-    B = magpy.core.current_polyline_field(
-        field="B",
-        observers=np.array([(0, 0, 1)] * 3),
-        segment_start=vert[:-1],
-        segment_end=vert[1:],
-        current=np.array([1, 1, 1]),
-    )
-    H = magpy.core.current_polyline_field(
-        field="H",
-        observers=np.array([(0, 0, 1)] * 3),
-        segment_start=vert[:-1],
-        segment_end=vert[1:],
-        current=np.array([1, 1, 1]),
-    )
-    np.testing.assert_allclose(B, MU0 * H)
+
+    kw = {
+        "observers": np.array([(0, 0, 1)] * 3),
+        "current": np.array([1, 1, 1]),
+        "segment_start": vert[:-1],
+        "segment_end": vert[1:],
+    }
+    H, B, M, _ = helper_check_HBMJ_consistency(BHJM_current_polyline, **kw)
 
     Btest = (
         np.array(
@@ -538,24 +520,19 @@ def test_current_polyline_field_BH():
     )
     np.testing.assert_allclose(H, Htest, rtol=0, atol=1e-7)
 
+    Mtest = [(0, 0, 0)] * 3
+    np.testing.assert_allclose(M, Mtest, rtol=1e-06)
+
 
-def test_dipole_field_BH():
+def test_BHJM_dipole():
     """Test of dipole field core function"""
-    obs = np.array([(1, 2, 3), (-1, -2, -3), (3, 3, -1)])
     pol = np.array([(0, 0, 1), (1, 0, 1), (-1, 0.321, 0.123)])
-    mom = pol * 4 * np.pi / 3 / MU0
 
-    B = magpy.core.dipole_field(
-        field="B",
-        observers=obs,
-        moment=mom,
-    )
-    H = magpy.core.dipole_field(
-        field="H",
-        observers=obs,
-        moment=mom,
-    )
-    np.testing.assert_allclose(B, MU0 * H)
+    kw = {
+        "observers": np.array([(1, 2, 3), (-1, -2, -3), (3, 3, -1)]),
+        "moment": pol * 4 * np.pi / 3 / MU0,
+    }
+    H, B, M, _ = helper_check_HBMJ_consistency(BHJM_dipole, **kw)
 
     Btest = [
         [4.09073329e-03, 8.18146659e-03, 5.90883698e-03],
@@ -571,6 +548,9 @@ def test_dipole_field_BH():
     ]
     np.testing.assert_allclose(H, Htest)
 
+    Mtest = [(0, 0, 0)] * 3
+    np.testing.assert_allclose(M, Mtest, rtol=1e-06)
+
 
 #######################################################################################
 #######################################################################################
@@ -579,37 +559,6 @@ def test_dipole_field_BH():
 # OTHER TESTS AND V4 TESTS
 
 
-def test_line_deprecation():
-    """test line deprecation"""
-    kw = {
-        "field": "B",
-        "observers": np.array([(0, 0, 1)]),
-        "segment_start": np.array([(0, 0, 1)]),
-        "segment_end": np.array([(0, 0, 1)]),
-        "current": np.array([1]),
-    }
-    B1 = current_polyline_field(**kw)
-    with pytest.warns(MagpylibDeprecationWarning):
-        B2 = current_line_field(**kw)
-
-    np.testing.assert_allclose(B1, B2)
-
-
-def test_loop_deprecation():
-    """test loop deprecation"""
-    kw = {
-        "field": "B",
-        "observers": np.array([(0, 0, 1)]),
-        "diameter": np.array([1]),
-        "current": np.array([1]),
-    }
-    B1 = current_circle_field(**kw)
-    with pytest.warns(MagpylibDeprecationWarning):
-        B2 = current_loop_field(**kw)
-
-    np.testing.assert_allclose(B1, B2)
-
-
 def test_field_loop_specials():
     """test loop special cases"""
     cur = np.array([1, 1, 1, 1, 0, 2])
@@ -617,7 +566,7 @@ def test_field_loop_specials():
     obs = np.array([(0, 0, 0), (1, 0, 0), (0, 0, 0), (1, 0, 0), (1, 0, 0), (0, 0, 0)])
 
     B = (
-        current_circle_field(
+        BHJM_circle(
             field="B",
             observers=obs,
             diameter=dia,
@@ -646,7 +595,7 @@ def test_field_line_special_cases():
 
     # only normal
     B1 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po1,
             current=c1,
@@ -661,7 +610,7 @@ def test_field_line_special_cases():
     # only on_line
     po1b = np.array([(1, 1, 1)])
     B2 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po1b,
             current=c1,
@@ -675,7 +624,7 @@ def test_field_line_special_cases():
 
     # only zero-segment
     B3 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po1,
             current=c1,
@@ -693,7 +642,7 @@ def test_field_line_special_cases():
     pe2 = np.array([(0, 0, 0), (2, 2, 2)])
     po2 = np.array([(1, 2, 3), (1, 1, 1)])
     B4 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po2,
             current=c2,
@@ -708,7 +657,7 @@ def test_field_line_special_cases():
     # normal + zero_segment
     po2b = np.array([(1, 2, 3), (1, 2, 3)])
     B5 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po2b,
             current=c2,
@@ -723,7 +672,7 @@ def test_field_line_special_cases():
     # normal + on_line
     pe2b = np.array([(2, 2, 2)] * 2)
     B6 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po2,
             current=c2,
@@ -741,7 +690,7 @@ def test_field_line_special_cases():
     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 = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=po4,
             current=c4,
@@ -759,7 +708,7 @@ def test_field_loop2():
     curr = np.array([1])
     dia = np.array([2])
     obs = np.array([[0, 0, 0]])
-    B = current_circle_field(
+    B = BHJM_circle(
         field="B",
         observers=obs,
         current=curr,
@@ -769,7 +718,7 @@ def test_field_loop2():
     curr = np.array([1] * 2)
     dia = np.array([2] * 2)
     obs = np.array([[0, 0, 0]] * 2)
-    B2 = current_circle_field(
+    B2 = BHJM_circle(
         field="B",
         observers=obs,
         current=curr,
@@ -811,7 +760,7 @@ def test_field_line_from_vert():
         cu = np.array([curr] * (len(vert) - 1))
         B += [
             np.sum(
-                current_polyline_field(
+                BHJM_current_polyline(
                     field="B",
                     observers=po,
                     current=cu,
@@ -843,7 +792,7 @@ def test_field_line_v4():
         ]
     )
     B = (
-        current_polyline_field(
+        BHJM_current_polyline(
             field="B",
             observers=obs,
             current=cur,
@@ -902,7 +851,7 @@ def test_triangle5():
     ver = np.array([[(0, 0, 0), (0, 5, 0), (5, 0, 0)]] * n)
 
     b1 = (
-        magpy.core.triangle_field(
+        BHJM_triangle(
             field="H",
             observers=obs1,
             polarization=mag,
@@ -912,7 +861,7 @@ def test_triangle5():
     )
     np.testing.assert_allclose(btest1, b1)
     b2 = (
-        magpy.core.triangle_field(
+        BHJM_triangle(
             field="H",
             observers=obs2,
             polarization=mag,
@@ -930,19 +879,19 @@ def test_triangle6():
     obs3 = np.array([(5, 0, 0)])
     mag = np.array([(1, 2, 3)])
     ver = np.array([[(0, 0, 0), (0, 5, 0), (5, 0, 0)]])
-    b1 = magpy.core.triangle_field(
+    b1 = BHJM_triangle(
         field="B",
         observers=obs1,
         polarization=mag,
         vertices=ver,
     )
-    b2 = magpy.core.triangle_field(
+    b2 = BHJM_triangle(
         field="B",
         observers=obs2,
         polarization=mag,
         vertices=ver,
     )
-    b3 = magpy.core.triangle_field(
+    b3 = BHJM_triangle(
         field="B",
         observers=obs3,
         polarization=mag,
diff --git a/tests/test_core.py b/tests/test_core.py
new file mode 100644
index 000000000..befd37cfa
--- /dev/null
+++ b/tests/test_core.py
@@ -0,0 +1,15 @@
+# here all core functions should be tested properly - ideally against FEM
+import numpy as np
+
+from magpylib._src.fields.field_BH_sphere import magnet_sphere_Bfield
+
+
+def test_magnet_sphere_Bfield():
+    "magnet_sphere_Bfield test"
+    B = magnet_sphere_Bfield(
+        observers=np.array([(0, 0, 0)]),
+        diameters=np.array([1]),
+        polarizations=np.array([(0, 0, 1)]),
+    )
+    Btest = np.array([(0, 0, 2 / 3)])
+    np.testing.assert_allclose(B, Btest)
diff --git a/tests/test_core_physics_consistency.py b/tests/test_core_physics_consistency.py
index 610d90c72..ae29e6d4a 100644
--- a/tests/test_core_physics_consistency.py
+++ b/tests/test_core_physics_consistency.py
@@ -1,7 +1,16 @@
 import numpy as np
+from scipy.constants import mu_0 as MU0
+
+from magpylib._src.fields.field_BH_circle import BHJM_circle
+from magpylib._src.fields.field_BH_cuboid import BHJM_magnet_cuboid
+from magpylib._src.fields.field_BH_cylinder import BHJM_magnet_cylinder
+from magpylib._src.fields.field_BH_cylinder_segment import BHJM_cylinder_segment
+from magpylib._src.fields.field_BH_dipole import BHJM_dipole
+from magpylib._src.fields.field_BH_polyline import BHJM_current_polyline
+from magpylib._src.fields.field_BH_sphere import BHJM_magnet_sphere
+from magpylib._src.fields.field_BH_tetrahedron import BHJM_magnet_tetrahedron
+from magpylib._src.fields.field_BH_triangle import BHJM_triangle
 
-import magpylib as magpy
-from magpylib._src.utility import MU0
 
 # PHYSICS CONSISTENCY TESTING
 #
@@ -18,6 +27,8 @@
 # Geometric approximation testing should give similar results for different
 # implementations when one geometry is constructed from another
 #
+# Scaling invariance of solutions
+#
 # ----------> Circle          # webpage numbers
 # Circle   <> Dipole          # mom = I*A (far field approx)
 # Polyline <> Dipole          # mom = I*A (far field approx)
@@ -44,26 +55,26 @@ def test_core_phys_moment_of_current_circle():
     curr = np.array([1e3, 1e3])
     mom = ((dia / 2) ** 2 * np.pi * curr * np.array([(0, 0, 1)] * 2).T).T
 
-    B1 = magpy.core.current_circle_field(
+    B1 = BHJM_circle(
         field="B",
         observers=obs,
         diameter=dia,
         current=curr,
     )
-    B2 = magpy.core.dipole_field(
+    B2 = BHJM_dipole(
         field="B",
         observers=obs,
         moment=mom,
     )
     np.testing.assert_allclose(B1, B2, rtol=1e-02)
 
-    H1 = magpy.core.current_circle_field(
+    H1 = BHJM_circle(
         field="H",
         observers=obs,
         diameter=dia,
         current=curr,
     )
-    H2 = magpy.core.dipole_field(
+    H2 = BHJM_dipole(
         field="H",
         observers=obs,
         moment=mom,
@@ -84,7 +95,7 @@ def test_core_phys_moment_of_current_square():
     curr4 = np.array([curr1] * 4)
     mom = (4 * curr1 * np.array([(0, 0, 1)]).T).T
 
-    B1 = magpy.core.current_polyline_field(
+    B1 = BHJM_current_polyline(
         field="B",
         observers=obs4,
         segment_start=vert[:-1],
@@ -92,14 +103,14 @@ def test_core_phys_moment_of_current_square():
         current=curr4,
     )
     B1 = np.sum(B1, axis=0)
-    B2 = magpy.core.dipole_field(
+    B2 = BHJM_dipole(
         field="B",
         observers=obs1,
         moment=mom,
     )[0]
     np.testing.assert_allclose(B1, -B2, rtol=1e-03)
 
-    H1 = magpy.core.current_polyline_field(
+    H1 = BHJM_current_polyline(
         field="H",
         observers=obs4,
         segment_start=vert[:-1],
@@ -107,7 +118,7 @@ def test_core_phys_moment_of_current_square():
         current=curr4,
     )
     H1 = np.sum(H1, axis=0)
-    H2 = magpy.core.dipole_field(
+    H2 = BHJM_dipole(
         field="H",
         observers=obs1,
         moment=mom,
@@ -126,13 +137,13 @@ def test_core_phys_circle_polyline():
     obs99 = np.array([(1, 2, 3)] * 299)
     dia = np.array([2])
 
-    H1 = magpy.core.current_circle_field(
+    H1 = BHJM_circle(
         field="H",
         observers=obs,
         diameter=dia,
         current=curr,
     )[0]
-    H2 = magpy.core.current_polyline_field(
+    H2 = BHJM_current_polyline(
         field="H",
         observers=obs99,
         segment_start=vert[:-1],
@@ -142,13 +153,13 @@ def test_core_phys_circle_polyline():
     H2 = np.sum(H2, axis=0)
     np.testing.assert_allclose(H1, -H2, rtol=1e-4)
 
-    B1 = magpy.core.current_circle_field(
+    B1 = BHJM_circle(
         field="B",
         observers=obs,
         diameter=dia,
         current=curr,
     )[0]
-    B2 = magpy.core.current_polyline_field(
+    B2 = BHJM_current_polyline(
         field="B",
         observers=obs99,
         segment_start=vert[:-1],
@@ -171,31 +182,31 @@ def test_core_physics_dipole_sphere():
     pol = np.array([(1, 2, 3), (0, 0, 1), (-1, -2, 0), (1, -1, 0.1)])
     mom = np.array([4 * (d / 2) ** 3 * np.pi / 3 * p / MU0 for d, p in zip(dia, pol)])
 
-    B1 = magpy.core.magnet_sphere_field(
+    B1 = BHJM_magnet_sphere(
         field="B",
         observers=obs,
         diameter=dia,
         polarization=pol,
     )
-    B2 = magpy.core.dipole_field(
+    B2 = BHJM_dipole(
         field="B",
         observers=obs,
         moment=mom,
     )
     np.testing.assert_allclose(B1, B2, rtol=0, atol=1e-16)
 
-    H1 = magpy.core.magnet_sphere_field(
+    H1 = BHJM_magnet_sphere(
         field="H",
         observers=obs,
         diameter=dia,
         polarization=pol,
     )
-    H2 = magpy.core.dipole_field(
+    H2 = BHJM_dipole(
         field="H",
         observers=obs,
         moment=mom,
     )
-    np.testing.assert_allclose(H1, H2, rtol=0, atol=1e-11)
+    np.testing.assert_allclose(H1, H2, rtol=0, atol=1e-10)
 
 
 # -> Circle, Cylinder
@@ -221,7 +232,7 @@ def test_core_physics_long_solenoid():
             BHz_long *= MU0
 
         # SOLENOID TEST constructed from circle fields
-        BH = magpy.core.current_circle_field(
+        BH = BHJM_circle(
             field=field,
             observers=np.linspace((0, 0, -L / 2), (0, 0, L / 2), N),
             diameter=np.array([2 * R] * N),
@@ -238,7 +249,7 @@ def test_core_physics_long_solenoid():
         obs = np.array([(0, 0, 0)])
 
         # cylinder
-        BHz_cyl = magpy.core.magnet_cylinder_field(
+        BHz_cyl = BHJM_magnet_cylinder(
             field=field,
             observers=obs,
             dimension=np.array([(2 * R, L)]),
@@ -249,7 +260,7 @@ def test_core_physics_long_solenoid():
         np.testing.assert_allclose(BHz_long, BHz_cyl, rtol=1e-5)
 
         # cuboid
-        BHz_cub = magpy.core.magnet_cuboid_field(
+        BHz_cub = BHJM_magnet_cuboid(
             field=field,
             observers=obs,
             dimension=np.array([(2 * R, 2 * R, L)]),
@@ -275,7 +286,7 @@ def test_core_physics_current_replacement():
     obs = np.array([(1.5, -2, -1.123)])
 
     Jz = 1
-    Hz_cyl = magpy.core.magnet_cylinder_field(
+    Hz_cyl = BHJM_magnet_cylinder(
         field="H",
         observers=obs,
         dimension=np.array([(2 * R, L)]),
@@ -284,7 +295,7 @@ def test_core_physics_current_replacement():
 
     N = 1000  # current discretization
     I = Jz / MU0 / N * L
-    H = magpy.core.current_circle_field(
+    H = BHJM_circle(
         field="H",
         observers=np.linspace((0, 0, -L / 2), (0, 0, L / 2), N) + obs,
         diameter=np.array([2 * R] * N),
@@ -303,7 +314,7 @@ def test_core_physics_geometry_cylinder_from_segments():
     obs = np.array([(1, 2, 3), (0.23, 0.132, 0.123)])
     pol = np.array([(2, 0.123, 3), (-0.23, -1, 0.434)])
 
-    B_cyl = magpy.core.magnet_cylinder_field(
+    B_cyl = BHJM_magnet_cylinder(
         field="B",
         observers=obs,
         dimension=np.array([(2 * r, h)] * 2),
@@ -313,7 +324,7 @@ def test_core_physics_geometry_cylinder_from_segments():
 
     Bseg = np.zeros((2, 3))
     for phi1, phi2 in zip(sections[:-1], sections[1:]):
-        B_part = magpy.core.magnet_cylinder_segment_field(
+        B_part = BHJM_cylinder_segment(
             field="B",
             observers=obs,
             dimension=np.array([(0, r, h, phi1, phi2)] * 2),
@@ -330,7 +341,7 @@ def test_core_physics_dipole_approximation_magnet_far_field():
     obs = np.array([(100, 200, 300), (-200, -200, -200)])
 
     mom = np.array([(1e6, 2e6, 3e6)] * 2)
-    Bdip = magpy.core.dipole_field(
+    Bdip = BHJM_dipole(
         field="H",
         observers=obs,
         moment=mom,
@@ -339,7 +350,7 @@ def test_core_physics_dipole_approximation_magnet_far_field():
     dim = np.array([(2, 2, 2)] * 2)
     vol = 8
     pol = mom / vol * MU0
-    Bcub = magpy.core.magnet_cuboid_field(
+    Bcub = BHJM_magnet_cuboid(
         field="H",
         observers=obs,
         dimension=dim,
@@ -350,7 +361,7 @@ def test_core_physics_dipole_approximation_magnet_far_field():
     dim = np.array([(0.5, 0.5)] * 2)
     vol = 0.25**2 * np.pi * 0.5
     pol = mom / vol * MU0
-    Bcyl = magpy.core.magnet_cylinder_field(
+    Bcyl = BHJM_magnet_cylinder(
         field="H",
         observers=obs,
         dimension=dim,
@@ -361,7 +372,7 @@ def test_core_physics_dipole_approximation_magnet_far_field():
     vert = np.array([[(0, 0, 0), (0, 0, 0.1), (0.1, 0, 0), (0, 0.1, 0)]] * 2)
     vol = 1 / 6 * 1e-3
     pol = mom / vol * MU0
-    Btetra = magpy.core.magnet_tetrahedron_field(
+    Btetra = BHJM_magnet_tetrahedron(
         field="H",
         observers=obs,
         vertices=vert,
@@ -372,7 +383,7 @@ def test_core_physics_dipole_approximation_magnet_far_field():
     dim = np.array([(0.1, 0.2, 0.1, -25, 25)] * 2)
     vol = 3 * np.pi * (50 / 360) * 1e-3
     pol = mom / vol * MU0
-    Bcys = magpy.core.magnet_cylinder_segment_field(
+    Bcys = BHJM_cylinder_segment(
         field="H",
         observers=obs + np.array((0.15, 0, 0)),
         dimension=dim,
@@ -395,7 +406,7 @@ def test_core_physics_circle_VS_webpage_numbers():
     Hz = [500, 176.8, 44.72, 15.81]
     Htest = [(0, 0, hz) for hz in Hz]
 
-    H = magpy.core.current_circle_field(
+    H = BHJM_circle(
         field="H",
         observers=obs,
         diameter=dia,
@@ -412,7 +423,7 @@ def test_core_physics_circle_VS_webpage_numbers():
     ]
     Btest = [(0, 0, bz) for bz in Bz]
 
-    B = magpy.core.current_circle_field(
+    B = BHJM_circle(
         field="B",
         observers=obs,
         diameter=dia,
@@ -429,7 +440,7 @@ def test_core_physics_cube_current_replacement():
     Jz = 1.23
     dim = np.array([(2, 2, h)] * 2)
     pol = np.array([(0, 0, Jz)] * 2)
-    Hcub = magpy.core.magnet_cuboid_field(
+    Hcub = BHJM_magnet_cuboid(
         field="H",
         observers=obs,
         dimension=dim,
@@ -451,7 +462,7 @@ def test_core_physics_cube_current_replacement():
 
     Hcurr = np.zeros((2, 3))
     for i, obss in enumerate([obs1, obs2]):
-        h = magpy.core.current_polyline_field(
+        h = BHJM_current_polyline(
             field="H",
             observers=obss,
             segment_start=start,
@@ -475,7 +486,7 @@ def test_core_physics_triangle_cube_geometry():
             [(1, 1, -1), (1, -1, -1), (-1, 1, -1)],  # bott2
         ]
     )
-    b = magpy.core.triangle_field(
+    b = BHJM_triangle(
         field="B",
         observers=obs,
         vertices=fac,
@@ -486,7 +497,7 @@ def test_core_physics_triangle_cube_geometry():
     obs = np.array([(3, 4, 5)])
     mag = np.array([(0, 0, 333)])
     dim = np.array([(2, 2, 2)])
-    bb = magpy.core.magnet_cuboid_field(
+    bb = BHJM_magnet_cuboid(
         field="B",
         observers=obs,
         dimension=dim,
@@ -505,7 +516,7 @@ def test_core_physics_triangle_VS_itself():
             [(0, 0, 0), (10, 0, 0), (0, 10, 0)],
         ]
     )
-    b = magpy.core.triangle_field(
+    b = BHJM_triangle(
         field="B",
         observers=obs,
         polarization=mag,
@@ -523,7 +534,7 @@ def test_core_physics_triangle_VS_itself():
             [(6, 0, 0), (10, 0, 0), (0, 10, 0)],
         ]
     )
-    bb = magpy.core.triangle_field(
+    bb = BHJM_triangle(
         field="B",
         observers=obs,
         polarization=mag,
@@ -577,13 +588,13 @@ def test_core_physics_Tetrahedron_VS_Cuboid():
         for obs in obss:
             obs6 = np.tile(obs, (6, 1))
             mag6 = np.tile(mag, (6, 1))
-            b = magpy.core.magnet_tetrahedron_field(
+            b = BHJM_magnet_tetrahedron(
                 field="B",
                 observers=obs6,
                 polarization=mag6,
                 vertices=ver,
             )
-            h = magpy.core.magnet_tetrahedron_field(
+            h = BHJM_magnet_tetrahedron(
                 field="H",
                 observers=obs6,
                 polarization=mag6,
@@ -595,13 +606,13 @@ def test_core_physics_Tetrahedron_VS_Cuboid():
             obs1 = np.reshape(obs, (1, 3))
             mag1 = np.reshape(mag, (1, 3))
             dim = np.array([(2, 2, 2)])
-            bb = magpy.core.magnet_cuboid_field(
+            bb = BHJM_magnet_cuboid(
                 field="B",
                 observers=obs1,
                 polarization=mag1,
                 dimension=dim,
             )[0]
-            hh = magpy.core.magnet_cuboid_field(
+            hh = BHJM_magnet_cuboid(
                 field="H",
                 observers=obs1,
                 polarization=mag1,
diff --git a/tests/test_display_plotly.py b/tests/test_display_plotly.py
index f6e5256ef..b9f51abfd 100644
--- a/tests/test_display_plotly.py
+++ b/tests/test_display_plotly.py
@@ -317,7 +317,7 @@ def test_subplots():
 
     # bad output value
     with pytest.raises(
-        ValueError, match=r"The `output` parameter must start with 'B' or 'H'.*"
+        ValueError, match=r"The `output` parameter must start with 'B', 'H', 'M', 'J'.*"
     ):
         magpy.show(*objs, canvas=fig, output="bad_output")
 
diff --git a/tests/test_exceptions.py b/tests/test_exceptions.py
index 3c9c0ffb8..07628ec58 100644
--- a/tests/test_exceptions.py
+++ b/tests/test_exceptions.py
@@ -10,7 +10,13 @@
 from magpylib._src.utility import format_obj_input
 from magpylib._src.utility import format_src_inputs
 
-GETBH_KWARGS = {"sumup": False, "squeeze": True, "pixel_agg": None, "output": "ndarray"}
+GETBH_KWARGS = {
+    "sumup": False,
+    "squeeze": True,
+    "pixel_agg": None,
+    "output": "ndarray",
+    "in_out": "auto",
+}
 
 
 def getBHv_unknown_source_type():
@@ -36,6 +42,7 @@ def getBH_level2_bad_input1():
         sumup=False,
         squeeze=True,
         pixel_agg=None,
+        in_out="auto",
         field="B",
         output="ndarray",
     )
diff --git a/tests/test_field_cylinder.py b/tests/test_field_cylinder.py
index 8aa13eb5c..b634feeb2 100644
--- a/tests/test_field_cylinder.py
+++ b/tests/test_field_cylinder.py
@@ -5,14 +5,17 @@
 import pytest
 
 import magpylib as magpy
-from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_core
-from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_field
+from magpylib._src.fields.field_BH_cylinder import BHJM_magnet_cylinder
+from magpylib._src.fields.field_BH_cylinder_segment import BHJM_cylinder_segment
+from magpylib._src.fields.field_BH_cylinder_segment import (
+    magnet_cylinder_segment_Hfield,
+)
 
 # pylint: disable="pointless-string-statement"
 # creating test data
 """ import os
 import numpy as np
-from magpylib._src.fields.field_BH_cylinder_tile import magnet_cylinder_segment_core
+from magpylib._src.fields.field_BH_cylinder_tile import magnet_cylinder_segment_Hfield
 
 N = 1111
 null = np.zeros(N)
@@ -38,9 +41,13 @@
 r = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(
+    magnetizations=MAG,
+    dimensions=dim,
+    observers=obs_pos,
+)
 DATA["cases [112, 212, 132, 232]"] = {
-    "inputs": {"mag": MAG, "dim": dim, "obs_pos": obs_pos},
+    "inputs": {"magnetizations": MAG, "dimensions": dim, "observers": obs_pos},
     "H_expected": H1,
 }
 
@@ -52,7 +59,7 @@
 r = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -66,7 +73,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -81,7 +88,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -95,7 +102,7 @@
 r = r2
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -110,7 +117,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -124,7 +131,7 @@
 r = r1
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -140,7 +147,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -155,7 +162,7 @@
 r1 = null
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -169,7 +176,7 @@
 r = r2
 obs_pos = np.array([r, phi, z]).T
 dim = np.array([r1, r2, phi1, phi2, z1, z2]).T
-H1 = magnet_cylinder_segment_core(mag=MAG, dim=dim, obs_pos=obs_pos)
+H1 = magnet_cylinder_segment_Hfield(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,
@@ -191,8 +198,15 @@
 )
 def test_cylinder_tile_slanovc(inputs, H_expected):
     "testing precomputed cylinder test cases"
-    H = magnet_cylinder_segment_core(**inputs)
-    assert np.allclose(np.nan_to_num(H), np.nan_to_num(H_expected))
+    inputs_mod = {
+        "magnetizations": inputs["mag"],
+        "observers": inputs["obs_pos"],
+        "dimensions": inputs["dim"],
+    }
+    H = (
+        magnet_cylinder_segment_Hfield(**inputs_mod) / 4 / np.pi * 1e7
+    )  # factors come from B <->H change
+    np.testing.assert_allclose(H, H_expected)
 
 
 def test_cylinder_field1():
@@ -206,7 +220,7 @@ def test_cylinder_field1():
     eins = np.ones(N)
     d, h, _ = dim.T
     dim5 = np.array([nulll, d / 2, h, nulll, eins * 360]).T
-    B1 = magnet_cylinder_segment_field(
+    B1 = BHJM_cylinder_segment(
         field="B", observers=poso, polarization=magg, dimension=dim5
     )
 
@@ -219,50 +233,47 @@ def test_cylinder_slanovc_field2():
         polarization=(22, 33, 44), dimension=(0.5, 1, 2, 0, 90)
     )
 
-    binn = (5.52525937, 13.04561569, 40.11111556)
-    bout = (0.0177018, 0.1277188, 0.27323195)
-    nulll = (0, 0, 0)
+    r_in = (0.5, 0.6, 0.3)
+    r_out = (1, 2, 3)
+    r_corn = (1, 0, 0)
+
+    b_in = (5.52525937, 13.04561569, 40.11111556)
+    b_out = (0.0177018, 0.1277188, 0.27323195)
+    b_corn = (0, 0, 0)
 
     # only inside
-    btest = np.array([binn] * 3)
-    B = src.getB([[0.5, 0.6, 0.3]] * 3)
+    btest = np.array([b_in] * 3)
+    B = src.getB([r_in] * 3)
     assert np.allclose(B, btest)
 
-    # only surf
-    btest = np.array([nulll] * 3)
-    B = src.getB([[1, 0, 0]] * 3)
+    # only edge
+    btest = np.array([b_corn] * 3)
+    B = src.getB([r_corn] * 3)
     assert np.allclose(B, btest)
 
     # only outside
-    btest = np.array([bout] * 3)
-    B = src.getB([[1, 2, 3]] * 3)
+    btest = np.array([b_out] * 3)
+    B = src.getB([r_out] * 3)
     assert np.allclose(B, btest)
 
-    # surf + out
-    btest = np.array([nulll, nulll, bout])
-    B = src.getB([0.6, 0, 1], [1, 0, 0.5], [1, 2, 3])
+    # edge + out
+    btest = np.array([b_corn, b_corn, b_out])
+    B = src.getB([r_corn, r_corn, r_out])
     assert np.allclose(B, btest)
 
     # surf + in
-    btest = np.array([nulll, nulll, binn])
-    B = src.getB([0, 0.5, 1], [1, 0, 0.5], [0.5, 0.6, 0.3])
+    btest = np.array([b_corn, b_corn, b_in])
+    B = src.getB(r_corn, r_corn, r_in)
     assert np.allclose(B, btest)
 
     # in + out
-    btest = np.array([bout, binn])
-    B = src.getB([1, 2, 3], [0.5, 0.6, 0.3])
+    btest = np.array([b_out, b_in])
+    B = src.getB(r_out, r_in)
     assert np.allclose(B, btest)
 
     # in + out + surf
-    btest = np.array([nulll, nulll, binn, bout, nulll, nulll])
-    B = src.getB(
-        [0.5, 0.5, 1],
-        [0, 1, 0.5],
-        [0.5, 0.6, 0.3],
-        [1, 2, 3],
-        [0.5, 0.6, -1],
-        [0, 1, -0.3],
-    )
+    btest = np.array([b_corn, b_corn, b_in, b_out, b_corn, b_corn])
+    B = src.getB([r_corn, r_corn, r_in, r_out, r_corn, r_corn])
     assert np.allclose(B, btest)
 
 
@@ -468,7 +479,7 @@ 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(
+    B = BHJM_magnet_cylinder(
         field="B",
         observers=np.array([(x, 0, 3) for x in np.logspace(-1.4, -1.2, 1000)]),
         polarization=np.array([(1, 1, 0)] * 1000),
@@ -490,13 +501,13 @@ def test_cyl_vs_cylseg_axial_H_inside_mask():
     dims = np.array([(1, 1)])
     dims_cs = np.array([(0, 0.5, 1, 0, 360)])
 
-    Bc = magpy.core.magnet_cylinder_field(
+    Bc = BHJM_magnet_cylinder(
         field=field,
         observers=obs,
         dimension=dims,
         polarization=pols,
     )
-    Bcs = magpy.core.magnet_cylinder_segment_field(
+    Bcs = BHJM_cylinder_segment(
         field=field,
         observers=obs,
         dimension=dims_cs,
diff --git a/tests/test_getBH_level2.py b/tests/test_getBH_level2.py
index 43dad817f..d34e07076 100644
--- a/tests/test_getBH_level2.py
+++ b/tests/test_getBH_level2.py
@@ -1,3 +1,6 @@
+import unittest
+import warnings
+
 import numpy as np
 import pytest
 
@@ -505,3 +508,66 @@ def test_pixel_agg3():
                     rtol=1e-5,
                     atol=1e-8,
                 )
+
+
+##############################################################
+def warnme1():
+    """test if in_out warning is thrown"""
+    sp = magpy.magnet.Sphere(
+        polarization=(1, 2, 3),
+        diameter=1,
+    )
+    sp.getB((1, 1, 1), in_out="inside")
+
+
+def warnme2():
+    """test if in_out warning is thrown"""
+    sp = magpy.magnet.Sphere(
+        polarization=(1, 2, 3),
+        diameter=1,
+    )
+    magpy.getH([sp, sp], (1, 1, 1), in_out="inside")
+
+
+class TestExceptions(unittest.TestCase):
+    """test class for exception testing"""
+
+    def test_warning(self):
+        """whatever"""
+        self.assertWarns(UserWarning, warnme1)
+        self.assertWarns(UserWarning, warnme2)
+
+
+##############################################################
+
+
+def do_not_warnme1():
+    """test if in_out warning is thrown"""
+    sp = magpy.magnet.Sphere(
+        polarization=(1, 2, 3),
+        diameter=1,
+    )
+    tetra = magpy.magnet.Tetrahedron(
+        polarization=(1, 2, 3),
+        vertices=[(1, 2, 3), (0, 0, 0), (1, 0, 0), (0, 1, 0)],
+    )
+    magpy.getH([sp, tetra], (1, 1, 1), in_out="inside")
+
+
+def do_not_warnme2():
+    """test if in_out warning is thrown"""
+    tetra = magpy.magnet.Tetrahedron(
+        polarization=(1, 2, 3),
+        vertices=[(1, 2, 3), (0, 0, 0), (1, 0, 0), (0, 1, 0)],
+    )
+    magpy.getH(tetra, (1, 1, 1), in_out="inside")
+
+
+def test_do_not_warn():
+    """do not warn"""
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        do_not_warnme1()
+        do_not_warnme2()
+        if len(w) > 0:
+            pytest.fail("WARNING SHOULD NOT HAVE BEEN RAISED")
diff --git a/tests/test_input_checks.py b/tests/test_input_checks.py
index 1def8caf1..f84cf0f4d 100644
--- a/tests/test_input_checks.py
+++ b/tests/test_input_checks.py
@@ -6,7 +6,7 @@
 from magpylib._src.exceptions import MagpylibBadUserInput
 from magpylib._src.exceptions import MagpylibDeprecationWarning
 from magpylib._src.exceptions import MagpylibMissingInput
-
+from magpylib._src.fields.field_BH_dipole import BHJM_dipole
 
 # pylint: disable=unnecessary-lambda-assignment
 
@@ -936,7 +936,7 @@ def test_input_getBH_field_good(field):
     """good getBH field inputs"""
     moms = np.array([[1, 2, 3]])
     obs = np.array([[1, 2, 3]])
-    B = magpy.core.dipole_field(field=field, observers=obs, moment=moms)
+    B = BHJM_dipole(field=field, observers=obs, moment=moms)
     assert isinstance(B, np.ndarray)
 
 
@@ -960,7 +960,7 @@ def test_input_getBH_field_bad(field):
     moms = np.array([[1, 2, 3]])
     obs = np.array([[1, 2, 3]])
     with pytest.raises(MagpylibBadUserInput):
-        magpy.core.dipole_field(field=field, observers=obs, moment=moms)
+        BHJM_dipole(field=field, observers=obs, moment=moms)
 
 
 def test_sensor_handedness():
diff --git a/tests/test_obj_BaseGeo.py b/tests/test_obj_BaseGeo.py
index c48038cc6..a31d91f4d 100644
--- a/tests/test_obj_BaseGeo.py
+++ b/tests/test_obj_BaseGeo.py
@@ -7,7 +7,6 @@
 import magpylib as magpy
 from magpylib._src.obj_classes.class_BaseGeo import BaseGeo
 
-
 # pylint: disable=no-member
 
 
@@ -410,104 +409,116 @@ def test_copy_order():
     assert desc_after == desc_before
 
 
-def test_describe():
-    """testing descibe method"""
+def match_string_up_to_id(s1, s2, /):
+    """Checks if 2 inputs (first is list, second string) match as long as id= follows a
+    series of numbers"""
+    patt = "id=[0-9]*[0-9]", "id=Regex"
+    assert re.sub(*patt, "\n".join(s1)).split("\n") == re.sub(*patt, s2).split("\n")
+
+
+def test_describe_with_label():
+    """testing describe method"""
+    # print("test = [\n    " + '",\n    '.join(f'"{s}' for s in desc.split("\n")) + '",\n]')
     # pylint: disable=protected-access
-    x1 = magpy.magnet.Cuboid(style_label="x1")
-    x2 = magpy.magnet.Cylinder(
-        style_label="x2", dimension=(1, 3), polarization=(2, 3, 4)
-    )
-    s1 = magpy.Sensor(position=[(1, 2, 3)] * 3, pixel=[(1, 2, 3)] * 15)
+    x = magpy.magnet.Cuboid(style_label="x1")
 
-    desc = x1.describe()
+    # describe calls print by default -> no return value
+    desc = x.describe()
     assert desc is None
 
-    test = (
-        "<pre>Cuboid(id=REGEX, label='x1')<br>  • parent: None <br>  • "
-        "position: [0. 0. 0.] m<br>  • orientation: [0. 0. 0.] deg<br>  • "
-        "dimension: None m<br>  • magnetization: None A/m<br>  • polarization: None T</pre>"
-    )
-    rep = x1._repr_html_()
-    rep = re.sub("id=[0-9]*[0-9]", "id=REGEX", rep)
-    assert test == rep
-
-    magpy.Collection(x1, x2)
+    # describe string
     test = [
-        "Cuboid(id=REGEX, label='x1')",
-        "  • parent: Collection(id=REGEX) ",  # INVISIBLE SPACE
+        "Cuboid(id=2743358534352, label='x1')",
+        "  • parent: None",
         "  • position: [0. 0. 0.] m",
         "  • orientation: [0. 0. 0.] deg",
         "  • dimension: None m",
         "  • magnetization: None A/m",
         "  • polarization: None T",
     ]
-    desc = x1.describe(return_string=True)
-    desc = re.sub("id=*[0-9]*[0-9]", "id=REGEX", desc)
-    assert test == desc.split("\n")
+    match_string_up_to_id(test, x.describe(return_string=True))
+
+    # describe html string
+    test_html = ("<pre>" + "\n".join(test) + "</pre>").split("\n")
+    match_string_up_to_id(test_html, x._repr_html_().replace("<br>", "\n"))
 
+
+def test_describe_with_parent():
+    """testing describe method"""
+    # print("test = [\n    " + '",\n    '.join(f'"{s}' for s in desc.split("\n")) + '",\n]')
+    x = magpy.magnet.Cuboid(style_label="x1")
+    magpy.Collection(x)  # add parent
     test = [
-        "Cylinder(id=REGEX, label='x2')",
-        "  • parent: Collection(id=REGEX) ",  # INVISIBLE SPACE
+        "Cuboid(id=1687262797456, label='x1')",
+        "  • parent: Collection(id=1687262859280)",
         "  • position: [0. 0. 0.] m",
         "  • orientation: [0. 0. 0.] deg",
-        "  • dimension: [1. 3.] m",
-        "  • magnetization: [1591549.43091895 2387324.14637843 3183098.86183791] A/m",
-        "  • polarization: [2. 3. 4.] T",
+        "  • dimension: None m",
+        "  • magnetization: None A/m",
+        "  • polarization: None T",
     ]
-    desc = x2.describe(return_string=True)
-    desc = re.sub("id=*[0-9]*[0-9]", "id=REGEX", desc)
-    assert test == desc.split("\n")
+    match_string_up_to_id(test, x.describe(return_string=True))
 
+
+def test_describe_with_path():
+    """testing describe method"""
+    # print("test = [\n    " + '",\n    '.join(f'"{s}' for s in desc.split("\n")) + '",\n]')
+    x = magpy.Sensor(position=[(1, 2, 3)] * 3)
     test = [
-        "Sensor(id=REGEX)",
-        "  • parent: None ",  # INVISIBLE SPACE
+        "Sensor(id=2743359152656)",
+        "  • parent: None",
         "  • path length: 3",
         "  • position (last): [1. 2. 3.] m",
         "  • orientation (last): [0. 0. 0.] deg",
-        "  • handedness: right ",
-        "  • pixel: 15 ",  # INVISIBLE SPACE
+        "  • handedness: right",
+        "  • pixel: None",
     ]
-    desc = s1.describe(return_string=True)
-    desc = re.sub("id=*[0-9]*[0-9]", "id=REGEX", desc)
-    assert test == desc.split("\n")
-
-    # exclude=None test
-    s = magpy.Sensor()
-    desc = s.describe(exclude=None, return_string=True)
-    test = (
-        "Sensor(id=REGEX)\n"
-        + "  • parent: None \n"
-        + "  • position: [0. 0. 0.] m\n"
-        + "  • orientation: [0. 0. 0.] deg\n"
-        + "  • handedness: right \n"
-        + "  • pixel: None \n"
-        + "  • style: SensorStyle(arrows=ArrowCS(x=ArrowSingle(color=None, show=True), "
-        + "y=ArrowSingle(color=None, show=True), z=ArrowSingle(color=None, show=True)),"
-        + " color=None, description=Description(show=None, text=None), label=None, "
-        + "legend=Legend(show=None, text=None), "
-        + "model3d=Model3d(data=[], showdefault=True), opacity=None, path=Path(frames=None,"
-        + " line=Line(color=None, style=None, width=None), marker=Marker(color=None,"
-        + " size=None, symbol=None), numbering=None, show=None), pixel=Pixel(color=None,"
-        + " size=1, sizemode=None, symbol=None), size=None, sizemode=None) "
-    )
-    desc = re.sub("id=*[0-9]*[0-9]", "id=REGEX", desc)
-    assert desc == test
-
-    # lots of sensor pixel
-    s = magpy.Sensor(pixel=[[[(1, 2, 3)] * 5] * 5] * 3, handedness="left")
-    desc = s.describe(return_string=True)
-    test = (
-        "Sensor(id=REGEX)\n"
-        + "  • parent: None \n"
-        + "  • position: [0. 0. 0.] m\n"
-        + "  • orientation: [0. 0. 0.] deg\n"
-        + "  • handedness: left \n"
-        + "  • pixel: 75 (3x5x5) "
-    )
-    desc = re.sub("id=*[0-9]*[0-9]", "id=REGEX", desc)
-    assert desc == test
+    match_string_up_to_id(test, x.describe(return_string=True))
+
 
-    # describe tringularmesh
+def test_describe_with_exclude_None():
+    """testing describe method"""
+    # print("test = [\n    " + '",\n    '.join(f'"{s}' for s in desc.split("\n")) + '",\n]')
+    x = magpy.Sensor()
+    test = [
+        "Sensor(id=1687262758416)",
+        "  • parent: None",
+        "  • position: [0. 0. 0.] m",
+        "  • orientation: [0. 0. 0.] deg",
+        "  • handedness: right",
+        "  • pixel: None",
+        (
+            "  • style: SensorStyle(arrows=ArrowCS(x=ArrowSingle(color=None, show=True),"
+            " y=ArrowSingle(color=None, show=True), z=ArrowSingle(color=None, show=True)),"
+            " color=None, description=Description(show=None, text=None), label=None,"
+            " legend=Legend(show=None, text=None), model3d=Model3d(data=[], showdefault=True),"
+            " opacity=None, path=Path(frames=None, line=Line(color=None, style=None, width=None),"
+            " marker=Marker(color=None, size=None, symbol=None), numbering=None, show=None),"
+            " pixel=Pixel(color=None, size=1, sizemode=None, symbol=None), size=None,"
+            " sizemode=None)"
+        ),
+    ]
+    match_string_up_to_id(test, x.describe(exclude=None, return_string=True))
+
+
+def test_describe_with_many_pixels():
+    """testing describe method"""
+    # print("test = [\n    " + '",\n    '.join(f'"{s}' for s in desc.split("\n")) + '",\n]')
+    x = magpy.Sensor(pixel=[[[(1, 2, 3)] * 5] * 5] * 3, handedness="left")
+    test = [
+        "Sensor(id=1687262996944)",
+        "  • parent: None",
+        "  • position: [0. 0. 0.] m",
+        "  • orientation: [0. 0. 0.] deg",
+        "  • handedness: left",
+        "  • pixel: 75 (3x5x5)",
+    ]
+    match_string_up_to_id(test, x.describe(return_string=True))
+
+
+def test_describe_with_triangularmesh():
+    """testing describe method"""
+    # print("test = [\n    " + '",\n    '.join(f'"{s}' for s in desc.split("\n")) + '",\n]')
     points = [
         (-1, -1, 0),
         (-1, 1, 0),
@@ -515,34 +526,32 @@ def test_describe():
         (1, 1, 0),
         (0, 0, 2),
     ]
-    s = magpy.magnet.TriangularMesh.from_ConvexHull(
+    x = magpy.magnet.TriangularMesh.from_ConvexHull(
         polarization=(0, 0, 1),
         points=points,
         check_selfintersecting="skip",
     )
-    desc = s.describe(return_string=True)
-    test = (
-        "TriangularMesh(id=REGEX)\n"
-        "  • parent: None \n"
-        "  • position: [0. 0. 0.] m\n"
-        "  • orientation: [0. 0. 0.] deg\n"
-        "  • magnetization: [     0.              0.         795774.71545948] A/m\n"
-        "  • polarization: [0. 0. 1.] T\n"
-        "  • barycenter: [0.         0.         0.46065534] \n"
-        "  • faces: shape(6, 3) \n"
-        "  • mesh: shape(6, 3, 3) \n"
-        "  • status_disconnected: False \n"
-        "  • status_disconnected_data: 1 part \n"
-        "  • status_open: False \n"
-        "  • status_open_data: [] \n"
-        "  • status_reoriented: True \n"
-        "  • status_selfintersecting: None \n"
-        "  • status_selfintersecting_data: None \n"
-        "  • vertices: shape(5, 3) "
-    )
-    desc = re.sub("id=*[0-9]*[0-9]", "id=REGEX", desc)
-    # to create test: print('\\n"\n'.join(f'"{s}' for s in desc.split("\n")) + '"')
-    assert desc == test
+    test = [
+        "TriangularMesh(id=1687257413648)",
+        "  • parent: None",
+        "  • position: [0. 0. 0.] m",
+        "  • orientation: [0. 0. 0.] deg",
+        "  • magnetization: [     0.              0.         795774.71545948] A/m",
+        "  • polarization: [0. 0. 1.] T",
+        "  • barycenter: [0.         0.         0.46065534]",
+        "  • faces: shape(6, 3)",
+        "  • mesh: shape(6, 3, 3)",
+        "  • status_disconnected: False",
+        "  • status_disconnected_data: 1 part",
+        "  • status_open: False",
+        "  • status_open_data: []",
+        "  • status_reoriented: True",
+        "  • status_selfintersecting: None",
+        "  • status_selfintersecting_data: None",
+        "  • vertices: shape(5, 3)",
+    ]
+
+    match_string_up_to_id(test, x.describe(return_string=True))
 
 
 def test_unset_describe():
diff --git a/tests/test_obj_Circle.py b/tests/test_obj_Circle.py
index 711474b75..5875560f8 100644
--- a/tests/test_obj_Circle.py
+++ b/tests/test_obj_Circle.py
@@ -15,7 +15,7 @@ def test_Circle_basic_B():
     assert np.allclose(B, Btest)
 
 
-def test_current_circle_field():
+def test_BHJM_circle():
     """test explicit field output values"""
     s = magpy.current.Circle(current=1, diameter=1)
 
diff --git a/tests/test_obj_Collection.py b/tests/test_obj_Collection.py
index f50121533..8db0a32bb 100644
--- a/tests/test_obj_Collection.py
+++ b/tests/test_obj_Collection.py
@@ -121,15 +121,15 @@ def test_Collection_basics():
 def test_col_getB(test_input, expected):
     """testing some Collection stuff with getB"""
     src1 = magpy.magnet.Cuboid(
-        polarization=(1, 0, 1), dimension=(8, 4, 6), position=(0, 0, 0)
+        polarization=(1, 0, 1), dimension=(1, 1, 1), position=(0, 0, 0)
     )
     src2 = magpy.magnet.Cylinder(
-        polarization=(0, 1, 0), dimension=(8, 5), position=(-15, 0, 0)
+        polarization=(0, 1, 0), dimension=(1, 1), position=(-1, 0, 0)
     )
-    sens1 = magpy.Sensor(position=(0, 0, 6))
-    sens2 = magpy.Sensor(position=(0, 0, 6))
-    sens3 = magpy.Sensor(position=(0, 0, 6))
-    sens4 = magpy.Sensor(position=(0, 0, 6))
+    sens1 = magpy.Sensor(position=(0, 0, 1))
+    sens2 = magpy.Sensor(position=(0, 0, 1))
+    sens3 = magpy.Sensor(position=(0, 0, 1))
+    sens4 = magpy.Sensor(position=(0, 0, 1))
 
     sens_col = sens1 + sens2 + sens3 + sens4
     src_col = src1 + src2
diff --git a/tests/test_obj_Cuboid.py b/tests/test_obj_Cuboid.py
index ce8c7ef33..0cb55083b 100644
--- a/tests/test_obj_Cuboid.py
+++ b/tests/test_obj_Cuboid.py
@@ -4,9 +4,9 @@
 import numpy as np
 
 import magpylib as magpy
+from magpylib._src.fields.field_BH_cuboid import BHJM_magnet_cuboid
 from magpylib._src.obj_classes.class_Sensor import Sensor
 
-
 # # # """data generation for test_Cuboid()"""
 
 # # N = 100
@@ -99,12 +99,8 @@ def test_cuboid_object_vs_lib():
     pol = 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(
-        field="B", observers=pos, polarization=pol, dimension=dim
-    )
-    H0 = magpy.core.magnet_cuboid_field(
-        field="H", observers=pos, polarization=pol, dimension=dim
-    )
+    B0 = BHJM_magnet_cuboid(field="B", observers=pos, polarization=pol, dimension=dim)
+    H0 = BHJM_magnet_cuboid(field="H", observers=pos, polarization=pol, dimension=dim)
 
     src = magpy.magnet.Cuboid(polarization=pol[0], dimension=dim[0])
     B1 = src.getB(pos)
@@ -112,3 +108,52 @@ def test_cuboid_object_vs_lib():
 
     np.testing.assert_allclose(B0[0], B1)
     np.testing.assert_allclose(H0[0], H1)
+
+
+def test_getM():
+    """getM test"""
+    m0 = (0, 0, 0)
+    m1 = (10, 200, 3000)
+    cube = magpy.magnet.Cuboid(dimension=(2, 2, 2), magnetization=m1)
+    obs = [
+        (2, 2, 2),
+        (0, 0, 0),
+        (0.5, 0.5, 0.5),
+        (3, 0, 0),
+    ]
+    sens = magpy.Sensor(pixel=obs)
+
+    M1 = cube.getM(obs)
+    M2 = magpy.getM(cube, sens)
+    M3 = sens.getM(cube)
+    Mtest = np.array([m0, m1, m1, m0])
+
+    np.testing.assert_allclose(M1, Mtest)
+    np.testing.assert_allclose(M2, Mtest)
+    np.testing.assert_allclose(M3, Mtest)
+
+
+def test_getJ():
+    """getM test"""
+    j0 = (0, 0, 0)
+    j1 = (0.1, 0.2, 0.3)
+    cube = magpy.magnet.Cuboid(
+        dimension=(2, 2, 2),
+        polarization=j1,
+    )
+    obs = [
+        (-2, 2, -2),
+        (0, 0, 0),
+        (-0.5, -0.5, 0.5),
+        (-3, 0, 0),
+    ]
+    sens = magpy.Sensor(pixel=obs)
+
+    J1 = cube.getJ(obs)
+    J2 = magpy.getJ(cube, sens)
+    J3 = sens.getJ(cube)
+
+    Jtest = np.array([j0, j1, j1, j0])
+    np.testing.assert_allclose(J1, Jtest)
+    np.testing.assert_allclose(J2, Jtest)
+    np.testing.assert_allclose(J3, Jtest)
diff --git a/tests/test_obj_Cylinder.py b/tests/test_obj_Cylinder.py
index c7a9cef69..93adb6e6d 100644
--- a/tests/test_obj_Cylinder.py
+++ b/tests/test_obj_Cylinder.py
@@ -110,3 +110,53 @@ def test_Cylinder_getBH():
         assert np.allclose(H1, H2)
         assert np.allclose(H1, H3)
         assert np.allclose(H1, H0)
+
+
+def test_getM():
+    """getM test"""
+    m0 = (0, 0, 0)
+    m1 = (10, 200, 3000)
+    cyl = magpy.magnet.Cylinder(dimension=(2, 2), magnetization=m1)
+    obs = [
+        (2, 2, 2),
+        (0, 0, 0),
+        (0.5, 0.5, 0.5),
+        (3, 0, 0),
+    ]
+    sens = magpy.Sensor(pixel=obs)
+
+    M1 = cyl.getM(obs)
+    M2 = magpy.getM(cyl, sens)
+    M3 = sens.getM(cyl)
+
+    Mtest = np.array([m0, m1, m1, m0])
+
+    np.testing.assert_allclose(M1, Mtest)
+    np.testing.assert_allclose(M2, Mtest)
+    np.testing.assert_allclose(M3, Mtest)
+
+
+def test_getJ():
+    """getM test"""
+    j0 = (0, 0, 0)
+    j1 = (0.1, 0.2, 0.3)
+    cyl = magpy.magnet.Cylinder(
+        dimension=(2, 2),
+        polarization=j1,
+    )
+    obs = [
+        (-2, 2, -2),
+        (0, 0, 0),
+        (-0.5, -0.5, 0.5),
+        (-3, 0, 0),
+    ]
+    sens = magpy.Sensor(pixel=obs)
+
+    J1 = cyl.getJ(obs)
+    J2 = magpy.getJ(cyl, sens)
+    J3 = sens.getJ(cyl)
+    Jtest = np.array([j0, j1, j1, j0])
+
+    np.testing.assert_allclose(J1, Jtest)
+    np.testing.assert_allclose(J2, Jtest)
+    np.testing.assert_allclose(J3, Jtest)
diff --git a/tests/test_obj_Polyline.py b/tests/test_obj_Polyline.py
index f8224ae84..ce734ad3c 100644
--- a/tests/test_obj_Polyline.py
+++ b/tests/test_obj_Polyline.py
@@ -89,6 +89,7 @@ def test_line_position_bug():
 
 def test_discontinous_line():
     """test discontinuous line"""
+
     line_1 = magpy.current.Polyline(
         current=1,
         vertices=[
@@ -117,6 +118,7 @@ def test_discontinous_line():
 
     B1 = magpy.getB((line_1, line_2), (0, 0, 0), sumup=True)
     B2 = line_12.getB((0, 0, 0))
+
     np.testing.assert_allclose(B1, B2)
 
 
diff --git a/tests/test_obj_Sphere.py b/tests/test_obj_Sphere.py
index a72915650..3066e24a3 100644
--- a/tests/test_obj_Sphere.py
+++ b/tests/test_obj_Sphere.py
@@ -4,6 +4,7 @@
 import numpy as np
 
 import magpylib as magpy
+from magpylib._src.fields.field_BH_sphere import BHJM_magnet_sphere
 
 # # """data generation for test_Sphere()"""
 
@@ -94,9 +95,7 @@ def test_sphere_object_vs_lib():
     pol = np.array([(10, 20, 30)])
     dia = np.array([1])
     pos = np.array([(2, 2, 2)])
-    B1 = magpy.core.magnet_sphere_field(
-        field="B", observers=pos, polarization=pol, diameter=dia
-    )[0]
+    B1 = BHJM_magnet_sphere(field="B", observers=pos, polarization=pol, diameter=dia)[0]
 
     src = magpy.magnet.Sphere(polarization=pol[0], diameter=dia[0])
     B2 = src.getB(pos)
diff --git a/tests/test_obj_Tetrahedron.py b/tests/test_obj_Tetrahedron.py
index d24bf402e..250514994 100644
--- a/tests/test_obj_Tetrahedron.py
+++ b/tests/test_obj_Tetrahedron.py
@@ -61,3 +61,21 @@ def test_tetra_barycenter():
     vert = [(1, 1, -1), (1, 1, 1), (-1, 1, 1), (1, -1, 1)]
     tetra = magpy.magnet.Tetrahedron(polarization=pol, vertices=vert)
     np.testing.assert_allclose(tetra.barycenter, (0.5, 0.5, 0.5))
+
+
+def test_tetra_in_out():
+    """test inside and outside"""
+    pol = (0.111, 0.222, 0.333)
+    vert = [(-1, -1, -1), (0, 1, -1), (1, 0, -1), (0, 0, 1)]
+    tetra = magpy.magnet.Tetrahedron(polarization=pol, vertices=vert)
+
+    obs_in = [(0, 0, 0), (0.1, 0.1, 0.1), (0.2, 0.2, 0.2)]
+    obs_out = [(2, 0, 0), (2, 2, 2), (0.2, 2, 0.2)]
+
+    Bauto = magpy.getB(tetra, obs_in)
+    Bin = magpy.getB(tetra, obs_in, in_out="inside")
+    np.testing.assert_allclose(Bauto, Bin)
+
+    Bauto = magpy.getB(tetra, obs_out)
+    Bout = magpy.getB(tetra, obs_out, in_out="outside")
+    np.testing.assert_allclose(Bauto, Bout)
diff --git a/tests/test_obj_Triangle.py b/tests/test_obj_Triangle.py
index 7cecd60d3..8156d4007 100644
--- a/tests/test_obj_Triangle.py
+++ b/tests/test_obj_Triangle.py
@@ -2,6 +2,7 @@
 
 import magpylib as magpy
 from magpylib._src.exceptions import MagpylibMissingInput
+from magpylib._src.fields.field_BH_triangle import BHJM_triangle
 
 
 def test_Triangle_repr():
@@ -46,9 +47,7 @@ def test_triangle_input3():
             [(6, 0, 0), (10, 0, 0), (0, 10, 0)],
         ]
     )
-    b = magpy.core.triangle_field(
-        field="B", observers=obs, polarization=pol, vertices=vert
-    )
+    b = BHJM_triangle(field="B", observers=obs, polarization=pol, vertices=vert)
     b = np.sum(b, axis=0)
 
     tri1 = magpy.misc.Triangle(polarization=pol[0], vertices=vert[0])
diff --git a/tests/test_obj_TriangularMesh.py b/tests/test_obj_TriangularMesh.py
index bd59b0bbc..a78d4894c 100644
--- a/tests/test_obj_TriangularMesh.py
+++ b/tests/test_obj_TriangularMesh.py
@@ -9,9 +9,9 @@
 
 import magpylib as magpy
 from magpylib._src.exceptions import MagpylibBadUserInput
+from magpylib._src.fields.field_BH_triangularmesh import BHJM_magnet_trimesh
 from magpylib._src.fields.field_BH_triangularmesh import fix_trimesh_orientation
 from magpylib._src.fields.field_BH_triangularmesh import lines_end_in_trimesh
-from magpylib._src.fields.field_BH_triangularmesh import magnet_trimesh_field
 
 
 def test_TriangularMesh_repr():
@@ -141,7 +141,7 @@ def test_magnet_trimesh_func():
     pts_inside = np.array([[0, 0, 1]])
     B0 = cube.getB(pts_inside)
     B1 = tmesh_cube.getB(pts_inside)
-    B2 = magnet_trimesh_field(
+    B2 = BHJM_magnet_trimesh(
         field="B",
         observers=pts_inside,
         polarization=np.array([pol]),