function [Fp,H9] = hexnet(NB=20,ND=60)

# Usage: [Fp,H9] = hexnet(NB=20,ND=60)
# construct a foveated grid of mostly hexagonal points
# at locations Fp and adjacency matrix H
# NB is the side length of the inner part
# ND is the number of outer circles (should be even)

NC=6*NB;
#global NB; global NC; global ND;
#NB = 20; NC=6*NB; ND = 3*NB; #ND=5; # ND==60 ist ca. Fovea-Größe
Nx = ceil(NB*sqrt(4/3)); Nx2 = 2*Nx+1;
Nx = NB; Nx2 = Nx+1;
Nh = 3*Nx*(Nx+1)+1;
[Dj,Dk] = meshgrid( 0:Nx, 0:Nx);
sel = Dj>Dk;
# lattice generators
v1=1; v2=-0.5+0.5i*sqrt(3);
Dp = Dj*v1+Dk*v2;

H7d = reshape(1:(Nx2*Nx2),[Nx2,Nx2]);
H71 = sparse( 1:(Nx2*Nx2), vec(shift(H7d,1,1)),ones(1,Nx2*Nx2))(sel,sel);
H72 = sparse( 1:(Nx2*Nx2), vec(shift(H7d,1,2)),ones(1,Nx2*Nx2))(sel,sel);
H73 = sparse( 1:(Nx2*Nx2), vec(shift(shift(H7d,1,1),1,2)),ones(1,Nx2*Nx2))(sel,sel);
H7 = H71 + H72 + H73;

Dp = Dp(sel);
#gplot((H7),[real(Dp),imag(Dp)])

[Ej,Ek] = meshgrid(0:(NB-1),0:(ND-1));
#Ep=vec(exp(-1i*(Ej'*v1+Ek'*v2+0.5-NB)/NC*2*pi));
Ep=vec(exp(-1i*((Ej'+ceil(Ek'/2))*v1+Ek'*v2+0.5-NB)/NC*2*pi));

H8d = reshape(1:NB*ND,[NB,ND]);
r1 = 1:NB*ND;
r2  =shift(r1,1); r2(1+NB*(0:ND-1)) = 0;
H81 = sparse( r2(r2>0),r1(r2>0),1,NB*ND,NB*ND);
H82 = sparse( NB+1:NB*ND, 1:(ND-1)*NB,1,NB*ND,NB*ND);
H84 = shift(H8d,1,1);
H84(:,1:2:ND) = shift(H84(:,1:2:ND),-1,1);
H83 = sparse(H84(2:end,1:end-1),H84(2:end,2:end),1,NB*ND,NB*ND);
r3 = H84(1,:);

#H81 = sparse( 1:NB*ND,vec(shift(H8d,1,1)), kron(ones(1,ND),[0,ones(1,NB-1)]));
#H82 = sparse( 1:NB*ND,vec(shift(H8d,1,2)), [zeros(1,NB),ones(1,NB*(ND-1))]);
#H83 = sparse( 1:NB*ND,vec(shift(shift(H8d,1,1),1,2)),kron([0,ones(1,ND-1)],[0,ones(1,NB-1)]),NB*ND,NB*ND);

H8 = H81 + H82 + H83;
#gplot(H8,[real(Ep),imag(Ep)])

Fp = [Dp/NB;vec(Ep)]*exp(0.01i);
#Fp *= sqrt(2/sqrt(3)); # is falsch
lF = length(Fp);
H91=blkdiag(H73,H81);
H92=blkdiag(H72',H82);
H93=blkdiag(H71,H83);
# links zw. innen und außen herstellen

H92(sub2ind(size(H92),length(Dp)+1-(1:NB),length(Dp)+(1:NB))) = 1;
H93(sub2ind(size(H92),length(Dp)+1-(1:NB-1),length(Dp)+1+(1:NB-1))) = 1;

H9=H91 + H92 + H93;
#H9=kron(eye(6),H9); # leider ist das Ergebnis nicht sparse
H9=blkdiag(H9,H9,H9,H9,H9,H9);
Fp = kron(exp(-1i*pi/3*(0:5))',Fp);
r1=((2:NB+1) .*(1:NB))/2;
r2=((1:NB) .*(0:NB-1))/2+1;
q1=(0:ND)*NB+length(Dp);
r3 += length(Dp);
for j=0:4;
H9(sub2ind(size(H9),r2(NB)+lF,1+q1(1))) = 1;
H9(sub2ind(size(H9),r1,r2+lF)) = 1;
H9(sub2ind(size(H9),r1(2:NB),r2(1:NB-1)+lF)) = 1;
H9(sub2ind(size(H9),1+q1(1:ND),NB+q1(1:ND)+lF)) = 1;
H9(sub2ind(size(H9),r3(1:2:end),r3(2:2:end)+lF)) =1;
H9(sub2ind(size(H9),r3(3:2:end),r3(2:2:end-2)+lF)) =1;
#H9(sub2ind(size(H9),1+q1(2:ND),q1(2:ND)+lF)) = 1;
r1 += lF;
r2 += lF;
r3 += lF;
q1 += lF;
endfor;
H9(sub2ind(size(H9),r2(NB)-5*lF,1+q1(1))) = 1;
H9(sub2ind(size(H9),r1,r2-5*lF)) = 1;
H9(sub2ind(size(H9),r1(2:NB),r2(1:NB-1)-5*lF)) = 1;
H9(sub2ind(size(H9),1+q1(1:ND),NB+q1(1:ND)-5*lF)) = 1;
H9(sub2ind(size(H9),r3(1:2:end),r3(2:2:end)-5*lF)) =1;
H9(sub2ind(size(H9),r3(3:2:end),r3(2:2:end-2)-5*lF)) =1;
#H9(sub2ind(size(H9),1+q1(2:ND),q1(2:ND)-5*lF)) = 1;

#Fpw = Fp(end-(0:(NC-1)));
#for j=1:50; 
#Fp = diag(1./sum(H9)) * H9 * Fp; Fp(end-(0:(NC-1))) = Fpw;
#endfor;
lz = length(Fp)+1;
H9 = H9 +H9';
H9(lz,1+lF*(0:5))=1;
H9(1+lF*(0:5),lz)=1;
Fp(lz) = 0;

sFp = find(abs(Fp) < abs(Fp(lz-1)-0.1));
for j=1:50;
Fp(sFp) = (diag(1./sum(H9)) * H9 * Fp)(sFp);
endfor;

endfunction
#gplot(H9,[real(Fp),imag(Fp)])
