% WORKS
function h = etaijX300(P,m,n,i,j,X,y1,y2)
global A B C mobs mhid
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
Di  = repmat([m],[1 i]);
D1n = repmat([m],[1 n]);
hD  = repmat([mhid],[1 n]);
hDjn = repmat([mhid],[1 n-j+1]);
hDij = repmat([mhid],[1 j-i-1]);
hDi  = repmat([mhid],[1 i]);

Py1 = getPX(P,m,n,[X y1],1:i);
Py2 = getPX(P,m,n,[X y2],1:i);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    for hxi = 1:prod(hDjn)
      hxjn = ind2coord(hxi,hDjn);
      for hzi=1:prod(hDij)
        hzij = ind2coord(hzi,hDij);
        for hwi=1:prod(hDi)
          hw1i = ind2coord(hwi,hDi);
          Ph = Phid0([hw1i hzij hxjn]);
          Po = Pobs(hw1i(1:i-1),X,1,i-1);
          Po = Po * Pobs(hzij,zij,i+1,j-1);
          Po = Po * Pobs(hxjn,xjn,j,n);
          p1 = Ph * Po * B(y1,hw1i(i),i);
          p2 = Ph * Po * B(y2,hw1i(i),i);          
          
          Sum = Sum + p1/Py1 - p2/Py2;
        end
      end
    end
  end
  h = h + .5*abs(Sum);
end
return


function h = etaijX3(P,m,n,i,j,X,y1,y2)
global A B C mobs mhid
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
Di  = repmat([m],[1 i]);
D1n = repmat([m],[1 n]);
hD  = repmat([mhid],[1 n]);
hDjn = repmat([mhid],[1 n-j+1]);
hDij = repmat([mhid],[1 j-i-1]);
hDi  = repmat([mhid],[1 i]);

Py1 = getPX(P,m,n,[X y1],1:i);
Py2 = getPX(P,m,n,[X y2],1:i);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    for hxi = 1:prod(hDjn)
      hxjn = ind2coord(hxi,hDjn);
      for hzi=1:prod(hDij)          
        hzij = ind2coord(hzi,hDij);
        for hwi=1:prod(hDi)          
          hw1i = ind2coord(hwi,hDi);
          
          p = 1;
          %Ph = Phid([hw1i hzij hxjn],1,n);
          if 1
            if isempty(hzij)
              %p = p * Phid([hw1i hxjn],1,n);
              p = p * A(hxjn(1),hw1i(i),i); % [w x] tfira
            else
              %p = p * Phid([hw1i hzij hxjn],1,n);
              p = p * A(hzij(1),hw1i(i),i); % [w z] tfira
              p = p * A(hxjn(1),hzij(end),j-1); % [z x] tfira
            end
            p = p * Phid(hw1i,1,i);
            p = p * Phid(hxjn,j,n);
            p = p * Phid(hzij,i+1,j-1);
            Ph = p;
          end
          p = p * Pobs(hw1i(1:i-1),X,1,i-1);
          p = p * Pobs(hzij,zij,i+1,j-1);
          p = p * Pobs(hxjn,xjn,j,n);
          p1 = p * B(y1,hw1i(i),i);
          p2 = p * B(y2,hw1i(i),i);
          
          %Po = Pobs(hw1i(1:i-1),X,1,i-1)*Pobs(hzij,zij,i+1,j-1)*Pobs(hxjn,xjn,j,n);
          %p1 = Ph * Po * B(y1,hw1i(i),i);
          %p2 = Ph * Po * B(y2,hw1i(i),i);          
               
          

          
          Sum = Sum + p1/Py1 - p2/Py2;
        end
      end
    end
  end
  h = h + .5*abs(Sum);
end
return





%%%%%%%%%%%WORKS%%%%%%%%%%%%
function h = abcjunk(P,m,n,i,j,X,y1,y2)
global A B C mobs mhid
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
Di  = repmat([m],[1 i]);
D1n = repmat([m],[1 n]);
hD  = repmat([mhid],[1 n]);
hDjn = repmat([mhid],[1 n-j+1]);
hDij = repmat([mhid],[1 j-i-1]);
hDi  = repmat([mhid],[1 i]);

Py1 = getPX(P,m,n,[X y1],1:i);
Py2 = getPX(P,m,n,[X y2],1:i);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    for hxi = 1:prod(hDjn)
      hxjn = ind2coord(hxi,hDjn);
      for hzi=1:prod(hDij)
        hzij = ind2coord(hzi,hDij);
        for hwi=1:prod(hDi)
          hw1i = ind2coord(hwi,hDi);
          

          %Pro = pi(hxjn,xjn,j,n);      
          %if isempty(hzij)
          %  Pro = Pro * A(hxjn(1),hw1i(end),i); % [w x] tfira
          %else
          %  Pro = Pro * A(hzij(1),hw1i(end),i); % [w z] tfira
          %  Pro = Pro * A(hxjn(1),hzij(end),j-1); % [z x] tfira
          %  Pro = Pro * pi(hzij,zij,i+1,j-1);
          %end
          %p1 = Pro * pi(hw1i,[X y1],1,i);
          %p2 = Pro * pi(hw1i,[X y2],1,i);         
          
          p = 1;
          p = p * Phid(hw1i,1,i);
          if isempty(hzij)
            %p = p * Phid([hw1i hxjn],1,n);
            p = p * A(hxjn(1),hw1i(i),i); % [w x] tfira
          else
            %p = p * Phid([hw1i hzij hxjn],1,n);
            p = p * A(hzij(1),hw1i(i),i); % [w z] tfira
            p = p * A(hxjn(1),hzij(end),j-1); % [z x] tfira
          end
          p = p * Phid(hxjn,j,n);
          p = p * Pobs(hw1i(1:i-1),X,1,i-1);
          P = p * Pobs(hzij,zij,i+1,j-1);
          p = p * Pobs(hxjn,xjn,j,n);
          
          p1 = p * B(y1,hw1i(i),i);
          p2 = p * B(y2,hw1i(i),i);
          
          %Ph = Phid([hw1i hzij hxjn]);
          %Po = Pobs(hw1i(1:i-1),X,1,i-1);
          %Po = Po * Pobs(hzij,zij,i+1,j-1);
          %Po = Po * Pobs(hxjn,xjn,j,n);
          %p1 = Ph * Po * B(y1,hw1i(i),i);
          %p2 = Ph * Po * B(y2,hw1i(i),i);          
          
          Sum = Sum + p1/Py1 - p2/Py2;
        end
      end
    end
  end
  h = h + .5*abs(Sum);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




function h = etaijX30(P,m,n,i,j,X,y1,y2)
global A B C mobs mhid
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
Di  = repmat([m],[1 i]);
D1n = repmat([m],[1 n]);
hD  = repmat([mhid],[1 n]);
hDjn = repmat([mhid],[1 n-j+1]);
hDij = repmat([mhid],[1 j-i-1]);
hDi  = repmat([mhid],[1 i]);

Py1 = getPX(P,m,n,[X y1],1:i);
Py2 = getPX(P,m,n,[X y2],1:i);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    for hxi = 1:prod(hDjn)
      hxjn = ind2coord(hxi,hDjn);
      Pro = pi(hxjn,xjn,j,n);
      for hzi=1:prod(hDij)
        hzij = ind2coord(hzi,hDij);
        if ~isempty(hzij)
          Pro = Pro * A(hxjn(1),hzij(end),j-1); % [z x] tfira
          Pro = Pro * pi(hzij,zij,i+1,j-1);
        end
        for hwi=1:prod(hDi)
          hw1i = ind2coord(hwi,hDi);
          if isempty(hzij)
            Pro = Pro * A(hxjn(1),hw1i(end),i); % [w x] tfira
          else
            Pro = Pro * A(hzij(1),hw1i(end),i); % [w z] tfira
          end
          p1 = Pro * pi(hw1i,[X y1],1,i);
          p2 = Pro * pi(hw1i,[X y2],1,i);          
          
          Sum = Sum + p1/Py1 - p2/Py2;
        end
      end
    end
    h = h + .5*abs(Sum);
  end
end
return



function h = etaijX300(P,m,n,i,j,X,y1,y2)
global A B C mobs mhid
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
Di  = repmat([m],[1 i]);
D1n = repmat([m],[1 n]);
hD  = repmat([mhid],[1 n]);
hDjn = repmat([mhid],[1 n-j+1]);
hDij = repmat([mhid],[1 j-i-1]);
hDi  = repmat([mhid],[1 i]);

Py1 = getPX(P,m,n,[X y1],1:i);
Py2 = getPX(P,m,n,[X y2],1:i);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    for hxi = 1:prod(hDjn)
      hxjn = ind2coord(hxi,hDjn);
      Pro = pi(hxjn,xjn,j,n);      
      for hzi=1:prod(hDij)
        hzij = ind2coord(hzi,hDij);
        for hwi=1:prod(hDi)
          hw1i = ind2coord(hwi,hDi);
          

          if isempty(hzij)
            Pro = Pro * A(hxjn(1),hw1i(end),i); % [w x] tfira
          else
            Pro = Pro * A(hzij(1),hw1i(end),i); % [w z] tfira
            Pro = Pro * A(hxjn(1),hzij(end),j-1); % [z x] tfira
            Pro = Pro * pi(hzij,zij,i+1,j-1);
          end
          
          
          p1 = Pro * pi(hw1i,[X y1],1,i);
          p2 = Pro * pi(hw1i,[X y2],1,i);         
          
          
          %Ph = Phid([hw1i hzij hxjn]);
          %Po = Pobs(hw1i(1:i-1),X,1,i-1);
          %Po = Po * Pobs(hzij,zij,i+1,j-1);
          %Po = Po * Pobs(hxjn,xjn,j,n);
          %p1 = Ph * Po * B(y1,hw1i(i),i);
          %p2 = Ph * Po * B(y2,hw1i(i),i);          
          
          Sum = Sum + p1/Py1 - p2/Py2;
        end
      end
    end
  end
  h = h + .5*abs(Sum);
end
return