/* @(#)randev.c	17.1.1.1 (ES0-DMD) 01/25/02 17:34:55 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT   (c)  1995  European Soutern Observatory
.IDENT       randev.c
.LANGUAGE    C
.AUTHOR      P.Grosbol,  IPG/ESO
.ENVIRON     UNIX
.KEYWORDS    Random numbers, exponential, gaussian
.COMMENT     Algorithm is adapted from 'Numerical Recipes in C' p.214
.VERSION     1.0  1990-Nov-07 : Creation,  PJG
.VERSION     1.1  1994-May-28 : Change name of constant PI, PJG
.VERSION     1.2  1995-Mar-09 : Correct error + double gammln, PJG
-----------------------------------------------------------------------*/
#include <math.h>                    /* Standard mathematical library  */

#define PICONST       3.141592654

double expdev(pseed)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE   compute random number taken from an Exponential disribution
           with unit mean.
.RETURN    exponential distributed random number
-----------------------------------------------------------------------*/
int    *pseed;                /* pointer to seed if <0 initiate seq.   */
{
  double  x;
  float   randm();

  x = randm(pseed);
  return -log(x);
}

static  int     gflag = 0;
static  double  gval;

double gaudev(pseed)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE   compute random number taken from a Gaussian disribution
           with zero mean and unit standard deviation.
.RETURN    gaussian distributed random number
-----------------------------------------------------------------------*/
int    *pseed;                /* pointer to seed if <0 initiate seq.   */
{
  double  r, fac, v1, v2;
  float   randm();

  if (gflag) { gflag = 0; return gval; }
  else {
     do {
	v1 = 2.0*randm(pseed) - 1.0;
	v2 = 2.0*randm(pseed) - 1.0;
	r  = v1*v1 + v2*v2;
      } while (1.0<=r);
     fac  = sqrt( -2.0*log(r)/r );
     gval = v1*fac;
     gflag = 1;
     return v2*fac;
   }
}

double gamdev(ia,pseed)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE   compute random number taken from a Gamma  disribution of
           interger order of 'ia'
.RETURN    gamma distributed random number
-----------------------------------------------------------------------*/
int    ia;                    /* integer order of gamma function       */
int    *pseed;                /* pointer to seed if <0 initiate seq.   */
{
  int   j;
  float randm();
  double am, e, s, v1, v2, x, y;

  if (ia < 1) return 0.0;
  if (ia < 6) {
     x = 1.0;
     for (j=0; j<ia; j++) x *= randm(pseed);
     x = -log(x);
   } 
  else {
     do {
	do {
	   do {
	      v1 = 2.0*randm(pseed) - 1.0;
	      v2 = 2.0*randm(pseed) - 1.0;
	    } while (v1*v1+v2*v2 > 1.0);
	   y  = v2/v1;
	   am = ia - 1;
	   s  = sqrt(2.0*am+1.0);
	   x  = s*y + am;
	 } while (x <= 0.0);
	e = (1.0+y*y)*exp(am*log(x/am)-s*y);
      } while (randm(pseed) > e);
   }

  return x;
}

static double sq, alxm, g, oldm=(-1.0);

double poidev(xm,pseed)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE   compute random number taken from a Poisson disribution
           with mean 'xm'.
.RETURN    poisson distributed random number
-----------------------------------------------------------------------*/
double     xm;                /* mean value of poisson distribution    */
int        *pseed;            /* pointer to seed if <0 initiate seq.   */
{
  float    randm();
  double   em, t, y, gammln();
  
  if (xm < 12.0) {
     if (xm != oldm) {
	oldm = xm;
	g = exp(-xm);
      }
     em = -1;
     t = 1.0;
     do {
	em += 1.0;
	t  *= randm(pseed);
      } while (t > g);
   } else {
      if (xm != oldm) {
	 oldm = xm;
	 sq = sqrt(2.0*xm);
	 alxm = log(xm);
	 g = xm*alxm - gammln(xm+1.0);
       }
      do {
	 do {
	    y = tan(PICONST*randm(pseed));
	    em = sq*y + xm;
	  } while (em < 0.0);
	 em = floor(em);
	 t = 0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
       } while (randm(pseed) > t);
    }
  return em;
}