function index = n_hcomp2index (N) if N < 1 then N = 1 end m = ceil(N/2) index=[0] xs=3 for j = 0:(m-1) for i = 0:(m-1) x = xs + 3*i index = [index round(max(x^2+3/4*j^2,(x-0.5)^2 +3/4*(j+1)^2 ))] end if modulo(j,3) then xs=xs+1.5 else xs=xs-1.5 end end index = unique(gsort(index,"g","i")) index=index(1:N) endfunction