func countincell1d (xx,L=,nl=,nn=)
/* DOCUMENT 2D countincell1d(xx,L)
 EXAMPLE
   uu=[];for(i=1;i<=35;i++){xx=random(80);u=countincell1d(xx,L=1,nn=30,nl=50); grow,uu,[u];}
   plb,uu(,::5,avg),indgen(0:29),(1./indgen(50))(::5),width=4
*/
{
  res=[];
  if (is_void(nn)) nn=1+numberof(xx);
  if (is_void(nl)) nl=15;
  if (is_void(L)) L=1.;
  for(i=1;i<=nl;i++)
    {
      rg=indgen(0:i)*double(L)/i;
      u=histo1d(xx,rg);
      //           error,"ah";
      grow,res, [histo1d(u(*),indgen(0:nn))];
    }
  res=double(res);
  norm=res(sum,); w= where(!norm); if(numberof(w)) norm(w)=1;
  return res/norm(-,);
}


func countincell2d (xx,L=,nl=,nn=)
/* DOCUMENT 2D countincell2d(xx,L)
   EXAMPLE
   uu=[];for(i=1;i<=15;i++){xx=random(8000,2);u=countincell2d(xx,L=1,nn=30,nl=50); grow,uu,[u];} pli,uu(,,avg)
 */
{
    res=[];
  if (is_void(nn)) nn=1+numberof(xx);
  if (is_void(nl)) nl=15;
  if (is_void(L)) L=1.;
  for(i=1;i<=nl;i++)
    {
      rg=indgen(0:i)*double(L)/i;
      u=histo2d(xx,rg,rg);
      //      error,"ah";
      grow,res, [histo1d(u(*),indgen(0:nn))];
    }
  res=double(res);
  return res/res(sum,)(-,);
}



func countincell3d (xx,L=,nl=,nn=)
/* DOCUMENT 3D countincell
   EXAMPLE   uu=[];for(i=1;i<=15;i++){xx=random(8000,3);u=countincell3d(xx,L=1,nn=30,nl=50); grow,uu,[u];} pli,uu(,,avg)
 */
{
    res=[];
  if (is_void(nn)) nn=1+numberof(xx);
  if (is_void(nl)) nl=15;
  if (is_void(L)) L=1.;
  for(i=1;i<=nl;i++)
    {
      rg=indgen(0:i)*double(L)/i;
      u=histo3d(xx,rg,rg,rg);
      //      error,"ah";
      grow,res, [histo1d(u(*),indgen(0:nn))];
    }
  res=double(res);
  return res/res(sum,)(-,);
}



func binomialPDF (x,p,n)
/* DOCUMENT 
     PDF of binomial ditribution
   SEE ALSO:
 */
{
if(numberof(x)>1)
  {
    qq=[];
    for(i=1;i<=numberof(x);i++)
      grow,qq,binomialPDF(x(i),p,n);
    return qq;
  }

  return bico(n,x)* p^x*(1-p)^(n-x);

}

func binomialCDF (x,p,n)
/* DOCUMENT 
     CDF of binomial ditribution
   SEE ALSO:
 */
{
  return binomialPDF(indgen(0:x),p,n)(sum);
}


func binomialMoments (q,p,n)
/* DOCUMENT 
     Moments of binomial ditribution
   SEE ALSO:
 */
{
  if(numberof(q)>1)
  {
    qq=[];
    for(i=1;i<=numberof(q);i++)
      grow,qq,binomialMoments(q(i),p,n);
    return qq;
  }
  qq=1;
  for(i=0;i<=q-1;i++) qq *= p*(n-i);
  return qq;
}




func count2moment (counts,&err,k=,error=)
/* DOCUMENT
   returns the normalized factorial moments
   F_k= < N(N-1)..(N-k+1) >/< N(N-1)..(N-k+1) >_binomial
   note that the first element corresponds to the void 
   EXAMPLE
   uu=[];for(i=1;i<=1000;i++){xx=random(30);u=countincell1d(xx,L=1,nn=30,nl=25); grow,uu,[u];}
   tt=count2moment(uu,er,k=4,error=1);
   plb,tt(2:20,2:4)-1,marker=1,width=4;
   plb,(tt+er/sqrt(1000))(2:20,2:4)-1,type=2
   plb,(tt-er/sqrt(1000))(2:20,2:4)-1,type=3
   SEE ALSO:
 */
{
 
  if(is_void(k)) k=1;
  d=dimsof(counts);
  if(d(1)==2) { counts=[counts];   d=dimsof(counts);}
  p=d(0);
  n=d(2)-1;
  nl=d(3);
  res=[];
  if(error) err=[];
  vec=indgen(0:n);
  for(i=0;i<=k;i++)
    {
      grow,res,[(counts*vec(,-,-))(sum,,avg)];
       if(error) grow,err,[(counts*vec(,-,-))(sum,,rms)];
      vec *= max(0,indgen(0:n)-i-1);
    }
  res=double(res);
  for(i=1;i<=nl;i++)
    {
      tt=binomialMoments(indgen(k+1),1./i,n);
      w=where(tt==0); if(numberof(w)) tt(w)=1;
      res(i,) /= tt;
       if(error) err(i,) /= tt;
    }
  
  return res;
}


/* pp=1000;nn=15;
x=random(pp);y=random(pp); z=random(pp)
u=histo2d([x,y],span(0,1,nn),span(0,1,nn));
pli,u
*/