wolffd@0: function pot = extension_pot(oldpot, ddom_u, dsizes, ctaildom_u, csizes)
wolffd@0: % EXTENSION_POT Extense a stable CG potential.
wolffd@0: % pot = extension_pot(oldpot, ddom_u, ctaildom_u, dsizes, csizes)
wolffd@0: % ddom_u Added discrete nodes
wolffd@0: % ctaildom_u Added continuous tail nodes
wolffd@0: % csizes is the size of the tail nodes.
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % A CG potential can be extended by adding discrete variables to its %
wolffd@0: % domain of continuous variables to its tail                         %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: ddom = myunion(oldpot.ddom, ddom_u);
wolffd@0: ctaildom = myunion(oldpot.ctaildom, ctaildom_u);
wolffd@0: cheaddom = oldpot.cheaddom;
wolffd@0: udom = myunion(ddom_u, ctaildom_u);
wolffd@0: domain = myunion(oldpot.domain, udom);
wolffd@0: 
wolffd@0: ns = zeros(1,max(domain));
wolffd@0: ns(ddom_u) = dsizes;
wolffd@0: ns(ctaildom_u) = csizes;
wolffd@0: ns(oldpot.ddom) = oldpot.dsizes;
wolffd@0: ns(oldpot.cheaddom) = oldpot.cheadsizes;
wolffd@0: ns(oldpot.ctaildom) = oldpot.ctailsizes;
wolffd@0: 
wolffd@0: dsizes = ns(ddom);
wolffd@0: dsize = prod(ns(ddom));
wolffd@0: cheadsizes = ns(cheaddom);
wolffd@0: cheadsize = sum(ns(cheaddom));
wolffd@0: ctailsizes = ns(ctaildom);
wolffd@0: ctailsize = sum(ns(ctaildom));
wolffd@0: 
wolffd@0: BZ = zeros(cheadsize, ctailsize);
wolffd@0: potarray = cell(1, dsize);
wolffd@0: mask = find_equiv_posns(oldpot.ddom, ddom);
wolffd@0: 
wolffd@0: tmask = find_equiv_posns(oldpot.ctaildom, ctaildom);
wolffd@0: tu = block(tmask, ctailsizes);
wolffd@0: 
wolffd@0: for i=1:dsize
wolffd@0:     sub1 = ind2subv(dsizes, i);
wolffd@0:     sub2 = sub1(mask);
wolffd@0:     ind = subv2ind(oldpot.dsizes, sub2);
wolffd@0:     if isempty(ind)
wolffd@0:         ind = 1;
wolffd@0:     end
wolffd@0:     potc = struct(oldpot.scgpotc{ind});
wolffd@0:     p = potc.p;
wolffd@0:     B = BZ;
wolffd@0:     if ~isempty(B)
wolffd@0:         B(:, tu) = potc.B;
wolffd@0:     end
wolffd@0:     potarray{i} = scgcpot(cheadsize, ctailsize, p, potc.A, B, potc.C);
wolffd@0: end
wolffd@0: 
wolffd@0: pot = scgpot(ddom, cheaddom, ctaildom, ns,potarray);