function H9x = intJ(H9,Fp,lim=0.70)

# Usage: H9x = intJ(H9,Fp,lim)
# this function computes the part of an adjacency matrix for
# a refinement of the grid between the old and new graph nodes
# If present, parameter determines the size threshold w.r.t. first
# distance in Fp

[i,j,v]= find(H9);
# Inzidenzmatrix aus Adjazenz berechnen
trii = sparse(1:sum(i<j),i(i<j),1,sum(sum(i<j)),size(H9)(1));
trii += sparse(1:sum(i<j),j(i<j),1,sum(sum(i<j)),size(H9)(1));
# H9 == trii'*trii -diag(diag(trii'*trii))
tri0 = trii*trii' -diag(diag(trii*trii'));
tmp = ((tri0*trii)==2);
tmp .*= spones(tmp)-trii;
[i,j] = find(tmp);
# j (i==220) ergibt die beiden Dreiecke zu find(trii(220,:))
[k,l]=find(trii);
l= reshape(sortrows([k,l])(:,2),2,[]);
z1(i) = j;
ij = [i,j](z1(i)' != j,:);
z2(ij(:,1)) = ij(:,2);
#ij = ij(z2(ij(:,1))' != ij(:,2),:);
#z3(ij(:,1)) = ij(:,2);

tmp = [z1, z2(z2>0); l,l(:,z2>0)];
tmp = unique(sort(tmp)',"rows");

tribar = mean(Fp(tmp),2);
tstd = std(Fp(tmp),1,2);
Fsel = find(tstd > lim*abs(Fp(2)-Fp(1)));
tribar = tribar(Fsel); 
tmp=tmp(Fsel,:);
#if(nargout>2) tristd =tstd(Fsel); endif
H9x = sparse(1:length(tribar),tmp(:,1),1,length(tribar),length(Fp));
H9x += sparse(1:length(tribar),tmp(:,2),1,length(tribar),length(Fp));
H9x += sparse(1:length(tribar),tmp(:,3),1,length(tribar),length(Fp));
endfunction
