如何用圆方程和周期边界条件画圆[Python]？

Question

我要写一个二维伊辛模型模拟，在那里我不忽略远处邻居的影响，所以我想计算一个圆中的自旋。在

我写了一个简单的函数，它可以得到一个网格的元素，这些元素在一个圆中。在

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

它就像一个符咒，但是如果它是网格的一部分，它会把圆的一些部分切割掉。对于周期性边界条件，我该如何求解？在

提前谢谢！在

Answer 1

这里有一个解决方案。它使用numpy数组并以这样一种方式移动网格：可以很容易地将给定点周围半径r内的所有元素相加。一些有用的问题是

• How to apply a disc shaped mask to an array
• How to select a window from a numpy array with periodic boundary-conditions

以下代码中有一个限制，即2*半径+1必须小于或等于网格的最小形状。在

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