From c9834b8b3d6209a9f41f3d56e1ca85828327fd9d Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Fri, 17 Sep 2021 12:55:16 -0500
Subject: [PATCH] Add Matt's sphere_sample
---
pytools/__init__.py | 51 ++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 50 insertions(+), 1 deletion(-)
diff --git a/pytools/__init__.py b/pytools/__init__.py
index ad21b25..64a896b 100644
--- a/pytools/__init__.py
+++ b/pytools/__init__.py
@@ -2,7 +2,10 @@
# (Yes, it has a point!)
-__copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2009-2013 Andreas Kloeckner
+Copyright (C) 2020 Matt Wala
+"""
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -165,6 +168,11 @@ Hashing
.. autofunction:: unordered_hash
+Sampling
+--------
+
+.. autofunction:: sphere_sample
+
Type Variables Used
-------------------
@@ -2691,6 +2699,47 @@ def unordered_hash(hash_instance, iterable, hash_constructor=None):
# }}}
+# {{{ sphere_sample
+
+def sphere_sample(npoints_approx, r=1):
+ """Generate points regularly distributed on a sphere.
+ Based on: https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
+ Returns an array of shape (npoints, 3). npoints does not generally equally
+ npoints_approx.
+ """
+
+ import numpy as np
+ points = []
+
+ count = 0
+ a = 4 * np.pi / npoints_approx
+ d = a ** 0.5
+ M_theta = int(np.ceil(np.pi / d)) # noqa: N806
+ d_theta = np.pi / M_theta
+ d_phi = a / d_theta
+
+ for m in range(M_theta):
+ theta = np.pi * (m + 0.5) / M_theta
+ M_phi = int(np.ceil((2 * np.pi * np.sin(theta) / d_phi))) # noqa: N806
+ for n in range(M_phi):
+ phi = 2 * np.pi * n / M_phi
+ points.append([r * np.sin(theta) * np.cos(phi),
+ r * np.sin(theta) * np.sin(phi),
+ r * np.cos(theta)])
+ count += 1
+
+ # Add poles.
+ for i in range(3):
+ for sign in [-1, +1]:
+ pole = np.zeros(3)
+ pole[i] = r * sign
+ points.append(pole)
+
+ return np.array(points)
+
+# }}}
+
+
def _test():
import doctest
doctest.testmod()
--
GitLab