From b33a8e7ea8cd50815c79290bc60f98419ba1dfd6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Robert?= <cr52@protonmail.com>
Date: Sun, 6 Feb 2022 22:37:44 +0100
Subject: [PATCH] BUG: fix a bug in pixelize_cylinder where some pixels were
 skipped around theta=PI/2

---
 yt/utilities/lib/pixelization_routines.pyx | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx
index e628e260bc..0ef395b2d0 100644
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -523,9 +523,10 @@ def pixelize_cylinder(np.float64_t[:,:] buff,
                       np.float64_t[:] field,
                       extents):
 
-    cdef np.float64_t x, y, dx, dy, r0, theta0
+    cdef np.float64_t x, y, dx, dy, mdx, r0, theta0
     cdef np.float64_t rmax, x0, y0, x1, y1
-    cdef np.float64_t r_i, theta_i, dr_i, dtheta_i, dthetamin
+    cdef np.float64_t r_i, theta_i, dr_i, dtheta_i
+    cdef np.float64_t r_inc, theta_inc
     cdef np.float64_t costheta, sintheta
     cdef int i, pi, pj
 
@@ -535,6 +536,7 @@ def pixelize_cylinder(np.float64_t[:,:] buff,
     x0, x1, y0, y1 = extents
     dx = (x1 - x0) / buff.shape[1]
     dy = (y1 - y0) / buff.shape[0]
+    mdx = fmin(dx, dy)
     cdef np.float64_t rbounds[2]
     cdef np.float64_t corners[8]
     # Find our min and max r
@@ -558,9 +560,9 @@ def pixelize_cylinder(np.float64_t[:,:] buff,
         rbounds[0] = 0.0
     if y0 < 0 and y1 > 0:
         rbounds[0] = 0.0
-    dthetamin = dx / rmax
-    for i in range(radius.shape[0]):
 
+    r_inc = mdx
+    for i in range(radius.shape[0]):
         r0 = radius[i]
         theta0 = theta[i]
         dr_i = dradius[i]
@@ -569,15 +571,15 @@ def pixelize_cylinder(np.float64_t[:,:] buff,
         if r0 + dr_i < rbounds[0] or r0 - dr_i > rbounds[1]:
             continue
         theta_i = theta0 - dtheta_i
-        # Buffer of 0.5 here
-        dthetamin = 0.5*dx/(r0 + dr_i)
+        theta_inc = mdx/(r0 + dr_i)
+
         while theta_i < theta0 + dtheta_i:
             r_i = r0 - dr_i
             costheta = math.cos(theta_i)
             sintheta = math.sin(theta_i)
             while r_i < r0 + dr_i:
                 if rmax <= r_i:
-                    r_i += 0.5*dx
+                    r_i += r_inc
                     continue
                 y = r_i * costheta
                 x = r_i * sintheta
@@ -586,8 +588,8 @@ def pixelize_cylinder(np.float64_t[:,:] buff,
                 if pi >= 0 and pi < buff.shape[0] and \
                    pj >= 0 and pj < buff.shape[1]:
                     buff[pi, pj] = field[i]
-                r_i += 0.5*dx
-            theta_i += dthetamin
+                r_i += r_inc
+            theta_i += theta_inc
 
 cdef int aitoff_Lambda_btheta_to_xy(np.float64_t Lambda, np.float64_t btheta,
                                np.float64_t *x, np.float64_t *y) except -1: