function [C12,C34]=FamilySize(J,MAX)
% clear all;close all;[C12,C34]=FamilySize(10^6,25)
D=2509;
q12=[103*231/(D^2) 2/317 26/317 9/D];q12=q12/2;
q34=[103*231/(D^2) 80/317 34/D ];q34=q34/2;X=[1 2 1 0];
X=[1 2 1 0];
C12=zeros(MAX,1);C34=C12;
for N=1:MAX
    N
    C12(N)=Multinomial(N,q12,0)*(J/N);
    C34(N)=Multinomial(N,q34,X)*(J/N);
end;
semilogy(C12);hold on;semilogy(C34,'r')
   