function Z = Tprod_c(X,Y,sx,sy)
% a custom - defined tensor product

DX = size(X);
DY = size(Y);
[rX,sX,rx] = decompDX(DX,sx);
[rY,sY,ry] = decompDX(DY,sy);
Ds = DX(sx);
lds = length(Ds);
DZ = [rX,rY];
ldx = length(DX);
ldy = length(DY);
lrx = ldx-lds;
lry = ldy-lds;
ldz = lrx + lry;
%lrx = length(rX);
%lry = length(rY);

%NX = prod(DX);
%NY = prod(DY);
NZ = prod(DZ);
Ns = prod(Ds); %this will be 1 if Ds is empty


Z = zeros(NZ,1);
X = X(:);
Y = Y(:);


Cx = zeros(1,ldx);
Cy = zeros(1,ldy);

for iz = 1:NZ
  Cz = ind2coord(iz,DZ);
  Crx = Cz(1:lrx);
  Cry = Cz(lrx+1:lrx+lry);
  Cx(rx) = Crx;
  Cy(ry) = Cry;
  for is = 1:Ns
    Cs = ind2coord(is,Ds);
    Cx(sx) = Cs;
    ix = coord2ind(Cx,DX);
    Cy(sy) = Cs;
    iy = coord2ind(Cy,DY);
    Z(iz) = Z(iz) + X(ix) * Y(iy);
  end
end

Z = reshape(Z,DZ);

return



function [rX,sX,iy] = decompDX(DX,ix)
sX = DX(ix);
DX(ix) = zeros(1,length(ix));
iy = find(DX);
rX = DX(iy);
return

function c = ind2coord(ind,dim)
ind = ind-1;
c = zeros(size(dim));
for ij = 1:length(dim)
  c(ij) = mod(ind,dim(ij));
  ind = (ind - c(ij))/dim(ij);
end
c = c + 1;
return


function ind = coord2ind(c,dim)
ind = 0;
c = c - 1;
f = 1;
for ij = 1:length(dim)
  ind = ind + c(ij)*f;
  f = f*dim(ij);
end
ind = ind+1;
return


function junk1
%function DZ=DXYZ(DX,DY,ix,iy)
DZ = apdimZ(DX,ix);
DZ = [DZ,apdimZ(DY,iy)];
%return
%function DZ = apdimZ(DX,ix)
DZ = DX(1:ix(1)-1);
for k=1:length(ix)-1
  DZ = [DZ,DX(ix(k)+1:ix(k+1)-1)];
end
DZ = [DZ,DX(ix(end)+1:end)];
return