function h = etaijXdtree23(P,m,n,i,j,y,w1,w2)
global A C T
global LEVELS
h = 0;
LEVELS = repmat({[]},n,1);
Ti = getSubt(T,i);
jj = Ti(Ti>=j);
if isempty(jj)
  return
end
getLevels(T,1,0);
j0 = min(jj);
levi0 = getLevu(i);
levj0 = getLevu(j0);

Ej = getEj(T,j0);
IC = intersect(Ej(:,2),Ti);
DC = repmat([m],[1 length(C)]);
IZ = intersect(Ti,i+1:j0-1);
DZ = repmat([m],[1 length(IZ)]);

C0 = intersect(IC,LEVELS{levj0});
C1 = setdiff(IC,C0);
Z0 = setdiff(LEVELS{levj0},C0);
Z1 = setdiff(IZ,Z0);

DC0 = repmat([m],[1 length(C0)]);
DC1 = repmat([m],[1 length(C1)]);
DZ0 = repmat([m],[1 length(Z0)]);
DZ1 = repmat([m],[1 length(Z1)]);

for ci0=1:prod(DC0)
  xC0 = ind2coord(ci0,DC0);
  for ci1=1:prod(DC1)
    xC1 = ind2coord(ci1,DC1);
    % begin zeta
    Sum = 0;
    for zi0=1:prod(DZ0)
      xZ0 = ind2coord(zi0,DZ0);
      for zi1=1:prod(DZ1)
        xZ1 = ind2coord(zi1,DZ1);
        
        if isempty(Z0)
          pC1Z0 = 1;
        else
          pC1Z0 = getcondP0(C1,Z0,xC1,xZ0);
        end
        
        if isempty(Z1)
          % no buffer zone
          pcond1 = getcondP0([Z0 C0],[i],[xZ0 xC0],[w1]);
          pcond2 = getcondP0([Z0 C0],[i],[xZ0 xC0],[w2]);
          Sum = Sum + pC1Z0*(pcond1 - pcond2);
        else
          pC0Z1 = getcondP0(C0,Z1,xC0,xZ1);
          pcond1 = getcondP0([Z1 Z0],[i],[xZ1 xZ0],[w1]);
          pcond2 = getcondP0([Z1 Z0],[i],[xZ1 xZ0],[w2]);
          Sum = Sum + pC1Z0*pC0Z1*(pcond1 - pcond2);
        end

      end
    end
    h = h + .5*abs(Sum);
  end
end
return



function getLevels(T,u,ell)
global LEVELS
ell = ell+1;
LEVELS{ell}(end+1) = u;
ukids = find(T(u,:));
for vi=1:length(ukids)
  v = ukids(vi);
  getLevels(T,v,ell);
end
return

function [ell0,LEVu] = getLevu(u)
global LEVELS
for ell=1:length(LEVELS)
  if ismember(u,LEVELS{ell})
    LEVu = LEVELS{ell};
    ell0 = ell;
    return
  end
end
return


function Ej = getEj(T,j)
Ej = [];
[nodesi,nodesj]=find(T);
for e=1:length(nodesi)
  u = nodesi(e);
  v = nodesj(e);
  if ((u<j) & (v>=j))
  %if ((u<j))
    Ej = [Ej;[u v]];
  end
end
return


function Pc = getcondP0(jj,ii,x,y)
global P m n
% P = P{ X[jj]=x | X[ii]=y }
Pjoin = getPX(P,m,n,[y x],[ii jj]);
Pmarg = getPX(P,m,n,[y  ],[ii   ]);
Pc = Pjoin/(Pmarg+realmin);
return

function sT = getSubt(T,v)
% a list of v's descendents in T
% includes v itself
global SUBT
SUBT = [];
Pi0(T,v);
sT = SUBT;
return

function Pi0(T,v)
global SUBT
SUBT(end+1) = v;
vkids = find(T(v,:));
for vi=1:length(vkids)
  v1 = vkids(vi);
  Pi0(T,v1);
end