#基爾霍夫衍射,衍射屏坐標范圍-dGrid:dGrid,光源坐標(0,0)
#簡單的矩孔衍射,dSource為光源到小孔的距離；dScreen為衍射屏到小孔距離
#dHole為矩孔網(wǎng)格尺寸；dGrid為衍射屏網(wǎng)格尺寸；nGrid為網(wǎng)格數(shù)目
def squareDiff(dSource=1,dScreen=1.4,dWave=1e-6,
            dHole=3e-5,nGrid=100,dGrid=1e-5):
    nX,nY = nGrid*np.array([1,1])
    axisX = np.arange(-nX,nX+1)*dGrid       
    axisY = np.arange(-nY,nY+1)*dGrid
    xAxis,yAxis = np.meshgrid(axisX,axisY)  #此為衍射屏的x坐標
    axisX = np.arange(-nX,nX+1)*dHole       
    axisY = np.arange(-nY,nY+1)*dHole
    xHole,yHole = np.meshgrid(axisX,axisY)  #此為矩孔的x坐標
    dArrS = np.sqrt(xHole**2+yHole**2+dSource**2)   #孔平面到光源的距離
    nSide = int(nGrid*2+1)              #格點個數(shù)
    pane = np.zeros([nSide,nSide])      #衍射屏強度
    for m in range(nSide):
        for n in range(nSide):
            dArr = np.sqrt((xHole-xAxis[m,n])**2+(yHole-yAxis[m,n])**2+dScreen**2)
            pane[m,n] = np.sum(
                np.cos(np.pi*2*(dArr+dArrS)/dWave)/dArr/dArrS)
    pane = np.abs(pane)
    pane = pane/np.max(np.array(pane))
    fig = plt.figure()
    ax = axd(fig)
    ax.plot_surface(xAxis,yAxis,pane)
    plt.show()
    return pane, xAxis, yAxis

#輸入對于M*M矩陣的第一個值到N*N矩陣的距離，返回M(m,n)的距離矩陣
def getDisMat(dMat,N,m,n):
    dMat = np.mat(dMat)
    A = dMat[1:m,1:n]
    B = dMat[1:m,0:N-n+1]
    C = dMat[0:N-m+1,1:n]
    D = dMat[0:N-m+1,0:N-n+1]
    return np.vstack((np.hstack((np.flip(A),np.flipud(B))),np.hstack((np.fliplr(C),D))))    #stack矩陣拼接，flip翻轉(zhuǎn)矩陣

    dArrS = np.sqrt(xHole**2+yHole**2+dSource**2) #孔平面到光源的距離
    dScreen = np.sqrt(                             #衍射平上第(0,0)個點的距離矩陣
        (xHole-xAxis[0,0])**2+(yHole-yAxis[0,0])**2+dScreen**2)
    nSide = int(nGrid*2+1)              #格點個數(shù)
    pane = np.zeros([nSide,nSide])      #衍射屏強度
    for m in range(nSide):
        for n in range(nSide):
            dArr = getDisMat(dScreen,nSide,m,n)