From 62dd66be7938f23bda322a88d863c07d2251eab2 Mon Sep 17 00:00:00 2001
From: "Paul M. Sutter" <pmsutter@gmail.com>
Date: Sat, 8 Jun 2024 16:18:33 +0100
Subject: [PATCH] added vornoi plotting routine

---
 python_source/voidUtil/plotUtil.py | 89 +++++++++++++++++++++++++++++-
 1 file changed, 88 insertions(+), 1 deletion(-)

diff --git a/python_source/voidUtil/plotUtil.py b/python_source/voidUtil/plotUtil.py
index 1132b59..c6cf85b 100644
--- a/python_source/voidUtil/plotUtil.py
+++ b/python_source/voidUtil/plotUtil.py
@@ -17,13 +17,16 @@
 #   with this program; if not, write to the Free Software Foundation, Inc.,
 #   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #+
-__all__=['plotNumberFunction','plotEllipDist','plotVoidCells']
+__all__=['plotNumberFunction','plotEllipDist','plotVoidCells','plotVoronoi']
 
 from backend.classes import *
 from .plotDefs import *
 import numpy as np
 import os
 import pylab as plt
+from scipy.spatial import ConvexHull
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+rng = np.random.default_rng(11)
 import backend.cosmologyTools as vp
 import backend.surveyTools as sp
 from voidUtil import getArray, shiftPart, getVoidPart
@@ -298,3 +301,87 @@ def plotVoidCells(catalog,
   plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
 
   return
+
+
+# -----------------------------------------------------------------------------
+def plotVoronoi(catalog,
+              voidID,
+              figDir="./",
+              plotName="voronoi",
+              plotCenter = (0,0,0),
+              plotWidth = (100,100,25),
+              sliceWidth=250):
+
+# plots the voronoi cells belonging to a void (or a slice if no void given)
+#   catalog: void catalog
+#   voidID: ID of void to plot (-1 to use all particles)
+#   figDir: output directory for figures
+#   plotName: name to prefix to all outputs
+
+  print("Plotting voronoi cells...")
+  sample = catalog.sampleInfo
+
+  plt.clf()
+  fig = plt.figure()
+  ax = fig.add_subplot(111, projection='3d')
+
+
+  print("  Selecting particles...")
+  if (voidID == -1):
+    xmin = plotCenter[0]-plotWidth[0]/2.
+    xmax = plotCenter[0]+plotWidth[0]/2.
+    ymin = plotCenter[1]-plotWidth[1]/2.
+    ymax = plotCenter[1]+plotWidth[1]/2.
+    zmin = plotCenter[2]-plotWidth[2]/2.
+    zmax = plotCenter[2]+plotWidth[2]/2.
+
+    partList = []
+    for p in catalog.part:
+      #print(p.x, p.y, p.z)
+      #print(xmin,xmax,ymin,ymax,zmin,zmax)
+      if (p.x >= xmin and p.x <= xmax and \
+          p.y >= ymin and p.y <= ymax and \
+          p.z >= zmin and p.z <= zmax):
+        partList.append(p)
+
+    #ax.set_xlim(xmin,xmax)
+    #ax.set_ylim(ymin,ymax)
+    #ax_set.zlim(zmin,zmax)
+  else:
+    partList = getVoidPart(catalog, voidID)
+    plotCenter = catalog.voids[voidID].macrocenter
+    plotRadius = catalog.voids[voidID].radius
+
+    xmin = plotCenter[0]-plotRadius
+    xmax = plotCenter[0]+plotRadius
+    ymin = plotCenter[1]-plotRadius
+    ymax = plotCenter[1]+plotRadius
+    zmin = plotCenter[2]-plotRadius
+    zmax = plotCenter[2]+plotRadius
+
+    #ax.set_xlim(xmin,xmax)
+    #ax.set_ylim(ymin,ymax)
+    #ax_set_zlim(zmin,zmax)
+
+  print("  Working with %d particles..." % len(partList))
+  print("  Computing convex hulls and creating polygons...")
+  for part in partList:
+    adjsPart = []
+    for adj in part.adjs:
+      adjsPart.append(catalog.partPos[adj])
+
+    hull = ConvexHull(adjsPart)
+    polygon = Poly3DCollection(hull.points[hull.simplices], alpha=0.25,
+                               facecolors=rng.uniform(0,1,3),
+                               linewidths=0.01, edgecolors='gray')
+    ax.add_collection3d(polygon)
+    
+    ax.scatter(part.x, part.y, part.z, s=0.01, color='black')
+
+  print("  Saving...")
+  plt.savefig(figDir+"/fig_"+plotName+".pdf", bbox_inches="tight")
+  plt.savefig(figDir+"/fig_"+plotName+".eps", bbox_inches="tight")
+  plt.savefig(figDir+"/fig_"+plotName+".png", bbox_inches="tight")
+
+  print("  Done!")
+  return