function [X, P1] = MEBPexpd(X, P1, samplingftn, weightsftn)

% function X = MEBPexpd(X, samplingftn, weightsftn)
%
% adds new levels to the state of an MEBP process
%
% X in is the current state of the process at step k (levels 0 to nmax)
% P1 in is the vector of weights for the first crossing at each level 
% samplingftn is a function handle for a function that samples from 
%   the offspring sampling distribution n*p(n)/mu
% weightsftn is a function handle for a function that samples from 
%   the weights distribution, given the number of offspring
% X out is the expanded state of the process at step k (nmax may have increased)
% P1 out is the vector of weights for the first crossing at each level
%
% row n+1 of X is (kappa(0,n,k), S(0,n,k), Z(0,n+1,k), alpha(0,n,k), P(0,n+1,k))
% row nmax+1 of X out satisfies S(0,nmax,k) < Z(0,nmax+1,k)
%
% types are 0+ 0- 1+ 1- -1+ -1-, labelled 1 2 3 4 5 6
%
% Owen Jones 6 July 2003, 17 January 2007

nmax = length(X) - 1;
while X{nmax+1}{2} == X{nmax+1}{3}
    nmax = nmax + 1;
    X{nmax+1} = {1, 0, 0, 0, []};
    X{nmax+1}{3} = feval(samplingftn, []);
    X{nmax+1}{2} = ceil(rand*X{nmax+1}{3});
    X{nmax+1}{5} = feval(weightsftn, X{nmax+1}{3});
    if X{nmax}{4} == 3
        if X{nmax+1}{2} == X{nmax+1}{3}
            X{nmax+1}{4} = 3;
        elseif mod(X{nmax+1}{2},2) == 1
            X{nmax+1}{4} = 1;
        else
            X{nmax+1}{4} = 5;
        end
    else
        if X{nmax+1}{2} == X{nmax+1}{3}
            X{nmax+1}{4} = 6;
        elseif mod(X{nmax+1}{2},2) == 1
            X{nmax+1}{4} = 2;
        else
            X{nmax+1}{4} = 4;
        end
    end
    P1(nmax+1) = X{nmax+1}{5}(1);
end