Skip to content

Commit

Permalink
BUG: fix a bug in pixelize_cylinder where some pixels were skipped ar…
Browse files Browse the repository at this point in the history
…ound theta=PI/2
  • Loading branch information
neutrinoceros committed Feb 7, 2022
1 parent fe007fb commit b33a8e7
Showing 1 changed file with 11 additions and 9 deletions.
20 changes: 11 additions & 9 deletions yt/utilities/lib/pixelization_routines.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -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:
Expand Down

0 comments on commit b33a8e7

Please sign in to comment.