beam me up

BeamTutorial.ipynb

@@ -220,7 +220,8 @@
     "<hr>\n",
     "\n",
     "Alright, so let's see it in action!\n",
-    "\n"
+    "\n",
+    "...\n"
    ]
   },
   {

spheres/beams.py

@@ -1,29 +1,15 @@
 from spheres import *
 
 import numpy as np
-from numpy import sqrt, pi, exp
+from numpy import sqrt, pi, exp, abs, angle
 from math import factorial
-from cmath import phase
 from scipy.special import eval_genlaguerre
 
-import matplotlib.pyplot as plt
-
-from colorsys import hls_to_rgb
-def colorize(z):
-    n,m = z.shape
-    c = np.zeros((n,m,3))
-    c[np.isinf(z)] = (1.0, 1.0, 1.0)
-    c[np.isnan(z)] = (0.5, 0.5, 0.5)
-    idx = ~(np.isinf(z) + np.isnan(z))
-    A = (np.angle(z[idx]) + np.pi) / (2*np.pi)
-    A = (A + 0.5) % 1.0
-    B = 1.0 - 1.0/(1.0+abs(z[idx]))
-    c[idx] = [hls_to_rgb(a, b, 1) for a,b in zip(A,B)]
-    return c
-
-# N is an integer.
-# l runs from -N to N in steps of 2. 
-def laguerre_gauss_beam(N, l, coordinates="cartesian"):
+def laguerre_gauss_mode(N, l, coordinates="cartesian"):
+    """
+    N is an integer.
+    l runs from -N to N in steps of 2.
+    """
     def mode(r, phi, z): # cylindrical coordinates
         w0 = 1 # waist radius
         n = 1 # index of refraction
@@ -43,15 +29,14 @@ def mode(r, phi, z): # cylindrical coordinates
     elif coordinates == "cartesian":
         def cartesian(x, y, z):
             c = x+1j*y
-            r = abs(c)
-            phi = phase(c)
+            r, phi = abs(c), angle(c)
             return mode(r, phi, z)
         return cartesian
 
 def spin_beam(spin, coordinates="cartesian"):
     j = (spin.shape[0]-1)/2
     v = components(spin)
-    lg_basis = [laguerre_gauss_beam(int(2*j), int(2*m), coordinates=coordinates) for m in np.arange(-j, j+1)]
+    lg_basis = [laguerre_gauss_mode(int(2*j), int(2*m), coordinates=coordinates) for m in np.arange(-j, j+1)]
     if coordinates == "cartesian":
         def beam(x, y, z):
             return sum([v[int(m+j)]*lg_basis[int(m+j)](x, y, z) for m in np.arange(-j, j+1)])
@@ -61,10 +46,87 @@ def beam(r, phi, z):
             return sum([v[int(m+j)]*lg_basis[int(m+j)](r, phi, z) for m in np.arange(-j, j+1)])
         return beam
 
-def viz_beam(beam, size=5, n_samples=100):
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+from mpl_toolkits.mplot3d import Axes3D
+from colorsys import hls_to_rgb
+
+def colorize(z):
+    n,m = z.shape
+    c = np.zeros((n,m,3))
+    c[np.isinf(z)] = (1.0, 1.0, 1.0)
+    c[np.isnan(z)] = (0.5, 0.5, 0.5)
+    idx = ~(np.isinf(z) + np.isnan(z))
+    A = (np.angle(z[idx]) + np.pi) / (2*np.pi)
+    A = (A + 0.5) % 1.0
+    B = 1.0 - 1.0/(1.0+abs(z[idx]))
+    c[idx] = [hls_to_rgb(a, b, 1) for a,b in zip(A,B)]
+    return c
+
+def viz_beam(beam, size=3.5, n_samples=100):
+    x = np.linspace(-size, size, n_samples)
+    y = np.linspace(-size, size, n_samples)
+    X, Y = np.meshgrid(x, y)
+    plt.imshow(colorize(beam(X, Y, np.zeros(X.shape))), interpolation="none", extent=(-size,size,-size,size))
+    plt.show()
+
+def viz_spin_beam(spin, size=3.5, n_samples=100):
+    stars = spin_xyz(spin)
+    beam = spin_beam(spin)
+    fig = plt.figure(figsize=plt.figaspect(0.5))
+
+    bloch_ax = fig.add_subplot(1, 2, 1, projection='3d')
+    sphere = qt.Bloch(fig=fig, axes=bloch_ax)
+    sphere.point_size=[300]*(spin.shape[0]-1)
+    sphere.add_points(stars.T)
+    sphere.add_vectors(stars)
+    sphere.make_sphere()
+
+    beam_ax = fig.add_subplot(1, 2, 2)
     x = np.linspace(-size, size, n_samples)
     y = np.linspace(-size, size, n_samples)
     X, Y = np.meshgrid(x, y)
-    Z = np.array([[beam(x_, y_, 0) for y_ in y] for x_ in x])
-    plt.imshow(colorize(Z), interpolation="none", extent=(-size,size,-size,size))
-    plt.show()
+    beam_ax.imshow(colorize(beam(X, Y, np.zeros(X.shape))), interpolation="none", extent=(-size,size,-size,size))
+
+    plt.show()
+
+def animate_spin_beam(spin, H, dt=0.1, T=2*np.pi, size=3.5, n_samples=100, filename=None, fps=20):
+    fig = plt.figure(figsize=plt.figaspect(0.5))
+
+    bloch_ax = fig.add_subplot(1, 2, 1, projection='3d')
+    sphere = qt.Bloch(fig=fig, axes=bloch_ax)
+    sphere.point_size=[300]*(spin.shape[0]-1)
+    sphere.make_sphere()
+
+    beam_ax = fig.add_subplot(1, 2, 2)
+    x = np.linspace(-size, size, n_samples)
+    y = np.linspace(-size, size, n_samples)
+    X, Y = np.meshgrid(x, y)
+    Z = np.zeros(X.shape)
+
+    U = (-1j*H*dt).expm()
+    sphere_history = []
+    beam_history = []
+    steps = int(T/dt)
+    for t in range(steps):
+        sphere_history.append(spin_xyz(spin))
+        beam = spin_beam(spin)
+        beam_history.append(beam(X, Y, Z))
+        spin = U*spin
+
+    sphere.make_sphere()
+    im = beam_ax.imshow(colorize(beam_history[0]), interpolation="none", extent=(-size,size,-size,size))
+
+    def animate(t):
+        sphere.clear()
+        sphere.add_points(sphere_history[t].T)
+        sphere.add_vectors(sphere_history[t])
+        sphere.make_sphere()
+
+        im.set_array(colorize(beam_history[t])) 
+        return [bloch_ax, im]
+
+    ani = animation.FuncAnimation(fig, animate, range(steps), repeat=False)
+    if filename:
+        ani.save(filename, fps=fps)
+    return ani

spheres/polyhedra.py

@@ -1,9 +1,9 @@
+import sys
 import random
 from math import pi, asin, atan2, cos, sin, sqrt
-import sys
 
 # https://www.chiark.greenend.org.uk/~sgtatham/polyhedra/
-def make_poly_points(n, verbose=False):
+def equidistribute_points(n, verbose=False):
     points = []
     for i in range(n):
         # Invent a randomly distributed point.