%mkQ.m  makes RM as rate-matrix for birth-death MC. 
%	computes Q(KK) as transition-prob matrix

RM=0;Q=0;

RM=zeros(Sp+1,Sp+1);

%	rate matrix RM first row
RM(1,1)=-Sp*lam;RM(1,2)=Sp*lam;

%	rate matrix RM last row
RM(Sp+1,Sp)=Sp*theta;RM(Sp+1,Sp+1)=-Sp*theta;

%	all other rows
for l=2:Sp RM(l,l-1)=(l-1)*theta;
	RM(l,l)=-((Sp-l+1)*lam+(l-1)*theta);
	RM(l,l+1)=(Sp-l+1)*lam;
end

%evaluate Q(KK) as matrix-exp of (RM.*KK) for transition probs
%Cox & Miller THY STOC PROC 1965, p182 where it's called P(t)
%	and P(t) = expm(Qt) for constant rate matrix Q
		Q=expm(RM.*KK);

%compute q as long-term occupancies
q=0;q(1)=1;sum=1;
for j=2:(Sp+1) q(j) = q(j-1)*RM(j-1,j)/RM(j,j-1);
		sum=sum+q(j);
end
q=q./sum;