def dispersion():
    MAP_CENTER = [8192/2, 8192/2]
    param = [1,1]
    frame = 50
    axis='y'
    sim=sims[0]
    # Choose a random point (py,px) in a 8192x8192 map
    py=np.random.randint(8192)
    px=np.random.randint(8192)
    at = access_thing(sim,res,frame,param[0],param[1], machine=machine)

    # load stokes maps and calculate \phi in radians
    Q = at.get_fits_array('Q%s_n0-%04d_p-%d'%(axis,param[0],param[1]), axis)
    U = at.get_fits_array('U%s_n0-%04d_p-%d'%(axis,param[0],param[1]), axis)
    _phi = 0.5*np.arctan2(U,Q)

    SIZE = int(sys.argv[1])                 # size of cutout around (py,px)
    c = (SIZE/2,SIZE/2)                     # center of cutout
    max_l = SIZE/2
    l_map=np.zeros([SIZE,SIZE]).astype(int)
    l_list = range(1,max_l)
    l_masks = {}                            # dict {l:l_mask}
    nl = {}                                 # dict {l:number of pixels at a distance l from the center}
    cos_delta_phi_l = {}                    # dict {l: Sum of all cos(\Delta\phi) at a given l}
    dispersion_function = []                # 1 - <cos(\Delta\phi>

    # Make a map of the distance l from the center (py,px) to each pint (i,j) in the map
    # along with masks for each l that is True for all points (i,j) at a distance l from the center
    # and False for all points that are not a distance l away.
    # here is what l_map looks like when size = 15
    #[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    # [0, 0, 0, 0, 0, 6, 6, 6, 6, 6, 0, 0, 0, 0, 0],
    # [0, 0, 0, 6, 6, 5, 5, 5, 5, 5, 6, 6, 0, 0, 0],
    # [0, 0, 6, 6, 5, 4, 4, 4, 4, 4, 5, 6, 6, 0, 0],
    # [0, 0, 6, 5, 4, 4, 3, 3, 3, 4, 4, 5, 6, 0, 0],
    # [0, 6, 5, 4, 4, 3, 2, 2, 2, 3, 4, 4, 5, 6, 0],
    # [0, 6, 5, 4, 3, 2, 1, 1, 1, 2, 3, 4, 5, 6, 0],
    # [0, 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6, 0],
    # [0, 6, 5, 4, 3, 2, 1, 1, 1, 2, 3, 4, 5, 6, 0],
    # [0, 6, 5, 4, 4, 3, 2, 2, 2, 3, 4, 4, 5, 6, 0],
    # [0, 0, 6, 5, 4, 4, 3, 3, 3, 4, 4, 5, 6, 0, 0],
    # [0, 0, 6, 6, 5, 4, 4, 4, 4, 4, 5, 6, 6, 0, 0],
    # [0, 0, 0, 6, 6, 5, 5, 5, 5, 5, 6, 6, 0, 0, 0],
    # [0, 0, 0, 0, 0, 6, 6, 6, 6, 6, 0, 0, 0, 0, 0],
    # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
    y,x = np.mgrid[0:SIZE,0:SIZE]
    for l in l_list:
        l_mask = np.logical_and(((x-c[1])**2 + (y-c[0])**2 >= (l - 0.5)**2),((x-c[1])**2 + (y-c[0])**2 < (l + 0.5)**2))
        l_masks[l] = l_mask
        l_map[l_mask] = l
        nl[l] = len(l_map[l_map==l])
        cos_delta_phi_l[l] = 0


    print "calculating dispersion..."
    for i in range(SIZE):
        print "i =", i
        for j in range(SIZE):
            '''
            ## ROLL ARRAY
            # move point to center of box
            phi = np.roll(_phi, MAP_CENTER[0] - (py-i), axis=0)
            phi = np.roll( phi, MAP_CENTER[1] - (px-j), axis=1)

            # cutout a box in the stokes maps centered at (py-i,px-j)
            phi = phi[MAP_CENTER[0]-SIZE/2 : MAP_CENTER[0]+SIZE/2, MAP_CENTER[1]-SIZE/2 : MAP_CENTER[1]+SIZE/2]
            '''
            ## CUTOUT WITHOUT ROLLING ARRAY
            phi = _phi[(py-i)-SIZE/2 : (py-i)+SIZE/2, (px-j)-SIZE/2 : (px-j)+SIZE/2]
            phi_c = phi[c]     # phi_c = phi at the center of the cutout

            # calculate cos(\Delta\phi) for each pixel at each distance l
            # sum all the cos(\Delta\phi) in a given l bin
            for l in l_list:
                phi_l = phi[phi * l_masks[l] != 0]
                cos_delta_phi_l[l] += np.sum(np.cos(phi_c*np.ones(phi_l.shape) - phi_l))

    # calclate 1 - <cos(\Delta\phi>
    # there are nl*SIZE*SIZE pixels in each l bin
    for l in l_list:
        dispersion_function.append(1 - 1./(nl[l]*SIZE*SIZE)*cos_delta_phi_l[l])

    plt.scatter(l_list, dispersion_function)
    plt.xlabel("l")
    plt.ylabel(r'1 - $<cos(\Delta\phi)>$')
    plt.title("Dispersion Function")
    plt.savefig("dispersion_%s_lmax-%d_point-%d-%d"%(axis, max_l, py, px))