How to draw a circle with circle equation and periodic boundary condition [Python]?

Question

I have to write a 2D Ising model simulation, where I don't neglect the effect of the distant neighbors, so I want to count the spins in a circle.

I've written a simple function which can get the elements of a grid, which are in a circle.

def countInCircle(g, x, y, r):
    spinSum = 0
    for R in range(0, r + 1, 1):
        for i in range(0, g.shape[0], 1):
            for j in range(0, g.shape[1], 1):
                if ((i - x) ** 2 + (j - y) ** 2) == R:
                    spinSum = spinSum + g[i][j]

    return spinSum

It works like a charm, but it cuts down some parts of the circle if it's ouf of the grid. How should I solve this for periodic boundary condition?

Thanks in advance!

Unfortunately it's unclear what you are asking. What does "it cuts down some parts of the circle if it's ouf of the grid" mean? Where does the boundary condition come into play? Are `x`,`y` integers? What about the pixel `g[1][1]` which is sqrt(2) distant from `(0,0)`, will that ever be accounted for? — ImportanceOfBeingErnest, Dec 05 '16 at 19:21
Well, by "t cuts down some parts of the circle if it's ouf of the grid" I mean that if I want to count in a circle with r=10 radius, and 2,2 center, the function "cuts down" the top and the left of the circle(because that would be out of the grid), and that's a problem, because: I should use periodic boundary condition: the top row of the grid is "adjacent" to the bottom of the grid, and the left side of the grid is "adjacent" to the right side of the grid. So the top of the circle(which is cutted down currently) should be "come in" at the bottom of the grid, and the left side on the right — Martin, Dec 05 '16 at 19:29
It haven't fit in the 585 characters. Yes, x,y are integers. And yes, it will be accounted. Because that's on a circle with r=1, so if the circle has (0,0) origin, so the condition will be 1^2+1^2==1. But for R=2, it's 1^2+1^2==2. But now as I've written this down, I think I have a bug in my code. The condition should have R^2. Actually I don't know why, but it works if I use R, but it's not if I use R^2. But my main problem is the periodic boundary condition — Martin, Dec 05 '16 at 19:57
Indeed. What I mean is the following: If your radius is `r=2` you should sum over 13 sites, but your code sums over 9 sites only (the cross centered at (x,y)) and neglects the four intermediate pixels. Independent of if it's R or R^2 you only count those values which give an integer - which is counterintuitive to me, as the spins should interact isotropically. — ImportanceOfBeingErnest, Dec 05 '16 at 20:09
Now I see the problem! Thank you! Well, this can be solved, if the outer for runs only on the powers of two. So for the radius r=2 the if gets 4, and it is true for the other 4 pixels too. But unfortunately my main problem is still the periodic boundary condition. — Martin, Dec 05 '16 at 20:35

score 0 · Accepted Answer · edited May 23 '17 at 12:08

Here is a solution. It uses numpy arrays and shifts the grid in such a manner that one can easily sum up all elements within a radius r around the given point (x,y). Some useful questions are

There is one restriction in the following code which is that 2*radius+1 must be smaller or equal the minimum shape of the grid.

import numpy as np

def countInCircle(grid, x, y, r):
    #restriction: 2*r+1 < min(g.shape)
    shifted_grid=np.roll(np.roll(grid,shift=-x+r,axis=0),shift=-y+r,axis=1)
    Y,X = np.ogrid[-r: r+1, -r: r+1]    
    mask = X**2+Y**2 <= r**2
    return np.sum(shifted_grid[mask])

g  = np.ones((5,5))    
s  = countInCircle(g, 0, 0, 2)
print "s = ", s
# setting r=2 and summing all ones around (0,0) gives 13. Works fine. 

# setting some spin (-1,0,1) particles
g2 = np.random.randint(-1,2, size=(10,15))
print g2
s  = countInCircle(g2, 3,9, r=3)
print "s = ", s

It works, and it's much more faster than my code too. (in the meantime, I've modified my code too, and it seems to work, except the periodic boundary condition. I've removed the outer for, and changed the inner condition to math.sqrt((i - x) ** 2 + (j - y) ** 2) <= r. It gives 13 if I use r=2, etc.) Thank you for your answer! You've helped me a lot! — Martin, Dec 06 '16 at 16:52

How to draw a circle with circle equation and periodic boundary condition [Python]?

1 Answers1