let periodic distance =
    let length = 2*range +1 in
    let init x y =
      Femtolib.Array_simulate_2d.init length length (fun i j -> i+x-range, j+y-range)
      |> Array.to_list
      |> List.filter (fun (i,j) -> (i <> x || j <> y)  &&  distance (i,j) (x,y) <= range)
      |> Array.of_list
      |> Array.map (fun (i,j) -> Femtolib.Array_simulate_2d.to_1d width height
                      ((i+width) mod width) ((j+height) mod height)) in
    Femtolib.Array_simulate_2d.init width height init