function D2Ix = ortho( D2Ix, TS, lim=0.001 )

# Usage: D2Ix = ortho( D2Ix, TS)
# improve orthogonality for interpolation operator

ZZ2 = D2Ix' * D2Ix;
ZZ2 .*= (abs(ZZ2)>0.01);
D2Ix -= 0.5*(D2Ix*diag(TS)*ZZ2-D2Ix);
if(lim) ; D2Ix.*=(abs(D2Ix)>lim); endif
D2Ix = diag(1./sum(D2Ix,2))*D2Ix;
endfunction
