function [Mvars, Mpos, D] = mapsvarspos(n, P)

Np = factorial(n+P)/factorial(n)/factorial(P)-1;
%Np = prod(n+P:-1:n+1)/factorial(2)-1;

D = zeros(Np+1);
D(1,:) = 1:Np+1;

Mvars = cell(Np+1,1);
Mpos = cell(n,1);

Mvars{1} = [];

vars = 1;
for l=1:Np
    Mvars{l+1} = vars;
    Mpos = addrecur(Mpos, vars, l, n);
    vars = allnxtseq(vars,n);
end;    

for l1=1:Np
    for l2=l1:Np
        v1 = Mvars{l1+1};
        v2 = Mvars{l2+1};
        if length(v1)+length(v2) <= P
            D(l1+1,l2+1) = posrecur(Mpos,sort([v1 v2]))+1;
        end
    end
end;    
D = D+D'-diag(diag(D));


function Mc = addrecur(Mc, inds, val, n)
if length(inds) == 1
    Mc{inds(1)}{1} = val;
else
    if length(Mc{inds(1)}) == 2
        subMc = Mc{inds(1)}{2};
    else
        subMc = cell(n-inds(1));
    end;
    subMc = addrecur(subMc, inds(2:end), val, n);
    Mc{inds(1)}{2} = subMc;
end;    
