function grid(n,p) {
	return(colshape(J(1,p,1::n),1),J(n,1,1::p))
}

n = 20
a = grid(n,n)
a = a,sqrt(a[.,1]:^2+a[.,2]:^2)
a[selectindex(floor(a[.,3]):==a[.,3] :& a[.,3]:<=n),]