import numpy as np

x = np.array([0, 1, 2])
y = np.array([10, 20])

xx, yy = np.meshgrid(x, y, indexing="xy")
surface = xx + yy

ii, jj = np.meshgrid(x, y, indexing="ij")

xs, ys = np.meshgrid(x, y, indexing="xy", sparse=True)
sparse_surface = xs + ys

assert xx.shape == (2, 3)
assert yy.shape == (2, 3)
assert surface.shape == (2, 3)
assert ii.shape == (3, 2)
assert jj.shape == (3, 2)
assert xs.shape == (1, 3)
assert ys.shape == (2, 1)
assert np.array_equal(surface, sparse_surface)

print("xy grid shapes:", xx.shape, yy.shape)
print("xx:")
print(xx)
print("yy:")
print(yy)
print("surface:")
print(surface)
print("ij grid shapes:", ii.shape, jj.shape)
print("sparse grid shapes:", xs.shape, ys.shape)
print("sparse matches dense:", np.array_equal(surface, sparse_surface))