diff --git a/lectures-wegewijs/steering.py b/lectures-wegewijs/steering.py
new file mode 100644
index 0000000000000000000000000000000000000000..e77ca55f8c1ce6c9e1e61d99e2410d8ddf70e6b5
--- /dev/null
+++ b/lectures-wegewijs/steering.py
@@ -0,0 +1,206 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+from  matplotlib.colors import LinearSegmentedColormap
+
+
+
+### Vectors
+
+def M1non(c,l):
+    s0 = np.sqrt(l); s1 = np.sqrt(1-s0**2)
+    v0=np.cos(c); v1=np.sin(c)
+    return np.transpose(np.array([[s0*v0,s1*v1],[-s0*v1,s1*v0]]))
+
+def M1(c,l):
+    s0 = np.sqrt(l); s1 = np.sqrt(1-s0**2)
+    x = M1non(c,l)
+    return x/np.linalg.norm(x,ord=2,axis=0)
+
+def M2(c):
+    v0=np.cos(c); v1=np.sin(c);
+    return np.transpose(np.array([[v0,v1],[-v1,v0]]))
+
+
+
+class Steering:
+    # define red-green colormap
+    colorlist = ["red","red","limegreen"]
+    segmentlist = [0.,0.4, 1.]
+    cmap=LinearSegmentedColormap.from_list('rg',
+        list(zip(segmentlist,colorlist)), N=20)
+    arrow_params = dict(
+            head_width=0.05,
+            head_length=0.05,
+            length_includes_head=True
+            )
+    alice_arrow_params = dict(width=0.01, color='black', **arrow_params)
+    bob_arrow_params = dict(width=0.01, color=(0.,0.,1.), **arrow_params)
+    bob_mirror_arrow_params = dict(width=0.01, color=(0.2,0.2,1.,0.3), **arrow_params)
+    marginal_arrow_params = dict(width=0.001, color='lightgray',linestyle='--', **arrow_params)
+
+    def __init__(self, resolution=200):
+        self.fig = plt.figure(figsize=plt.figaspect(0.5))
+        self.fig.subplots_adjust(left=0.15, bottom=0.25) 
+
+        self.ax1 = self.fig.add_subplot(1, 2, 1)
+        self.ax2 = self.fig.add_subplot(1, 2, 2)
+
+        self.ax1.axis('equal')  # set the axes to the same scale
+        self.ax2.axis('equal')
+        self.ax1.plot(0, 0, 'ok') # plot a black point at the origin
+        self.ax2.plot(0, 0, 'ok') 
+        self.arrows_alice = []
+        self.arrows_bob = []
+        m = 0.15 # margin
+        self.ax1.set_xlim([-1.0-m,1.0+m])
+        self.ax1.set_ylim([-1.0-m,1.0+m])
+        self.ax2.set_xlim([-1.0-m,1.0+m])
+        self.ax2.set_ylim([-1.0-m,1.0+m])
+
+        self.drag_pointer = 0
+
+        # Hide axes
+        self.ax1.xaxis.set_visible(False)
+        self.ax1.yaxis.set_visible(False)
+        self.ax2.xaxis.set_visible(False)
+        self.ax2.yaxis.set_visible(False)
+
+        self.ax1.set_title('Alice: Ensemble state vectors\n'
+                      r'$\vert \psi_{Ab} \rangle = \langle b_B \vert \psi_{AB} \rangle / |\langle b_B \vert \psi_{AB} \rangle|$')
+        self.ax2.set_title('Bob: Observation vectors\n'
+                      r'$\vert b_{B} \rangle$',)
+        self.ax1.arrow(-1,0,2,0,**self.marginal_arrow_params,zorder=-2)
+        self.ax1.arrow(0,-1,0,2,**self.marginal_arrow_params,zorder=-2)
+        self.ax2.arrow(-1,0,2,0,**self.marginal_arrow_params,zorder=-2)
+        self.ax2.arrow(0,-1,0,2,**self.marginal_arrow_params,zorder=-2)
+
+        self.ax1.text(0.85,-0.2, 'eigen-\n' r'vectors $\rho_A$',
+                      horizontalalignment='center', verticalalignment='center',
+                      color = 'grey')
+        self.ax2.text(0.85,-0.2, 'eigen-\n' r'vectors $\rho_B$',
+                      horizontalalignment='center', verticalalignment='center',
+                      color = 'grey')
+        self.ax2.text(-0.75,0.925, 'Grab vectors\n to rotate',
+                      horizontalalignment='center', verticalalignment='center',
+                      color = 'grey')
+
+        self.angles = np.linspace(0, 2*np.pi, resolution)
+
+
+
+        #### Slider ###
+        # Define initial values.
+        l_min = 1e-5   # avoid ugly drawings for l=0 ...
+        c = 0*np.pi
+        l = 0.5
+
+        ax_slider_c =  plt.axes([0.3, 0.1, 0.55, 0.03])
+        self.slider_c = Slider(
+                ax_slider_c,
+                'Eigenbasis $\\rho_B$: angle $\\theta$ $\\in$ [-$\\pi$, 3$\\pi$]',
+                -1*np.pi,
+                3*np.pi,
+                valinit=c
+                )
+        ax_slider_l =  plt.axes([0.3, 0.05, 0.55, 0.03])
+        self.slider_l = Slider(
+                ax_slider_l,
+                'Eigenvalue $\\rho_B$: $\lambda_0$ $\\in$ [0,1]',
+                l_min,
+                1-l_min,
+                valinit=l
+                )
+
+        self.slider_c.on_changed(self.update)
+        self.slider_l.on_changed(self.update)
+
+        # Dragging the arrows requires handling mouse input events:
+        self.fig.canvas.mpl_connect('button_press_event', self.on_press)
+        self.fig.canvas.mpl_connect('button_release_event', self.on_release)
+        self.fig.canvas.mpl_connect('motion_notify_event', self.drag)
+
+    def on_press(self, event):
+        'Mouse button press event: check if this should start dragging an arrow'
+        if event.inaxes != self.ax2:
+            return
+        if self.arrows_bob[0].contains(event, 500)[0] or self.arrows_bob[1].contains(event, 500)[0]:
+            vec1 = (event.xdata - self.arrows_bob[0].xy[0,0], event.ydata - self.arrows_bob[0].xy[0,1])
+            vec2 = (event.xdata - self.arrows_bob[1].xy[0,0], event.ydata - self.arrows_bob[1].xy[0,1])
+            if (vec1[0]**2 + vec1[1]**2 < vec2[0]**2 + vec2[1]**2):
+                self.drag_pointer = 1
+            else:
+                self.drag_pointer = 2
+
+
+    def on_release(self, event):
+        'Stop dragging arrow'
+        self.drag_pointer = 0
+
+    def drag(self, event):
+        '''Drag one of Bob's arrows'''
+        if event.inaxes != self.ax2:
+            return
+        if self.drag_pointer == 1:
+            self.slider_c.set_val(np.arctan2(event.ydata, event.xdata))
+        elif self.drag_pointer == 2:
+            self.slider_c.set_val(np.arctan2(-event.xdata, event.ydata))
+
+
+    def run(self):
+        'Start the animation.'
+        self.update()
+        plt.show()
+
+
+    def steer(self, c, l):
+        # Alice (left panel)
+        for arr in self.arrows_alice:
+            arr.remove()
+        for arr in self.arrows_bob:
+            arr.remove()
+        self.arrows_alice.clear()
+        self.arrows_bob.clear()
+        x = np.sqrt(l)*np.cos(self.angles)
+        y = np.sqrt(1-l)*np.sin(self.angles)
+        try:
+            self.ellipse.set_data(x, y)
+        except AttributeError:
+            self.ellipse, = self.ax1.plot(x, y, color='lightgray', zorder=-1)
+        M = M2(c)
+        for i in range(0, 2):
+            self.arrows_alice.append(
+                    self.ax1.arrow(0, 0, M[0,i], M[1,i], **self.bob_mirror_arrow_params)
+                    )
+        M = M1(c,l)
+        for i in range(0, 2):
+            self.arrows_alice.append(
+                    self.ax1.arrow(0, 0, M[0,i], M[1,i], **self.alice_arrow_params)
+                    )
+        M = M1non(c,l)
+        for i in range(0, 2):
+            self.arrows_alice.append(
+                    self.ax1.arrow(0, 0, M[0,i], M[1,i], **self.arrow_params, width=0.05,
+                      color = self.cmap( np.sqrt(M[0,i]**2+M[1,i]**2) ) )
+                    )
+        # Bob (right panel)    
+        M = M2(c)
+        for i in range(0, 2):
+            self.arrows_bob.append(
+                    self.ax2.arrow(0, 0, M[0,i], M[1,i], **self.bob_arrow_params)
+                    )
+
+
+    def update(self, *args): 
+        'Read values from the sliders and update the figure'
+        c = self.slider_c.val
+        l = self.slider_l.val
+        self.steer(c, l)
+        self.fig.canvas.draw_idle()
+
+
+
+if __name__ == '__main__':
+    Steering().run()
diff --git a/lectures-wegewijs/steering_bloch.py b/lectures-wegewijs/steering_bloch.py
new file mode 100644
index 0000000000000000000000000000000000000000..784914c05c0dd70e289d3fa0a8efa8031155b2ad
--- /dev/null
+++ b/lectures-wegewijs/steering_bloch.py
@@ -0,0 +1,274 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+from matplotlib.colors import LinearSegmentedColormap
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib.patches import FancyArrowPatch
+from mpl_toolkits.mplot3d import proj3d
+
+class Arrow3D(FancyArrowPatch):
+    # https://stackoverflow.com/questions/42281966/how-to-plot-vectors-in-python-using-matplotlib
+    def __init__(self, xs, ys, zs, *args, **kwargs):
+        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
+        self._verts3d = xs, ys, zs
+    def draw(self, renderer):
+        xs3d, ys3d, zs3d = self._verts3d
+        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
+        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
+        FancyArrowPatch.draw(self, renderer)
+
+# Bloch vectors
+def M1non(c,t,l):
+    'Probability-normalized Bloch vectors of A-ensemble'
+    s0 = np.sqrt(l); s1 = np.sqrt(1-s0**2)
+    xB=np.cos(c)*np.sin(t); yB=np.sin(c)*np.sin(t); zB=np.cos(t)
+    xA1=          s0*s1*xB;     xA2=        s0*s1*(-xB)
+    yA1=        - s0*s1*yB;     yA2=       -s0*s1*(-yB)
+    zA1= ((s0**2-s1**2)+zB)/2;  zA2= ((s0**2-s1**2)-zB)/2
+    return np.transpose(np.array([[xA1,yA1,zA1],[xA2,yA2,zA2]]))
+
+def M1(c,t,l):
+    'Normalized Bloch vectors of A-ensemble'
+    x = M1non(c,t,l)
+    return x/np.linalg.norm(x,ord=2,axis=0)
+
+def M2(c,t):
+    'Normalized Bloch vectors of B-measurement'
+    xB=np.cos(c)*np.sin(t); yB=np.sin(c)*np.sin(t); zB=np.cos(t)
+    return np.transpose(np.array([[ xB, yB, zB],
+                                  [-xB,-yB,-zB]]))
+# Probabilities
+def prob(c,t,l):
+    'Probabilities of A-ensemble'
+    s0 = np.sqrt(l); s1 = np.sqrt(1-s0**2)
+    zB=np.cos(t)
+    return np.transpose(np.array([ (1 + (s0**2-s1**2)*  zB  )/2,
+                                   (1 + (s0**2-s1**2)*(-zB) )/2 ]))
+
+#########################
+
+
+class Steering:
+    # Red-Green colormap
+    colorlist = ["red","red","limegreen"]
+    segmentlist = [0.,0.2, 1.]
+    cmap=LinearSegmentedColormap.from_list('rg', list(zip(segmentlist,colorlist)), N=20)
+
+    # Bloch sphere/circle parameters
+    # stride = key performance parameter: < 6 is slow,  > 6 gives ugly spheres
+    stride=6
+    sphere_params1 = dict(color='blue',linewidth=0.3,alpha=0.01,rstride=stride,cstride=stride)
+    sphere_params2 = sphere_params1
+    ellipsoid_params1 = dict(color='yellow',linewidth=0.1, alpha=0.05,rstride=stride,cstride=stride)
+    circle_params1 = dict(color='gray',linewidth=0.2, alpha=1.0)
+    circle_params2 = circle_params1
+
+    def __init__(self, resolution=200):
+        # Initialize figures
+        self.fig = plt.figure(figsize=plt.figaspect(0.5))
+        self.ax1 = self.fig.add_subplot(1, 2, 1,projection='3d') #, adjustable='box')
+        self.ax2 = self.fig.add_subplot(1, 2, 2,projection='3d')
+
+        # Initialize what needs NO update
+        self.fig.subplots_adjust(left=0.03,right=1-0.05,bottom=0.15,top=1-0.03,
+                                 wspace=0.01)
+
+        self.ax1.xaxis.set_ticks([]); self.ax1.yaxis.set_ticks([]); self.ax1.zaxis.set_ticks([]);
+        self.ax2.xaxis.set_ticks([]); self.ax2.yaxis.set_ticks([]); self.ax2.zaxis.set_ticks([]);
+        self.ax1.set_xlim([-1,1]);self.ax1.set_ylim([-1,1]);self.ax1.set_zlim([-1,1]);
+        self.ax2.set_xlim([-1,1]);self.ax2.set_ylim([-1,1]);self.ax2.set_zlim([-1,1]);
+        
+        # Bloch spheres
+        self.u = np.linspace(0, 2 * np.pi, 100)
+        self.v = np.linspace(0, np.pi, 100)
+        self.xsphere = 1 * np.outer(np.cos(self.u), np.sin(self.v))
+        self.ysphere = 1 * np.outer(np.sin(self.u), np.sin(self.v))
+        self.zsphere = 1 * np.outer(np.ones(np.size(self.u)), np.cos(self.v))
+        self.ax1.plot_surface(self.xsphere, self.ysphere, self.zsphere,**self.sphere_params1, zorder=-1)
+        self.ax2.plot_surface(self.xsphere, self.ysphere, self.zsphere,**self.sphere_params2, zorder=-1)
+
+        # Circles
+        self.tline = np.linspace(0, 2*np.pi, 100)
+        self.cline = np.linspace(0, 2*np.pi, 100)
+
+        # Z-axis of Bob
+        self.ax2.add_artist(Arrow3D([0,0],[0,0],[-1,1],
+                                    mutation_scale=20,lw=2,arrowstyle="-",color='lightgrey',zorder=-1))
+
+        # Text
+        self.ax1.text2D(0.99, 0.99,
+                        'Alice: Ensemble state Bloch-vectors\n'
+                        r'for $\vert \psi_{Ab} \rangle = \langle b_B \vert \psi_{AB} \rangle / |\langle b_B \vert \psi_{AB} \rangle|$',
+                        horizontalalignment = 'right', verticalalignment = 'top',
+                        transform=self.ax1.transAxes)
+        self.ax2.text2D( 0.99, 0.99,
+                        'Bob: Observation Bloch-vectors\n'
+                         r'for $\vert b_{B} \rangle$',
+                        horizontalalignment = 'right', verticalalignment = 'top',
+                        transform=self.ax2.transAxes)
+        self.ax2.text2D(0.99, 0.01,
+                        r'Grab figure to rotate: mouse (dx,dy) = ($d\phi$,$d\theta$)',
+                        color = 'grey',
+                        horizontalalignment = 'right', verticalalignment = 'bottom',
+                        transform=self.ax2.transAxes)
+    
+        # Initialize what needs to be updated
+        self.arrows1 = []; self.arrows2 = [];
+        self.drag_pointer = 0
+
+        # Initialize interactive stuff
+
+        # Initial values of l = \lambda_0, c = phi, t=theta
+        l_min = 1e-5   # avoid ugly drawings for l=0 ...
+        self.l = 0.
+        l = 0.5
+        c = 0*np.pi
+        t = 0.0*np.pi
+
+        # Sliders: positions adjusted relative to left and right panel
+        xoffset=0.3; yoffset=0.05; yheight=0.03
+        xslider=self.ax1.get_position().x0 + xoffset;
+        xwidth =self.ax2.get_position().x1 - xslider
+        yslider=self.ax1.get_position().y0 - yoffset
+        ax_slider_c =  plt.axes([xslider, yslider,             xwidth, yheight])
+        ax_slider_t =  plt.axes([xslider, yslider-1.5*yheight, xwidth, yheight])
+        ax_slider_l =  plt.axes([xslider, yslider-3.0*yheight, xwidth, yheight])
+        self.slider_c = Slider(ax_slider_c,
+                               'Eigenbasis $\\rho_B$: angle $\\phi$ $\\in$ [-$\\pi$, 3$\\pi$]',
+                               -1*np.pi, 3*np.pi, valinit=c)
+        self.slider_t = Slider(ax_slider_t,
+                               'Eigenbasis $\\rho_B$: angle $\\theta$ $\\in$ [-$\\pi$, 3$\\pi$]',
+                               -0.5*np.pi, 1.5*np.pi, valinit=t)
+        self.slider_l = Slider(ax_slider_l,
+                               'Eigenvalue $\\rho_B$: $\lambda_0$ $\\in$ [0,1]',
+                               l_min, 1-l_min,valinit=l)
+        self.slider_c.on_changed(self.update)
+        self.slider_t.on_changed(self.update)
+        self.slider_l.on_changed(self.update)
+        # Dragging: requires handling mouse input events:
+        self.fig.canvas.mpl_connect('button_press_event', self.on_press)
+        self.fig.canvas.mpl_connect('button_release_event', self.on_release)
+        self.fig.canvas.mpl_connect('motion_notify_event', self.drag)
+        self.ax1.disable_mouse_rotation(); self.ax2.disable_mouse_rotation()
+
+    def update_plots(self, c,t, l):
+        for arr in self.arrows1:
+            arr.remove()
+        for arr in self.arrows2:
+            arr.remove()
+        self.arrows1.clear(); self.arrows2.clear();
+        # Circles: \theta
+        xline=np.cos(c)*np.sin(self.tline); yline=np.sin(c)*np.sin(self.tline); zline=np.cos(self.tline)
+        try:
+            self.tcircle1.set_data(xline,-yline); self.tcircle1.set_3d_properties(zline) # update x,y,z data
+            self.tcircle2.set_data(xline, yline); self.tcircle2.set_3d_properties(zline) # update x,y,z data
+        except AttributeError:
+            self.tcircle1, = self.ax1.plot3D(xline, -yline, zline, **self.circle_params1,zorder=-1)
+            self.tcircle2, = self.ax2.plot3D(xline,  yline, zline, **self.circle_params2,zorder=-1)
+        # Circles: \phi (perhaps not helpful)
+        xline=np.cos(self.cline)*np.sin(t); yline=np.sin(self.cline)*np.sin(t); zline=np.cos(t)
+        try:
+            self.ccircle2.set_data(xline, yline); self.ccircle2.set_3d_properties(zline) # update x,y,z data
+        except AttributeError:
+            #self.ccircle1, = self.ax1.plot3D(xline, -yline, zline, **self.circle_params1,zorder=-1)
+            self.ccircle2, = self.ax2.plot3D(xline,  yline, zline, **self.circle_params2,zorder=-1)
+        # Ellipsoid of probability-normalized ensemble vectors
+        # (update only if l=\lambda changed which is not mouse-driven)
+        if self.l != l:
+            self.l = l
+            x = np.sqrt(l*(1-l))* self.xsphere; y = np.sqrt(l*(1-l))* self.ysphere; z = (l-(1-l))/2+ (1/2)* self.zsphere
+            try:
+                self.ellipsoid1.remove()
+                self.ellipsoid1 = self.ax1.plot_surface(x, y, z,**self.ellipsoid_params1,zorder=-1)
+            except (AttributeError,IndexError):
+                self.ellipsoid1 = self.ax1.plot_surface(x, y, z,**self.ellipsoid_params1,zorder=-1)
+
+        # Right panel (B)
+        # B-measurement Bloch vectors (normalized)
+        M=M2(c,t)
+        arrows= [ "-|>", "-" ]
+        for i in range(0,1+1):
+            self.arrows2.append(
+                self.ax2.add_artist(Arrow3D([0,M[0,i]],[0,M[1,i]],[0,M[2,i]],
+                               mutation_scale=10,lw=2, arrowstyle=arrows[i] )) #"-|>"))
+                )
+
+        # Left panel (A)
+        # A-ensemble Bloch vectors (normalized)
+        M=M1(c,t,l)
+        for i in range(0,1+1):
+            self.arrows1.append(
+                self.ax1.add_artist(Arrow3D([0,M[0,i]],[0,M[1,i]],[0,M[2,i]],
+                                   mutation_scale=10,lw=2, arrowstyle="-|>",color='black'))
+                )
+        #
+        self.arrows1.append(
+            self.ax1.add_artist(Arrow3D([M[0,0],M[0,1]],[M[1,0],M[1,1]],[M[2,0],M[2,1]],
+                                        mutation_scale=10,lw=1, arrowstyle="-", color = 'blue' ))
+            )
+        # A-ensemble Bloch vectors (probability normalized)
+        M=M1non(c,t,l)
+        for i in range(0,1+1):
+            self.arrows1.append(
+                self.ax1.add_artist(Arrow3D([0,M[0,i]],[0,M[1,i]],[0,M[2,i]],
+                               mutation_scale=10,lw=2, arrowstyle="-|>"
+                               ,color = self.cmap( prob(c,t,l)[i] )
+                ))
+            )
+        # A-mixed state Bloch vector
+        self.arrows1.append(
+            self.ax1.add_artist(Arrow3D([0,M[0,0]+M[0,1]],[0,M[1,0]+M[1,1]],[0,M[2,0]+M[2,1]],
+                                    mutation_scale=10,lw=2, arrowstyle="-|>", color = 'blue' ))
+            )
+        # Parallelogram of mixture (works, but is perhaps confusing)
+        #self.arrows1.append(
+        #    self.ax1.add_artist(Arrow3D([M[0,0],M[0,0]+M[0,1]],[M[1,0],M[1,0]+M[1,1]],[M[2,0],M[2,0]+M[2,1]],
+        #                                mutation_scale=10,lw=1, arrowstyle="-"
+        #                                ,color = self.cmap( prob(c,t,l)[1] )  ))
+        #)
+        #self.arrows1.append(
+        #    self.ax1.add_artist(Arrow3D([M[0,1],M[0,0]+M[0,1]],[M[1,1],M[1,0]+M[1,1]],[M[2,1],M[2,0]+M[2,1]],
+        #                                mutation_scale=10,lw=1, arrowstyle="-"
+        #                                ,color = self.cmap( prob(c,t,l)[0] )  ))
+        #)
+        
+    # Interactive stuff
+    def on_press(self, event):
+        'Mouse button press event: check if this should start dragging'
+        if event.inaxes != self.ax2:
+            return
+        self.drag_pointer = 1
+
+    def on_release(self, event):
+        'Stop dragging'
+        self.drag_pointer = 0
+
+    def drag(self, event):
+        'Drag right panel (B)'
+        if self.drag_pointer == 1:
+            # Map mouse positions (x,y) in [0,1]x[0,1] to (\phi,\theta) angles
+            x,y = self.ax2.transLimits.transform((event.xdata,event.ydata))
+            # Make-shift parameterization ...
+            self.slider_c.set_val((x-0.5)*1.25*np.pi)
+            self.slider_t.set_val((0.75-y)*2*np.pi)
+
+    def run(self):
+        'Start the animation.'
+        self.update()
+        plt.show()
+        
+    def update(self, *args): 
+        'Read values from the sliders and update the plots'
+        c = self.slider_c.val
+        t = self.slider_t.val
+        l = self.slider_l.val
+        self.update_plots(c,t, l)
+        self.fig.canvas.draw_idle()
+
+
+
+
+if __name__ == '__main__':
+    Steering().run()