#!/usr/bin/env python
    # -*- coding: utf-8 -*-

    import numpy as np
    import matplotlib.pyplot as plt
    import itertools

    def get_points(a, nmax, theta_min=None, theta_max=None, remove_center=True):
        uvecs = np.array([[1, 0], [0.5, np.sqrt(3)/2.]])
        i = np.arange(-nmax, nmax+1)
        I, J = [thing.flatten() for thing in np.meshgrid(i, i)]
        keep = np.abs(I + J) <= nmax
        if remove_center:
            keep *= (~((I == 0) * (J == 0)))
        pairs = np.stack([thing[keep] for thing in (I, J)], axis=1)
        upoints = (pairs[..., None] * uvecs).sum(axis=1)
        thetas = np.arctan2(upoints[:, 1], upoints[:, 0])
        keep = np.ones_like(thetas, dtype=bool)
        if theta_min != None:
            keep *= thetas >= (np.pi/180) * theta_min
        if theta_max != None:
            keep *= thetas <= (np.pi/180) * theta_max
        return a * upoints[keep], pairs[keep]

    r3o2 = np.sqrt(3)/2.
    a, nmax = 1.0, 10

    all_points, all_pairs = get_points(a=a, nmax=nmax)
    points, pairs = get_points(a=a, nmax=nmax, theta_min=-0.01, theta_max=30.01)

    if True:
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        x, y = all_points.T
        ax.plot(x, y, '.k')
        x, y = points.T
        ax.plot(x, y, 'or')
        ax.plot([0, nmax*a], [0, 0], '-b')
        ax.plot([0, nmax*a*r3o2], [0, nmax*a*0.5], '-b')
        ax.set_aspect('equal')
        ax.axes.xaxis.set_visible(False)
        ax.axes.yaxis.set_visible(False)
        plt.show()