//// Plot netlet formula of Anninos et al.
clear; clf; 
h=0.; c_in=10; c_ex=10;
for theta=1:7;
    rec=0;
    for a=0.:0.01:1.; //fraction of active nodes
        rec=rec+1;
        tmp1=(1-h)*a*c_ex; tmp2=h*a*c_in;
        s1=0; s2=0;
        for n=0:theta-1; s1=s1+tmp1^n/prod(1:n); end;
        for n=0:20;      s2=s2+tmp2^n/prod(1:n); end;
        a0(rec)=a;
        a1(rec)=(1-a)*exp(-tmp2)*s2*(1-exp(-tmp1)*s1);
    end
    plot(a0,a1);
end
plot([0 1],[0 1],'k'); plot([0 1],[1 0],'k')
xlabel('a(t)'); ylabel('a(t+1)');