import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
to_degs, to_rads = 180/pi, pi/180
root2, r3o2 = np.sqrt(2.), np.sqrt(3)/2.

nmax = 16

i = np.arange(-nmax, nmax+1)
I, J = np.meshgrid(i, i)
keep = (np.abs(I + J) <= nmax) * ~((I==0) * (J==0)) # [1:] avoids (0, 0)
I, J = [thing[keep].flatten() for thing in (I, J)]
u_dots = np.vstack((I + 0.5*J, r3o2 * J))

k, l = 3, 1

nmax = 4
i = np.arange(-nmax, nmax+1)
I, J = np.meshgrid(i, i)
keep = (np.abs(I + J) <= nmax) * ~((I==0) * (J==0)) # [1:] avoids (0, 0)
I, J = [thing[keep].flatten() for thing in (I, J)]

i_super = I * k + J * (-l)
j_super = I * l + J * (k+l)

u_dots_super = np.vstack((i_super + 0.5*j_super, r3o2 * j_super))

plt.figure()
x, y = u_dots
plt.scatter(x, y, s=12)
plt.gca().set_aspect('equal')
xlims = plt.xlim()
ylims = plt.ylim()
x, y = u_dots_super
plt.scatter(x, y, s=32, c='red')
plt.gca().set_aspect('equal')
plt.show()