diff --git a/.github/workflows/python-testing-matrix.yml b/.github/workflows/python-testing-matrix.yml
index 047e69e..89a1609 100644
--- a/.github/workflows/python-testing-matrix.yml
+++ b/.github/workflows/python-testing-matrix.yml
@@ -17,7 +17,7 @@ jobs:
     strategy:
       fail-fast: true
       matrix:
-        os: [ubuntu-latest, windows-latest, macos-latest]
+        os: [ubuntu-latest, macos-latest]   # , windows-latest due to TcK install errors
         python-version: ["3.8", "3.9", "3.10", "3.11"]
 
     steps:
diff --git a/bsplines2d/_class01_select.py b/bsplines2d/_class01_select.py
index f3994d3..7d16cd0 100644
--- a/bsplines2d/_class01_select.py
+++ b/bsplines2d/_class01_select.py
@@ -133,8 +133,7 @@ def _select_ind(
                         f"Provided: {ind.shape}"
                     )
                     raise Exception(msg)
-                # make sure R varies first
-                ind_tup = ind.T.nonzero()[::-1]
+                ind_tup = ind.nonzero()
                 ind_bool = ind
 
             else:
@@ -197,18 +196,18 @@ def _select_ind(
                 ind_bool = ind_bool & cropiknots
 
                 # ind_tup is not 2d anymore
-                ind_tup = ind_bool.T.nonzero()[::-1]  # R varies first
+                ind_tup = ind_bool.nonzero()
                 # warnings.warn("ind is not 2d anymore!")
 
             elif ind_tup[0].shape == cropi.shape:
                 ind_bool = ind_bool & cropi
                 # ind_tup is not 2d anymore
-                ind_tup = ind_bool.T.nonzero()[::-1]  # R varies first
+                ind_tup = ind_bool.nonzero()
                 # warnings.warn("ind is not 2d anymore!")
 
             else:
                 ind_bool = ind_bool & cropi
-                ind_tup = ind_bool.T.nonzero()[::-1]
+                ind_tup = ind_bool.nonzero()
         else:
             ind_bool &= cropi
 
@@ -222,13 +221,14 @@ def _select_ind(
     elif returnas is tuple:
         out = ind_tup
     elif returnas == 'tuple-flat':
-        # make sure R is varying first
-        out = (ind_tup[0].T.ravel(), ind_tup[1].T.ravel())
+        # make sure R is varying first => no!
+        out = (ind_tup[0].ravel(), ind_tup[1].ravel())
     elif returnas is np.ndarray:
-        out = ind_tup[0] + ind_tup[1]*nR
+        # make sure R is varying first => no!
+        out = ind_tup[0]*nZ + ind_tup[1]
     elif returnas == 'array-flat':
-        # make sure R is varying first
-        out = (ind_tup[0] + ind_tup[1]*nR).T.ravel()
+        # make sure R is varying first => no!
+        out = (ind_tup[0]*nZ + ind_tup[1]).ravel()
     else:
         out = ind_bool
 
@@ -763,6 +763,7 @@ def _mesh2DRect_bsplines_knotscents(
             knots_per_bs_Z = knots_per_bs_Z[:, ind[1]]
 
         nknots = knots_per_bs_R.shape[0]
+        # TBC: reverse repeat and tile ?
         knots_per_bs_R = np.tile(knots_per_bs_R, (nknots, 1))
         knots_per_bs_Z = np.repeat(knots_per_bs_Z, nknots, axis=0)
 
@@ -779,6 +780,7 @@ def _mesh2DRect_bsplines_knotscents(
             cents_per_bs_Z = cents_per_bs_Z[:, ind[1]]
 
         ncents = cents_per_bs_R.shape[0]
+        # TBC: reverse repeat and tile ?
         cents_per_bs_R = np.tile(cents_per_bs_R, (ncents, 1))
         cents_per_bs_Z = np.repeat(cents_per_bs_Z, ncents, axis=0)
 
@@ -787,10 +789,11 @@ def _mesh2DRect_bsplines_knotscents(
 
     if return_knots is True and return_cents is True:
         out = (
-            (knots_per_bs_R, knots_per_bs_Z), (cents_per_bs_R, cents_per_bs_Z)
+            (knots_per_bs_R, knots_per_bs_Z),
+            (cents_per_bs_R, cents_per_bs_Z)
         )
     elif return_knots is True:
         out = (knots_per_bs_R, knots_per_bs_Z)
     else:
         out = (cents_per_bs_R, cents_per_bs_Z)
-    return out
\ No newline at end of file
+    return out
diff --git a/bsplines2d/_class02_bsplines_operators_rect.py b/bsplines2d/_class02_bsplines_operators_rect.py
index 1dd2516..0b7c1a6 100644
--- a/bsplines2d/_class02_bsplines_operators_rect.py
+++ b/bsplines2d/_class02_bsplines_operators_rect.py
@@ -7,7 +7,6 @@
 # Common
 import numpy as np
 import scipy.sparse as scpsp
-import datastock as ds
 
 
 from . import _utils_bsplines_operators as _operators
@@ -67,8 +66,8 @@ def get_mesh2dRect_operators(
     # prepare
 
     nx, ny = knotsx_per_bs.shape[1], knotsy_per_bs.shape[1]
-    kR = np.tile(knotsx_per_bs, ny)
-    kZ = np.repeat(knotsy_per_bs, nx, axis=1)
+    kR = np.repeat(knotsx_per_bs, ny, axis=1)
+    kZ = np.tile(knotsy_per_bs, nx)
     nbs = nx*ny
 
     if cropbs_flat is None:
@@ -219,7 +218,8 @@ def get_mesh2dRect_operators(
         for ir in range(nx):
             for iz in range(ny):
 
-                iflat = ir + iz*nx
+                # iflat = ir + iz*nx
+                iflat = iz + ir*ny
                 if cropbs_flat is not False and not cropbs_flat[iflat]:
                     continue
 
@@ -238,8 +238,8 @@ def get_mesh2dRect_operators(
                     if cropbs_flat is not False and not cropbs_flat[jflat]:
                         continue
 
-                    jr = jflat % nx
-                    jz = jflat // nx
+                    jr = jflat // ny
+                    jz = jflat % ny
 
                     # store (i, j) and (j, i) (symmetric matrix)
                     if jr >= ir:
@@ -277,17 +277,12 @@ def get_mesh2dRect_operators(
             )
 
             # surface elements
-            dZ = np.repeat(knotsy_mult[1:] - knotsy_mult[:-1], nx)
+            dZ = np.tile(knotsy_mult[1:] - knotsy_mult[:-1], nx)
             if geometry == 'linear':
-                dR = np.tile(
-                    knotsx_mult[1:] - knotsx_mult[:-1],
-                    ny,
-                )
+                dR = knotsx_mult[1:] - knotsx_mult[:-1]
             else:
-                dR = np.tile(
-                    0.5*(knotsx_mult[1:]**2 - knotsx_mult[:-1]**2),
-                    ny,
-                )
+                dR = 0.5*(knotsx_mult[1:]**2 - knotsx_mult[:-1]**2)
+            dR = np.repeat(dR, ny)
 
             dS = dR*dZ
             if cropbs_flat is not False:
@@ -328,7 +323,7 @@ def get_mesh2dRect_operators(
             for ir in range(nx):
                 for iz in range(ny):
 
-                    iflat = ir + iz*nx
+                    iflat = iz + ir*ny
                     if cropbs_flat is not False and not cropbs_flat[iflat]:
                         continue
 
@@ -348,8 +343,8 @@ def get_mesh2dRect_operators(
                         if cropbs_flat is not False and not cropbs_flat[jflat]:
                             continue
 
-                        jr = jflat % nx
-                        jz = jflat // nx
+                        jr = jflat // ny
+                        jz = jflat % ny
 
                         # store (i, j) and (j, i) (symmetric matrix)
                         if jr >= ir:
@@ -401,7 +396,7 @@ def get_mesh2dRect_operators(
         for ir in range(nx):
             for iz in range(ny):
 
-                iflat = ir + iz*nx
+                iflat = iz + ir*ny
                 if cropbs_flat is not False and not cropbs_flat[iflat]:
                     continue
 
@@ -422,8 +417,8 @@ def get_mesh2dRect_operators(
                     if cropbs_flat is not False and not cropbs_flat[jflat]:
                         continue
 
-                    jr = jflat % nx
-                    jz = jflat // nx
+                    jr = jflat // ny
+                    jz = jflat % ny
 
                     # store (i, j) and (j, i) (symmetric matrix)
                     if jr >= ir:
@@ -463,4 +458,4 @@ def get_mesh2dRect_operators(
 
         raise NotImplementedError("Integral D3N2 not implemented for deg=3!")
 
-    return opmat, operator, geometry
\ No newline at end of file
+    return opmat, operator, geometry
diff --git a/bsplines2d/_class02_bsplines_rect.py b/bsplines2d/_class02_bsplines_rect.py
index 0848657..3b32c09 100644
--- a/bsplines2d/_class02_bsplines_rect.py
+++ b/bsplines2d/_class02_bsplines_rect.py
@@ -34,7 +34,10 @@ def evaluate_spline(t, c, k, xp, nu, extrapolate, out):
         c1 = c.reshape(c.shape[:ndim] + (-1,))
         num_c_tr = c1.shape[-1]
         strides_c1 = [stride // c.dtype.itemsize for stride in c.strides]
-        indices_k1d = np.unravel_index(np.arange((k+1)**ndim), (k+1,)*ndim)[0][:, None]
+        indices_k1d = np.unravel_index(
+            np.arange((k+1)**ndim),
+            (k+1,)*ndim,
+        )[0][:, None]
         return scpinterp._bspl.evaluate_ndbspline(
             xp[:, None],
             t[None, :],
@@ -255,7 +258,7 @@ def ev_details(
         self,
         x0=None,
         x1=None,
-        #derivatives
+        # derivatives
         deriv=None,
         # others
         indbs_tf=None,
@@ -394,20 +397,20 @@ def _get_overlap(
 ):
     # nb of overlapping, inc. itself in 1d
     nbs0, nbs1 = shapebs
-    ind00 = np.tile(np.arange(0, nbs0), nbs1)
-    ind10 = np.repeat(np.arange(0, nbs1), nbs0)
+    ind00 = np.repeat(np.arange(0, nbs0), nbs1)
+    ind10 = np.tile(np.arange(0, nbs1), nbs0)
 
     # complete
     ntot = 2*deg + 1
 
-    add0= np.tile(np.arange(-deg, deg+1), ntot)
-    add1 = np.repeat(np.arange(-deg, deg+1), ntot)
+    add0 = np.repeat(np.arange(-deg, deg+1), ntot)
+    add1 = np.tile(np.arange(-deg, deg+1), ntot)
 
     inter0 = ind00[None, :] + add0[:, None]
     inter1 = ind10[None, :] + add1[:, None]
 
     # purge
-    inter = inter0 + inter1*nbs0
+    inter = inter1 + inter0*nbs1
     indneg = (
         (inter0 < 0) | (inter0 >= nbs0) | (inter1 < 0) | (inter1 >= nbs1)
     )
@@ -469,4 +472,4 @@ def get_bs_class(
         knots1=knots1,
         deg=deg,
         shapebs=shapebs,
-    )
\ No newline at end of file
+    )
diff --git a/bsplines2d/_class02_interpolate.py b/bsplines2d/_class02_interpolate.py
index 4c1d665..d90eec6 100644
--- a/bsplines2d/_class02_interpolate.py
+++ b/bsplines2d/_class02_interpolate.py
@@ -247,9 +247,9 @@ def interpolate(
 
             x0_str = isinstance(x0, str)
             # if ref_com is None:
-                # x0 = dout_temp[kd0[0]]['data']
-                # if len(kd0) > 1:
-                    # x1 = dout_temp[kd0[1]]['data']
+            # x0 = dout_temp[kd0[0]]['data']
+            # if len(kd0) > 1:
+            # x1 = dout_temp[kd0[1]]['data']
             # else:
             x0 = dout_temp[kd0[0]]['key']
             if len(kd0) > 1:
@@ -386,7 +386,6 @@ def interpolate(
         #     err.args = (msg,)
         #     raise err
 
-
     # ----------
     # store
 
@@ -565,7 +564,7 @@ def _check_keys(
 
                 if not (hasref and hasvect):
                     msg = (
-                        f"Provided ref_key[{ii}] not a valid ref or ref vector!\n"
+                        f"Provided ref_key[{ii}] invalid ref or ref vector!\n"
                         f"Provided: {rr}"
                     )
                     raise Exception(msg)
@@ -579,7 +578,10 @@ def _check_keys(
 
             lok_nobs = [
                 k0 for k0, v0 in coll.ddata.items()
-                if all([coll.ddata[rr]['ref'][0] in v0['ref'] for rr in ref_key])
+                if all([
+                    coll.ddata[rr]['ref'][0] in v0['ref']
+                    for rr in ref_key
+                ])
             ]
             lok_bs = []
 
@@ -588,7 +590,10 @@ def _check_keys(
 
     else:
         # binning only for non-bsplines or 1d bsplines
-        lok_nobs = [k0 for k0, v0 in coll.ddata.items() if k0 not in dbs.keys()]
+        lok_nobs = [
+            k0 for k0, v0 in coll.ddata.items()
+            if k0 not in dbs.keys()
+        ]
         lok_bs = [
             k0 for k0, v0 in coll.ddata.items()
             if (
@@ -757,17 +762,19 @@ def _check_params_bsplines(
             returnas = int
 
         # compute validated indbs array with appropriate form
-        indbs_tf = coll.select_ind(
+        indbs_tf_new = coll.select_ind(
             key=keybs,
             returnas=returnas,
             ind=indbs_tf,
             crop=crop,
         )
 
-        if isinstance(indbs_tf, tuple):
-            nbs = indbs_tf[0].size
+        if isinstance(indbs_tf_new, tuple):
+            nbs = indbs_tf_new[0].size
         else:
-            nbs = indbs_tf.size
+            nbs = indbs_tf_new.size
+    else:
+        indbs_tf_new = indbs_tf
 
     # ---------
     # submesh
@@ -780,7 +787,7 @@ def _check_params_bsplines(
     if coll.dobj[coll._which_mesh][keym]['subkey'] is None:
         submesh = False
 
-    return val_out, crop, cropbs, indbs_tf, nbs, submesh
+    return val_out, crop, cropbs, indbs_tf_new, nbs, submesh
 
 
 # ################################################################
@@ -861,8 +868,8 @@ def _submesh_addtemp(
 
     if isinstance(x0, str):
 
-       assert not any([rr is None for rr in dout_temp[kd0[0]]['ref']])
-       dr_add, lr_add = None, None
+        assert not any([rr is None for rr in dout_temp[kd0[0]]['ref']])
+        dr_add, lr_add = None, None
 
     else:
 
@@ -941,9 +948,9 @@ def _prepare_bsplines(
     # azone
 
     # azone = ds._generic_check._check_var(
-        # azone, 'azone',
-        # types=bool,
-        # default=True,
+    # azone, 'azone',
+    # types=bool,
+    # default=True,
     # )
 
     # -----------
@@ -1082,7 +1089,7 @@ def _interp(
 
                 # nan0
                 if nan0 is True:
-                      dout[k0]['data'][dout[k0]['data'] == 0.] = np.nan
+                    dout[k0]['data'][dout[k0]['data'] == 0.] = np.nan
 
         except Exception as err:
             raise err
@@ -1129,1006 +1136,3 @@ def _interp(
             ])
 
     return
-
-
-# #############################################################################
-# #############################################################################
-#                           Mesh2DRect - interp
-# #############################################################################
-
-
-# def _interp2d_check_RZ(
-    # R=None,
-    # Z=None,
-    # grid=None,
-# ):
-    # # R and Z provided
-    # if not isinstance(R, np.ndarray):
-        # try:
-            # R = np.atleast_1d(R).astype(float)
-        # except Exception as err:
-            # msg = "R must be convertible to np.arrays of floats"
-            # raise Exception(msg)
-    # if not isinstance(Z, np.ndarray):
-        # try:
-            # Z = np.atleast_1d(Z).astype(float)
-        # except Exception as err:
-            # msg = "Z must be convertible to np.arrays of floats"
-            # raise Exception(msg)
-
-    # # grid
-    # grid = ds._generic_check._check_var(
-        # grid, 'grid',
-        # default=R.shape != Z.shape,
-        # types=bool,
-    # )
-
-    # if grid is True and (R.ndim > 1 or Z.ndim > 1):
-        # msg = "If grid=True, R and Z must be 1d!"
-        # raise Exception(msg)
-    # elif grid is False and R.shape != Z.shape:
-        # msg = "If grid=False, R and Z must have the same shape!"
-        # raise Exception(msg)
-
-    # if grid is True:
-        # R = np.tile(R, Z.size)
-        # Z = np.repeat(Z, R.size)
-    # return R, Z, grid
-
-
-# def _interp2d_check(
-    # # ressources
-    # coll=None,
-    # # interpolation base, 1d or 2d
-    # key=None,
-    # # external coefs (optional)
-    # coefs=None,
-    # # interpolation points
-    # R=None,
-    # Z=None,
-    # radius=None,
-    # angle=None,
-    # grid=None,
-    # radius_vs_time=None,
-    # azone=None,
-    # # time: t or indt
-    # t=None,
-    # indt=None,
-    # indt_strict=None,
-    # # bsplines
-    # indbs=None,
-    # # parameters
-    # details=None,
-    # reshape=None,
-    # res=None,
-    # crop=None,
-    # nan0=None,
-    # val_out=None,
-    # imshow=None,
-    # return_params=None,
-    # store=None,
-    # inplace=None,
-# ):
-
-    # # -------------
-    # # keys
-
-    # dk = {
-        # k0: v0['bsplines']
-        # for k0, v0 in coll.ddata.items()
-        # if v0.get('bsplines') not in ['', None]
-        # and 'crop' not in k0
-    # }
-    # dk.update({kk: kk for kk in coll.dobj['bsplines'].keys()})
-    # if key is None and len(dk) == 1:
-        # key = list(dk.keys())[0]
-    # if key not in dk.keys():
-        # msg = (
-            # "Arg key must the key to a data referenced on a bsplines set\n"
-            # f"\t- available: {list(dk.keys())}\n"
-            # f"\t- provided: {key}\n"
-        # )
-        # raise Exception(msg)
-
-    # keybs = dk[key]
-    # keym = coll.dobj['bsplines'][keybs]['mesh']
-    # mtype = coll.dobj[coll._which_mesh][keym]['type']
-
-    # # -------------
-    # # details
-
-    # details = ds._generic_check._check_var(
-        # details, 'details',
-        # types=bool,
-        # default=False,
-    # )
-
-    # # -------------
-    # # crop
-
-    # crop = ds._generic_check._check_var(
-        # crop, 'crop',
-        # types=bool,
-        # default=mtype == 'rect',
-        # allowed=[False, True] if mtype == 'rect' else [False],
-    # )
-
-    # # -------------
-    # # nan0
-
-    # nan0 = ds._generic_check._check_var(
-        # nan0, 'nan0',
-        # types=bool,
-        # default=True,
-    # )
-
-    # # -------------
-    # # val_out
-
-    # val_out = ds._generic_check._check_var(
-        # val_out, 'val_out',
-        # default=np.nan,
-        # allowed=[False, np.nan, 0.]
-    # )
-
-    # # -----------
-    # # time
-
-    # # refbs and hastime
-    # refbs = coll.dobj['bsplines'][keybs]['ref']
-
-    # if mtype == 'polar':
-        # radius2d = coll.dobj[coll._which_mesh][keym]['radius2d']
-        # # take out key if bsplines
-        # lk = [kk for kk in [key, radius2d] if kk != keybs]
-
-        # #
-        # hastime, reft, keyt, t, dind = coll.get_time_common(
-            # keys=lk,
-            # t=t,
-            # indt=indt,
-            # ind_strict=indt_strict,
-        # )
-
-        # rad2d_hastime = radius2d in dind.keys()
-        # if key == keybs:
-            # kind = radius2d
-        # else:
-            # kind = key
-            # hastime = key in dind.keys()
-
-        # if hastime:
-
-            # indt = dind[kind].get('ind')
-
-            # # Special case: all times match
-            # if indt is not None:
-                # rtk = coll.get_time(key)[2]
-                # if indt.size == coll.dref[rtk]['size']:
-                    # if np.allclose(indt, np.arange(0, coll.dref[rtk]['size'])):
-                        # indt = None
-
-            # if indt is not None:
-                # indtu = np.unique(indt)
-                # indtr = np.array([indt == iu for iu in indtu])
-            # else:
-                # indtu, indtr = None, None
-        # else:
-            # indt, indtu, indtr = None, None, None
-
-    # elif key != keybs:
-        # # hastime, t, indit
-        # hastime, hasvect, reft, keyt, t, dind = coll.get_time(
-            # key=key,
-            # t=t,
-            # indt=indt,
-            # ind_strict=indt_strict,
-        # )
-        # if dind is None:
-            # indt, indtu, indtr = None, None, None
-        # else:
-            # indt, indtu, indtr = dind['ind'], dind['indu'], dind['indr']
-    # else:
-        # hastime, hasvect = False, False
-        # reft, keyt, indt, indtu, indtr = None, None, None, None, None
-
-    # # -----------
-    # # indbs
-
-    # if details is True:
-        # if mtype == 'rect':
-            # returnas = 'tuple-flat'
-        # elif mtype == 'tri':
-            # returnas = int
-        # elif mtype == 'polar':
-            # if len(coll.dobj['bsplines'][keybs]['shape']) == 2:
-                # returnas = 'array-flat'
-            # else:
-                # returnas = int
-
-        # # compute validated indbs array with appropriate form
-        # indbs_tf = coll.select_ind(
-            # key=keybs,
-            # returnas=returnas,
-            # ind=indbs,
-        # )
-    # else:
-        # indbs_tf = None
-
-    # # -----
-    # # store
-
-    # store = ds._generic_check._check_var(
-        # store, 'store',
-        # types=bool,
-        # default=False,
-        # allowed=[False] if details else [True, False],
-    # )
-
-    # # -----------
-    # # coordinates
-
-    # # (R, Z) vs (radius, angle)
-    # lc = [
-        # R is None and Z is None and radius is None,
-        # R is not None and Z is not None and store is False,
-        # (R is None and Z is None)
-        # and (radius is not None and mtype == 'polar') and store is False
-    # ]
-
-    # if not any(lc):
-        # msg = (
-            # "Please provide either:\n"
-            # "\t- R and Z: for any mesh type\n"
-            # "\t- radius (and angle): for polar mesh\n"
-            # "\t- None: for rect / tri mesh\n"
-        # )
-        # raise Exception(msg)
-
-    # # R, Z
-    # if lc[0]:
-
-        # if store is True:
-            # if imshow:
-                # msg = "Storing only works if imshow = False"
-                # raise Exception(msg)
-
-        # # no spec => sample mesh
-        # R, Z = coll.get_sample_mesh(
-            # key=keym,
-            # res=res,
-            # mode='abs',
-            # grid=True,
-            # R=R,
-            # Z=Z,
-            # imshow=imshow,
-        # )
-        # lc[1] = True
-
-    # if lc[1]:
-
-        # # check R, Z
-        # R, Z, grid = _interp2d_check_RZ(R=R, Z=Z, grid=grid)
-
-        # # special case if polar mesh => (radius, angle) from (R, Z)
-        # if mtype == 'polar':
-
-            # rad2d_indt = dind[radius2d].get('ind') if rad2d_hastime else None
-
-            # # compute radius2d at relevant times
-            # radius, _, _ = coll.interpolate_profile2d(
-                # # coordinates
-                # R=R,
-                # Z=Z,
-                # grid=False,
-                # # quantities
-                # key=radius2d,
-                # details=False,
-                # # time
-                # indt=rad2d_indt,
-            # )
-
-            # # compute angle2d at relevant times
-            # angle2d = coll.dobj[coll._which_mesh][keym]['angle2d']
-            # if angle2d is not None:
-                # angle, _, _ = coll.interpolate_profile2d(
-                    # # coordinates
-                    # R=R,
-                    # Z=Z,
-                    # grid=False,
-                    # # quantities
-                    # key=angle2d,
-                    # details=False,
-                    # # time
-                    # indt=rad2d_indt,
-                # )
-
-            # # simplify if not time-dependent
-            # if radius2d not in dind:
-                # assert radius.ndim == R.ndim
-                # if angle2d is not None:
-                    # assert angle.ndim == R.ndim
-
-            # # check consistency
-            # if rad2d_hastime != (radius.ndim == R.ndim + 1):
-                # msg = f"Inconsistency! {radius.shape}"
-                # raise Exception(msg)
-
-            # radius_vs_time = rad2d_hastime
-
-        # else:
-            # radius_vs_time = False
-
-    # else:
-        # radius_vs_time = ds._generic_check._check_var(
-            # radius_vs_time, 'radius_vs_time',
-            # types=bool,
-            # default=False,
-        # )
-
-    # # -------------
-    # # radius, angle
-
-    # if mtype == 'polar':
-
-        # # check same shape
-        # if not isinstance(radius, np.ndarray):
-            # radius = np.atleast_1d(radius)
-
-        # # angle vs angle2d
-        # angle2d = coll.dobj[coll._which_mesh][keym]['angle2d']
-        # if angle2d is not None and angle is None:
-            # msg = (
-                # f"Arg angle must be provided for bsplines {keybs}"
-            # )
-            # raise Exception(msg)
-
-        # # angle vs radius
-        # if angle is not None:
-            # if not isinstance(angle, np.ndarray):
-                # angle = np.atleast_1d(angle)
-
-            # if radius.shape != angle.shape:
-                # msg = (
-                    # "Args radius and angle must be np.ndarrays of same shape!\n"
-                    # f"\t- radius.shape: {radius.shape}\n"
-                    # f"\t- angle.shape: {angle.shape}\n"
-                # )
-                # raise Exception(msg)
-
-    # # -------------
-    # # coefs
-
-    # shapebs = coll.dobj['bsplines'][keybs]['shape']
-    # # 3 possible coefs shapes:
-    # #   - None (if details = True or key provided)
-    # #   - scalar or shapebs if details = False and key = keybs
-
-    # if details is True:
-        # coefs = None
-        # axis = None
-        # assert key == keybs, (key, keybs)
-
-    # else:
-        # if coefs is None:
-            # if key == keybs:
-                # if mtype == 'polar' and rad2d_hastime:
-                    # if rad2d_indt is None:
-                        # r2dnt = coll.dref[dind[radius2d]['ref']]['size']
-                    # else:
-                        # r2dnt = rad2d_indt.size
-                    # coefs = np.ones(tuple(np.r_[r2dnt, shapebs]), dtype=float)
-                    # axis = [ii + 1 for ii in range(len(shapebs))]
-                # else:
-                    # coefs = np.ones(shapebs, dtype=float)
-                    # axis = [ii for ii in range(len(shapebs))]
-            # else:
-                # coefs = coll.ddata[key]['data']
-                # axis = [
-                    # ii for ii in range(len(coll.ddata[key]['ref']))
-                    # if coll.ddata[key]['ref'][ii] in refbs
-                # ]
-
-        # elif key != keybs:
-            # msg = f"Arg coefs can only be provided if key = keybs!\n\t- key: {key}"
-            # raise Exception(msg)
-
-        # elif np.isscalar(coefs):
-            # coefs = np.full(shapebs, coefs)
-            # axis = [ii for ii in range(len(shapebs))]
-
-        # # consistency
-        # nshbs = len(shapebs)
-        # c0 = (
-            # coefs.shape[-nshbs:] == shapebs
-            # and (
-                # (hastime and coefs.ndim == nshbs + 1)           # pre-shaped
-                # or (not hastime and coefs.ndim == nshbs + 1)    # pre-shaped
-                # or (not hastime and coefs.ndim == nshbs)        # [None, ...]
-            # )
-            # and len(axis) == len(refbs)
-        # )
-        # if not c0:
-            # msg = (
-                # f"Inconsistency of '{key}' shape:\n"
-                # f"\t- shape: {coefs.shape}\n"
-                # f"\t- shapebs: {shapebs}\n"
-                # f"\t- hastime: {hastime}\n"
-                # f"\t- radius_vs_time: {radius_vs_time}\n"
-                # f"\t- refbs: {refbs}\n"
-                # f"\t- axis: {axis}"
-            # )
-            # raise Exception(msg)
-
-        # # Make sure coefs is time dependent
-        # if hastime:
-            # if indt is not None and (mtype == 'polar' or indtu is None):
-                # if coefs.shape[0] != indt.size:
-                    # # in case coefs is already provided with indt
-                    # coefs = coefs[indt, ...]
-
-            # if radius_vs_time and coefs.shape[0] != radius.shape[0]:
-                # msg = (
-                    # "Inconstistent coefs vs radius!\n"
-                    # f"\t- coefs.shape = {coefs.shape}\n"
-                    # f"\t- radius.shape = {radius.shape}\n"
-                # )
-                # raise Exception(msg)
-        # else:
-            # if radius_vs_time is True:
-                # sh = tuple([radius.shape[0]] + [1]*len(shapebs))
-                # coefs = np.tile(coefs, sh)
-                # axis = [aa + 1 for aa in axis]
-            # elif coefs.ndim == nshbs:
-                # coefs = coefs[None, ...]
-                # axis = [aa + 1 for aa in axis]
-
-    # # -------------
-    # # azone
-
-    # azone = ds._generic_check._check_var(
-        # azone, 'azone',
-        # types=bool,
-        # default=True,
-    # )
-
-    # # -------------
-    # # return_params
-
-    # return_params = ds._generic_check._check_var(
-        # return_params, 'return_params',
-        # types=bool,
-        # default=False,
-    # )
-
-    # # -------
-    # # inplace
-
-    # inplace = ds._generic_check._check_var(
-        # inplace, 'inplace',
-        # types=bool,
-        # default=store,
-    # )
-
-    # return (
-        # key, keybs,
-        # R, Z,
-        # radius, angle,
-        # coefs,
-        # axis,
-        # hastime,
-        # reft, keyt,
-        # shapebs,
-        # radius_vs_time,
-        # azone,
-        # indbs, indbs_tf,
-        # t, indt, indtu, indtr,
-        # details, crop,
-        # nan0, val_out,
-        # return_params,
-        # store, inplace,
-    # )
-
-
-# def interp2d(
-    # # ressources
-    # coll=None,
-    # # interpolation base, 1d or 2d
-    # key=None,
-    # # external coefs (instead of key, optional)
-    # coefs=None,
-    # # interpolation points
-    # R=None,
-    # Z=None,
-    # radius=None,
-    # angle=None,
-    # grid=None,
-    # radius_vs_time=None,        # if radius is provided, in case radius vs time
-    # azone=None,                 # if angle2d is interpolated, exclusion zone
-    # # time: t or indt
-    # t=None,
-    # indt=None,
-    # indt_strict=None,
-    # # bsplines
-    # indbs=None,
-    # # parameters
-    # details=None,
-    # reshape=None,
-    # res=None,
-    # crop=None,
-    # nan0=None,
-    # val_out=None,
-    # imshow=None,
-    # return_params=None,
-    # store=None,
-    # inplace=None,
-# ):
-
-    # # ---------------
-    # # check inputs
-
-    # (
-        # key, keybs,
-        # R, Z,
-        # radius, angle,
-        # coefs,
-        # axis,
-        # hastime,
-        # reft, keyt,
-        # shapebs,
-        # radius_vs_time,
-        # azone,
-        # indbs, indbs_tf,
-        # t, indt, indtu, indtr,
-        # details, crop,
-        # nan0, val_out,
-        # return_params,
-        # store, inplace,
-    # ) = _interp2d_check(
-        # # ressources
-        # coll=coll,
-        # # interpolation base, 1d or 2d
-        # key=key,
-        # # external coefs (optional)
-        # coefs=coefs,
-        # # interpolation points
-        # R=R,
-        # Z=Z,
-        # radius=radius,
-        # angle=angle,
-        # grid=grid,
-        # radius_vs_time=radius_vs_time,
-        # azone=azone,
-        # # time: t or indt
-        # t=t,
-        # indt=indt,
-        # indt_strict=indt_strict,
-        # # bsplines
-        # indbs=indbs,
-        # # parameters
-        # details=details,
-        # res=res,
-        # crop=crop,
-        # nan0=nan0,
-        # val_out=val_out,
-        # imshow=imshow,
-        # return_params=return_params,
-        # store=store,
-        # inplace=inplace,
-    # )
-    # keym = coll.dobj['bsplines'][keybs]['mesh']
-    # meshtype = coll.dobj['mesh'][keym]['type']
-
-    # # ----------------------
-    # # which function to call
-
-    # if details is True:
-        # fname = 'func_details'
-    # elif details is False:
-        # fname = 'func_sum'
-    # else:
-        # raise Exception("Unknown details!")
-
-    # # ---------------
-    # # cropping ?
-
-    # cropbs = coll.dobj['bsplines'][keybs]['crop']
-    # if cropbs not in [None, False]:
-        # cropbs = coll.ddata[cropbs]['data']
-        # nbs = cropbs.sum()
-    # else:
-        # nbs = np.prod(coll.dobj['bsplines'][keybs]['shape'])
-
-    # if indbs_tf is not None:
-        # if isinstance(indbs_tf, tuple):
-            # nbs = indbs_tf[0].size
-        # else:
-            # nbs = indbs_tf.size
-
-    # # -----------
-    # # Interpolate
-
-    # if meshtype in ['rect', 'tri']:
-
-        # # manage time
-        # if indtu is not None:
-            # val0 = np.full(tuple(np.r_[indt.size, R.shape]), np.nan)
-            # coefs = coefs[indtu, ...]
-
-        # val = coll.dobj['bsplines'][keybs][fname](
-            # R=R,
-            # Z=Z,
-            # coefs=coefs,
-            # axis=axis,
-            # crop=crop,
-            # cropbs=cropbs,
-            # indbs_tf=indbs_tf,
-            # val_out=val_out,
-        # )
-
-        # # manage time
-        # if indtu is not None:
-            # for ii, iu in enumerate(indtu):
-                # val0[indtr[ii], ...] = val[ii, ...]
-            # val = val0
-
-        # shape_pts = R.shape
-
-    # elif meshtype == 'polar':
-
-        # val = coll.dobj['bsplines'][keybs][fname](
-            # radius=radius,
-            # angle=angle,
-            # coefs=coefs,
-            # axis=axis,
-            # indbs_tf=indbs_tf,
-            # radius_vs_time=radius_vs_time,
-            # val_out=val_out,
-        # )
-
-        # shape_pts = radius.shape
-
-    # # ---------------
-    # # post-treatment for angle2d only (discontinuity)
-
-    # # if azone is True:
-        # # lkm0 = [
-            # # k0 for k0, v0 in coll.dobj[coll._which_mesh].items()
-            # # if key == v0.get('angle2d')
-        # # ]
-        # # if len(lkm0) > 0:
-            # # ind = angle2d_inzone(
-                # # coll=coll,
-                # # keym0=lkm0[0],
-                # # keya2d=key,
-                # # R=R,
-                # # Z=Z,
-                # # t=t,
-                # # indt=indt,
-            # # )
-            # # assert val.shape == ind.shape
-            # # val[ind] = np.pi
-
-    # # ---------------
-    # # post-treatment
-
-    # if nan0 is True:
-        # val[val == 0] = np.nan
-
-    # if not hastime and not radius_vs_time:
-        # c0 = (
-            # (
-                # details is False
-                # and val.shape == tuple(np.r_[1, shape_pts])
-            # )
-            # or (
-                # details is True
-                # and val.shape == tuple(np.r_[shape_pts, nbs])
-            # )
-        # )
-        # if not c0:
-            # import pdb; pdb.set_trace()     # DB
-            # pass
-        # if details is False:
-            # val = val[0, ...]
-            # reft = None
-
-    # # ------
-    # # store
-
-    # # ref
-    # ct = (
-        # (hastime or radius_vs_time)
-        # and (
-            # reft in [None, False]
-            # or (
-                # reft not in [None, False]
-                # and indt is not None
-                # and not (
-                    # indt.size == coll.dref[reft]['size']
-                    # or np.allclose(indt, np.arange(0, coll.dref[reft]['size']))
-                # )
-            # )
-        # )
-    # )
-    # if ct:
-        # reft = f'{key}-nt'
-
-    # # store
-    # if store is True:
-        # Ru = np.unique(R)
-        # Zu = np.unique(Z)
-        # nR, nZ = Ru.size, Zu.size
-
-        # knR, knZ = f'{key}-nR', f'{key}-nZ'
-        # kR, kZ = f'{key}-R', f'{key}-Z'
-
-        # if inplace is True:
-            # coll2 = coll
-        # else:
-            # coll2 = ds.DataStock()
-
-        # # add ref nR, nZ
-        # coll2.add_ref(key=knR, size=nR)
-        # coll2.add_ref(key=knZ, size=nZ)
-
-        # # add data Ru, Zu
-        # coll2.add_data(
-            # key=kR,
-            # data=Ru,
-            # ref=knR,
-            # dim='distance',
-            # name='R',
-            # units='m',
-        # )
-        # coll2.add_data(
-            # key=kZ,
-            # data=Zu,
-            # ref=knZ,
-            # dim='distance',
-            # name='Z',
-            # units='m',
-        # )
-
-        # # ref
-        # if hastime or radius_vs_time:
-            # if ct:
-                # coll2.add_ref(key=reft, size=t.size)
-                # coll2.add_data(
-                    # key=f'{key}-t',
-                    # data=t,
-                    # ref=reft,
-                    # dim='time',
-                    # units='s',
-                # )
-            # else:
-                # if reft not in coll2.dref.keys():
-                    # coll2.add_ref(key=reft, size=coll.dref[reft]['size'])
-                # if keyt is not None and keyt not in coll2.ddata.keys():
-                    # coll2.add_data(
-                        # key=keyt,
-                        # data=coll.ddata[keyt]['data'],
-                        # dim='time',
-                    # )
-
-            # ref = (reft, knR, knZ)
-        # else:
-            # ref = (knR, knZ)
-
-        # coll2.add_data(
-            # data=val,
-            # key=f'{key}-map',
-            # ref=ref,
-            # dim=coll.ddata[key]['dim'],
-            # quant=coll.ddata[key]['quant'],
-            # name=coll.ddata[key]['name'],
-            # units=coll.ddata[key]['units'],
-        # )
-
-    # else:
-
-        # ref = []
-        # c0 = (
-            # reft not in [None, False]
-            # and (hastime or radius_vs_time)
-            # and not (
-                # meshtype == 'polar'
-                # and R is None
-                # and coefs is None
-                # and radius_vs_time is False
-            # )
-        # )
-        # if c0:
-            # ref.append(reft)
-
-        # if meshtype in ['rect', 'tri']:
-            # for ii in range(R.ndim):
-                # ref.append(None)
-            # if grid is True:
-                # for ii in range(Z.ndim):
-                    # ref.append(None)
-        # else:
-            # for ii in range(radius.ndim):
-                # if radius_vs_time and ii == 0:
-                    # continue
-                # ref.append(None)
-            # if grid is True and angle is not None:
-                # for ii in range(angle.ndim):
-                    # ref.append(None)
-
-        # if details is True:
-            # refbs = coll.dobj['bsplines'][keybs]['ref_bs'][0]
-            # if crop is True:
-                # refbs = f"{refbs}-crop"
-            # ref.append(refbs)
-        # ref = tuple(ref)
-
-        # if ref[0] == 'emiss-nt':
-            # import pdb; pdb.set_trace()     # DB
-
-        # if len(ref) != val.ndim:
-            # msg = (
-                # "Mismatching ref vs val.shape:\n"
-                # f"\t- key = {key}\n"
-                # f"\t- keybs = {keybs}\n"
-                # f"\t- val.shape = {val.shape}\n"
-                # f"\t- ref = {ref}\n"
-                # f"\t- reft = {reft}\n"
-                # f"\t- hastime = {hastime}\n"
-                # f"\t- radius_vs_time = {radius_vs_time}\n"
-                # f"\t- details = {details}\n"
-                # f"\t- indbs_tf = {indbs_tf}\n"
-                # f"\t- key = {key}\n"
-                # f"\t- meshtype = {meshtype}\n"
-                # f"\t- grid = {grid}\n"
-            # )
-            # if coefs is not None:
-                # msg += f"\t- coefs.shape = {coefs.shape}\n"
-            # if R is not None:
-                # msg += (
-                    # f"\t- R.shape = {R.shape}\n"
-                    # f"\t- Z.shape = {Z.shape}\n"
-                # )
-            # if meshtype == 'polar':
-                # msg += f"\t- radius.shape = {radius.shape}\n"
-                # if angle is not None:
-                    # msg += f"\t- angle.shape = {angle.shape}\n"
-            # raise Exception(msg)
-
-    # # ------
-    # # return
-
-    # if store and inplace is False:
-        # return coll2
-    # else:
-        # if return_params is True:
-            # return val, t, ref, dparams
-        # else:
-            # return val, t, ref
-
-
-# #############################################################################
-# #############################################################################
-#                   2d points to 1d quantity interpolation
-# #############################################################################
-
-
-# def interp2dto1d(
-    # coll=None,
-    # key=None,
-    # R=None,
-    # Z=None,
-    # indbs=None,
-    # indt=None,
-    # grid=None,
-    # details=None,
-    # reshape=None,
-    # res=None,
-    # crop=None,
-    # nan0=None,
-    # imshow=None,
-    # return_params=None,
-# ):
-
-    # # ---------------
-    # # check inputs
-
-    # # TBD
-    # pass
-
-    # # ---------------
-    # # post-treatment
-
-    # if nan0 is True:
-        # val[val == 0] = np.nan
-
-    # # ------
-    # # return
-
-    # if return_params is True:
-        # return val, dparams
-    # else:
-        # return val
-
-
-# #############################################################################
-# #############################################################################
-#                           Mesh2DRect - operators
-# #############################################################################
-
-
-def get_bsplines_operator(
-    coll,
-    key=None,
-    operator=None,
-    geometry=None,
-    crop=None,
-    store=None,
-    returnas=None,
-    # specific to deg = 0
-    centered=None,
-    # to return gradR, gradZ, for D1N2 deg 0, for tomotok
-    returnas_element=None,
-):
-
-    # check inputs
-    lk = list(coll.dobj.get('bsplines', {}).keys())
-    key = ds._generic_check._check_var(
-        key, 'key',
-        types=str,
-        allowed=lk,
-    )
-
-    store = ds._generic_check._check_var(
-        store, 'store',
-        default=True,
-        types=bool,
-    )
-
-    returnas = ds._generic_check._check_var(
-        returnas, 'returnas',
-        default=store is False,
-        types=bool,
-    )
-
-    crop = ds._generic_check._check_var(
-        crop, 'crop',
-        default=True,
-        types=bool,
-    )
-
-    # cropbs
-    cropbs = coll.dobj['bsplines'][key]['crop']
-    keycropped = coll.dobj['bsplines'][key]['ref_bs'][0]
-    if cropbs not in [None, False] and crop is True:
-        cropbs_flat = coll.ddata[cropbs]['data'].ravel(order='F')
-        if coll.dobj['bsplines'][key]['deg'] == 0:
-            cropbs = coll.ddata[cropbs]['data']
-        keycropped = f"{keycropped}-crop"
-    else:
-        cropbs = False
-        cropbs_flat = False
-
-    # compute and return
-    (
-        opmat, operator, geometry, dim,
-    ) = coll.dobj['bsplines'][key]['class'].get_operator(
-        operator=operator,
-        geometry=geometry,
-        cropbs_flat=cropbs_flat,
-        # specific to deg=0
-        cropbs=cropbs,
-        centered=centered,
-        # to return gradR, gradZ, for D1N2 deg 0, for tomotok
-        returnas_element=returnas_element,
-    )
-
-    # cropping
-    if operator == 'D1':
-        ref = (keycropped, keycropped)
-    elif operator == 'D0N1':
-        ref = (keycropped,)
-    elif 'N2' in operator:
-        ref = (keycropped, keycropped)
-
-    return opmat, operator, geometry, dim, ref, crop, store, returnas, key
\ No newline at end of file
diff --git a/bsplines2d/_class02_operators.py b/bsplines2d/_class02_operators.py
index 1eb67af..35d2b6b 100644
--- a/bsplines2d/_class02_operators.py
+++ b/bsplines2d/_class02_operators.py
@@ -164,7 +164,7 @@ def _check(
     cropbs = coll.dobj['bsplines'][key]['crop']
     keycropped = coll.dobj['bsplines'][key]['ref_bs'][0]
     if cropbs not in [None, False] and crop is True:
-        cropbs_flat = coll.ddata[cropbs]['data'].ravel(order='F')
+        cropbs_flat = coll.ddata[cropbs]['data'].ravel()
         if coll.dobj['bsplines'][key]['deg'] == 0:
             cropbs = coll.ddata[cropbs]['data']
         keycropped = f"{keycropped}-crop"
@@ -289,7 +289,7 @@ def _dout(
     elif operator in ['D1N2']:
         kmat = [f'tMM{ii}' for ii in range(nnd)]
     elif operator in ['D2N2']:
-        lcomb = [] if nd == '1d' else [(0,1)]
+        lcomb = [] if nd == '1d' else [(0, 1)]
         kmat = (
             [f'tMM{ii}{ii}' for ii in range(nnd)]
             + [f'tMM{ii}{jj}' for ii, jj in lcomb]
@@ -434,7 +434,7 @@ def _units(u0=None, operator=None, geometry=None):
         if geometry == 'linear':
             units = u0
         else:
-            units = u0**2
+            units = u0**2 / asunits.Unit('rad')
 
     elif operator == 'D0N2':
         if geometry == 'linear':
@@ -525,6 +525,7 @@ def apply_operator(
         keybs=keybs,
         # for units
         operator=operator,
+        integ_op=integ_op,
     )
 
     # -----------------
@@ -533,20 +534,12 @@ def apply_operator(
 
     if operator == 'D0N1':
         ind = [-1]
-        if cropbs is None:
-            for k0 in key:
-                ddata[k0]['data'][...] = np.tensordot(
-                    integ_op['M']['data'],
-                    coll.ddata[k0]['data'],
-                    (ind, daxis[k0]['axis']),
-                )
-        else:
-            for k0 in key:
-                ddata[k0]['data'][...] = np.tensordot(
-                    integ_op['M']['data'],
-                    coll.ddata[k0]['data'][daxis[k0]['slice']],
-                    (ind, daxis[k0]['axis']),
-                )
+        for k0 in key:
+            ddata[k0]['data'][...] = np.tensordot(
+                integ_op['M']['data'],
+                coll.ddata[k0]['data'][daxis[k0]['slice']],
+                (ind, daxis[k0]['axis'][0]),
+            )
     else:
         raise NotImplementedError()
 
@@ -699,6 +692,7 @@ def _apply_operator_prepare(
     key=None,
     keybs=None,
     operator=None,
+    integ_op=None,
 ):
 
     # ----------
@@ -716,14 +710,6 @@ def _apply_operator_prepare(
 
     refbs = coll.dobj[wbs][keybs]['ref']
 
-    # ----------
-    # unitsbs
-
-    unitsbs = [
-        coll.ddata[k0]['units']
-        for k0 in coll.dobj[wbs][keybs]['apex']
-    ]
-
     # --------------
     # fill dict
 
@@ -753,13 +739,11 @@ def _apply_operator_prepare(
         # units
         units0 = asunits.Unit(coll.ddata[k0]['units'])
         if operator == 'D0N1':
-            units = units0
-            for uu in unitsbs:
-                units = units * uu
+            units = units0 * integ_op['M']['units']
         else:
             raise NotImplementedError()
 
-         # populate
+        # populate
         ddata[k0] = {
             'data': np.full(shape, np.nan),
             'ref': ref,
@@ -768,18 +752,18 @@ def _apply_operator_prepare(
 
         # slicing
         if cropbs is None:
-            sli = None
-        else:
-            sli = tuple([
-                cropbs if ii == axis[0]
-                else slice(None)
-                for ii in range(len(shape0))
-                if ii not in axisf[1:]
-            ])
+            shcrop = tuple([shape0[ii] for ii in axis])
+            cropbs = np.ones(shcrop, dtype=bool)
+        sli = tuple([
+            cropbs if ii == axis[0]
+            else slice(None)
+            for ii in range(len(shape0))
+            if ii not in axisf[1:]
+        ])
 
         daxis[k0] = {
             'slice': sli,
             'axis': axis,
         }
 
-    return ddata, daxis, cropbs
\ No newline at end of file
+    return ddata, daxis, cropbs
diff --git a/bsplines2d/_utils_bsplines_operators.py b/bsplines2d/_utils_bsplines_operators.py
index b2307f3..f111e18 100644
--- a/bsplines2d/_utils_bsplines_operators.py
+++ b/bsplines2d/_utils_bsplines_operators.py
@@ -12,14 +12,14 @@
 
 
 def _D0N1_Deg0(kn, geom='linear'):
-    if geom  == 'linear':
+    if geom == 'linear':
         return kn[1, :] - kn[0, :]
     else:
         return 0.5 * (kn[1, :]**2 - kn[0, :]**2)
 
 
 def _D0N1_Deg1(kn, geom='linear'):
-    if geom  == 'linear':
+    if geom == 'linear':
         return 0.5 * (kn[2, :] - kn[0, :])
     else:
         return (
@@ -865,4 +865,4 @@ def _D2N2_Deg2_2_linear(k1, k2, k3, k4):
 def _D2N2_Deg2_2_toroidal(k1, k2, k3, k4):
     intt = np.zeros((k1.size,))
     intt[:-2] += 2.*(k3 + k2)[:-2] / ((k4 - k2)*(k3 - k2)*(k3 - k1))[:-2]
-    return intt
\ No newline at end of file
+    return intt
diff --git a/pyproject.toml b/pyproject.toml
index 0dcc680..5fed50a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -41,7 +41,7 @@ keywords = [
 requires-python = ">=3.8"
 dependencies = [
     "contourpy",
-	'datastock>=0.0.54',
+    'datastock>=0.0.55',
 ]