Changeset 3364 in Sophya for trunk/SophyaLib/BaseTools

Timestamp:

Oct 29, 2007, 6:34:27 PM (18 years ago)

Author:

cmv

Message:

add PoissRandLimit pour tirage sur distrib poisson avec protection, cmv 29/10/2007

Location:

trunk/SophyaLib/BaseTools

Files:

• srandgen.c (modified) (1 diff)
• srandgen.h (modified) (1 diff)

trunk/SophyaLib/BaseTools/srandgen.c

@@ -254 +254 @@
   double pp,ppi,x;
   int n;
+  ppi = pp = exp(-mu);
+  x = drand01();
+  n = 0;
+  while (x > ppi) {
+    n++;
+    pp = mu*pp/(double)n;
+    ppi += pp;
+  }
+  return n;
+}
+
+/*! \ingroup  BaseTools
+  \brief Poisson random number generator with approximation for large mean.
+  
+  Return an integer value (>=0) corresponding a Poisson distribution with mean \b mu.
+  If \b mu is greater or equal than \b mumax,
+  the large number gaussian approximation is applied.
+  \b mumax=10 is a good value.
+*/
+unsigned long PoissRandLimit(double mu,double mumax)
+// p(n) = (mu^n / n!) exp(-mu)  for n = 0, 1, 2, ...
+{
+  double pp,ppi,x;
+  unsigned long n;
+
+  if(mu>=mumax) {
+    pp = sqrt(mu);
+    while( (x=pp*NorRand()) < -mu );
+    return (unsigned long)(mu+x+0.5);
+  }
+
   ppi = pp = exp(-mu);
   x = drand01();

trunk/SophyaLib/BaseTools/srandgen.h

@@ -74 +74 @@
 
 int PoissRand(double mu);
+unsigned long PoissRandLimit(double mu,double mumax /* =10. */ );
 
 TIREALEA *init_tirage_alea(int nbin,double xmin,double xmax,double (*fonc) (double));