Sinusoid coupled functions with natural frequency "w" with eight neighbors coupling Grid (except for the border that behaves like a bi-dimensional Grid)

kuramotoGrid.py

import numpy as np
import scipy
from scipy.integrate import *
import matplotlib.pyplot as plot
sin = np.sin

def kuramotoGrid(x,t0, K, w):
    l = [len(w[1,:]),len(w[:,1])];
    x = x.reshape((l[1],l[0]));
# l[0] is horizontal, columns, l[1] is vertical, rows, dimension
    N = int(l[0]*l[1]);
    y = np.zeros((l[1],l[0]));
####################
# Contourn Equations
####################
#   The 4 pixels that are in the corners
    y[0,0] = w[0,0] + K*( sin(x[0,1]-x[0,0]) + sin(x[1,0]-x[0,0]) + sin(x[1,1]-x[0,0]) )/3
    y[0,l[0]-1] = w[0,l[0]-1] + K*( sin(x[0,l[0]-2]-x[0,l[0]-1]) + sin(x[1,l[0]-1]-x[0,l[0]-1]) + sin(x[1,l[0]-2]-x[0,l[0]-1]) )/3
    y[l[1]-1,0] = w[l[1]-1,0] + K*(sin( x[l[1]-2,0] - x[l[1]-1,0]) + sin(x[l[1]-2,1] - x[l[1]-1,0]) + sin(x[l[1]-1,1]- x[l[1]-1,0]) )/3
    y[l[1]-1,l[0]-1] = w[l[1]-1,l[0]-1] + K*(sin(x[l[1]-1,l[0]-2] - x[l[1]-1,l[0]-1]) + sin(x[l[1]-2,l[0]-1] - x[l[1]-1,l[0]-1]) + sin(x[l[1]-2,l[0]-2] - x[l[1]-1,l[0]-1]))/3
#   The border pixes that are not in the corner
    y[0,1:l[0]-1] = w[0,1:l[0]-1] + K*( sin(x[0,0:l[0]-2]-x[0,1:l[0]-1]) + sin(x[0,2:l[0]]-x[0,1:l[0]-1]) + sin(x[1,0:l[0]-2]-x[0,1:l[0]-1]) + sin(x[1,1:l[0]-1]-x[0,1:l[0]-1] ) + sin(x[1,2:l[0]]-x[0,1:l[0]-1]) )/5 ;
    y[1:l[1]-1,0] = w[1:l[1]-1,0] + K*( sin(x[0:l[1]-2,0] -x[1:l[1]-1,0] ) + sin(x[2:l[1],0] -x[1:l[1]-1,0]  ) + sin(x[0:l[1]-2,1] -x[1:l[1]-1,0]  ) + sin(x[1:l[1]-1,1] -x[1:l[1]-1,0]  ) + sin(x[2:l[1],1] -x[1:l[1]-1,0]  ) )/5
    y[l[1]-1,1:l[0]-1] = w[l[1]-1,1:l[0]-1] + K*( sin(x[l[1]-1,0:l[0]-2] - x[l[1]-1,1:l[0]-1] ) + sin(x[l[1]-1,2:l[0]] - x[l[1]-1,1:l[0]-1] ) + sin(x[l[1]-2,0:l[0]-2]- x[l[1]-1,1:l[0]-1] ) + sin(x[l[1]-2,1:l[0]-1] - x[l[1]-1,1:l[0]-1] ) + sin(x[l[1]-2,2:l[0]] - x[l[1]-1,1:l[0]-1] ) )/5
    y[1:l[1]-1,l[0]-1] = w[1:l[1]-1,l[0]-1] + K*( sin(x[0:l[1]-2,l[0]-1] -x[1:l[1]-1,l[0]-1] ) + sin(x[2:l[1],l[0]-1] -x[1:l[1]-1,l[0]-1]) + sin(x[0:l[1]-2,l[0]-2]-x[1:l[1]-1,l[0]-1] ) + sin(x[1:l[1]-1,l[0]-2]-x[1:l[1]-1,l[0]-1] ) + sin(x[2:l[1],l[0]-2]-x[1:l[1]-1,l[0]-1] ) )/5
##############################################
# Central pixels equations
    y[1:l[1]-1,1:l[0]-1] = w[1:l[1]-1,1:l[0]-1] + K*(sin(x[0:l[1]-2,1:l[0]-1]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[0:l[1]-2,0:l[0]-2]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[0:l[1]-2,2:l[0]]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[1:l[1]-1,0:l[0]-2]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[1:l[1]-1,2:l[0]]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[2:l[1],0:l[0]-2]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[2:l[1],1:l[0]-1]-x[1:l[1]-1,1:l[0]-1] ) + sin(x[2:l[1],2:l[0]]-x[1:l[1]-1,1:l[0]-1]))/8

    return y.reshape(-1)

TODO: add heterogenities in the coupling