diff --git a/examples/tutorials/plot_custom_emissivity.py b/examples/tutorials/plot_custom_emissivity.py
new file mode 100644
index 000000000..349aa3d45
--- /dev/null
+++ b/examples/tutorials/plot_custom_emissivity.py
@@ -0,0 +1,91 @@
+"""
+Computing a camera image with custom emissivity
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This tutorial defines an emissivity that varies in space and computes the signal
+received by a camera using this emissivity.
+"""
+
+###############################################################################
+# We start by loading a built-in `tofu` configuration and define a 2D camera.
+
+import matplotlib.pyplot as plt
+import numpy as np
+import tofu as tf
+
+configB2 = tf.geom.utils.create_config("B2")
+
+cam2d = tf.geom.utils.create_CamLOS2D(
+    config=configB2,
+    P=[3.4, 0, 0],
+    N12=100,
+    F=0.1,
+    D12=0.1,
+    angs=[np.pi, np.pi/6, 0],
+    Name="",
+    Exp="",
+    Diag="",
+)
+
+###############################################################################
+# Now, we define an emissivity function that depends on r and z coordinates.
+# We can plot its profile in the (0, X, Z) plane.
+
+def emissivity(pts, t=None, vect=None):
+    """Custom emissivity as a function of geometry.
+
+    :param pts: ndarray of shape (3, n_points) (each column is a xyz coordinate)
+    :param t: optional, time parameter to add a time dependency to the
+        emissivity function
+    :param vect:
+    :return:
+        - emissivity -- 2D array holding the emissivity for each point in the
+            input grid
+    """
+    r, z = np.hypot(pts[0, :], pts[1, :]), pts[2, :]
+    e = np.exp(-(r - 2.4) ** 2 / 0.2 ** 2 - z ** 2 / 0.4 ** 2)
+    if t is not None:
+        e = np.cos(np.atleast_1d(t))[:, None] * e[None, :]
+    else:
+        # as stated in documentation of calc_signal, e.ndim must be 2
+        e = np.reshape(e, (1, -1))
+    return e
+
+
+y = np.linspace(2, 3, num=90)
+z = np.linspace(-0.5, 0.5, num=100)
+Y, Z = np.meshgrid(y, z)
+X = np.zeros_like(Y)
+pts = np.c_[X.ravel(), Y.ravel(), Z.ravel()].T
+emissivity_vals = emissivity(pts)
+emissivity_vals = emissivity_vals.reshape(X.shape)
+
+
+def project_to_2D(xyz):
+    """Projection to (0, X, Z) plane."""
+    return xyz[0], xyz[2]
+
+
+fig, ax = plt.subplots()
+ax.pcolormesh(Y, Z, emissivity_vals)
+ax.set_xlabel('y')
+ax.set_ylabel('z')
+configB2.plot(lax=ax, proj='cross')
+cam_center, = ax.plot(*project_to_2D(cam2d._dgeom['pinhole']), '*', ms=20)
+ax.legend(handles=[cam_center], labels=['camera pinhole'], loc='upper right')
+
+###############################################################################
+# Finally, we compute an image using the 2D camera and this emissivity.
+# If we provide a time vector, the field will vary in a cosinusoidal fashion
+# (see above definition) across time.
+
+time_vector = np.linspace(0, 2 * np.pi, num=100)
+
+sig, units = cam2d.calc_signal(emissivity,
+                               res=0.01,
+                               reflections=False,
+                               newcalc=True,
+                               plot=False,
+                               t=time_vector)
+sig.plot(ntMax=1)
+plt.show()
diff --git a/tofu/data/_plot.py b/tofu/data/_plot.py
index 67a885884..d99cb0abe 100644
--- a/tofu/data/_plot.py
+++ b/tofu/data/_plot.py
@@ -374,7 +374,7 @@ def _DataCam12D_plot(lData, key=None, nchMax=_nchMax, ntMax=_ntMax,
     lt = [dd.t for dd in lData]
     nt = lData[0].nt
     if nt == 1:
-        Dt = [t[0]-0.001,t[0]+0.001]
+        Dt = [lt[0][0]-0.001, lt[0][0]+0.001]
     else:
         Dt = np.array([[np.nanmin(t), np.nanmax(t)] for t in lt])
         Dt = [np.min(Dt[:,0]), np.max(Dt[:,1])]
diff --git a/tofu/utils.py b/tofu/utils.py
index 38a54167a..f727a83ca 100644
--- a/tofu/utils.py
+++ b/tofu/utils.py
@@ -546,7 +546,7 @@ def _get_load_npzmat_dict(out, pfe, mode='npz', exclude_keys=[]):
                 if out[k].ndim == 1:
                     dout[k] = out[k].tolist()
                 else:
-                    dout[k] = np.squeeze(out[k],axis=0).tolist()
+                    dout[k] = np.atleast_1d(np.squeeze(out[k],axis=0)).tolist()
                 if type(dout[k][0]) is str:
                     dout[k] = [kk.strip() for kk in dout[k]]
             else:
@@ -555,7 +555,7 @@ def _get_load_npzmat_dict(out, pfe, mode='npz', exclude_keys=[]):
                 dout[k] = tuple(dout[k])
         elif typ=='ndarray':
             if mode == 'mat':
-                dout[k] = np.squeeze(out[k])
+                dout[k] = np.atleast_1d(np.squeeze(out[k]))
                 if dout[k].shape == (0,0):
                     dout[k] = dout[k].ravel()
             else:
@@ -2484,11 +2484,15 @@ def func(li, val=val, n1=n1, n2=n2):
 
     else:
         assert type(val) is np.ndarray
-        val = val.squeeze()
+        # val = np.atleast_1d(val.squeeze())
         ndim = val.ndim
-        assert ndim >= len(lrids)
-        assert len(lrids) >= ninds
-        assert ndim >= ninds
+        c0 = ndim >= len(lrids) and len(lrids) >= ninds and ndim >= ninds
+        if not c0:
+            msg = "Wrong dimension / shape / references!\n"
+            msg += "    val.ndim  : {}\n".format(ndim)
+            msg += "    lrids     : {}\n".format(str(lrids))
+            msg += "    len(linds): {}\n".format(ninds)
+            raise Exception(msg)
 
         if ndim == ninds:
             if ndim == 1:
@@ -2566,13 +2570,17 @@ def func(val, ind0=None, refb=refb):
                         return np.nanargmin(np.abs(ref-val[1]))
 
             else:
-                refb = 0.5*(ref[1:]+ref[:-1])
-                if Type == 'x':
-                    def func(val, ind0=None, refb=refb):
-                        return np.digitize([val[0]], refb)[0]
+                if ref.size == 1:
+                    def func(val, ind0=None):
+                        return 0
                 else:
-                    def func(val, ind0=None, refb=refb):
-                        return np.digitize([val[1]], refb)[0]
+                    refb = 0.5*(ref[1:]+ref[:-1])
+                    if Type == 'x':
+                        def func(val, ind0=None, refb=refb):
+                            return np.digitize([val[0]], refb)[0]
+                    else:
+                        def func(val, ind0=None, refb=refb):
+                            return np.digitize([val[1]], refb)[0]
         elif indother is None:
             assert ref.ndim == 2
             if np.any(np.isnan(ref)):
diff --git a/tofu/version.py b/tofu/version.py
index 3a178e146..1fc557b2f 100644
--- a/tofu/version.py
+++ b/tofu/version.py
@@ -1,2 +1,2 @@
 # Do not edit, pipeline versioning governed by git tags!
-__version__ = '1.4.1-190-gd25622d'
+__version__ = '1.4.1-217-g754c8bb'