/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:57 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "isranp.h"
#include <float.h>
#include <stdlib.h>
		/* PARAMETER translations */
#define	M	97
#define	NMAX	84
#define	TWO	2.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
	/* COMMON translations */
struct t_rancs2 {
	float snums[M];
	}	rancs2;
struct t_rancs1 {
	long int sptr;
	LOGICAL32 sgflag;
	}	rancs1;
	/* end of COMMON translations */
long /*FUNCTION*/ isranp(
float xmean)
{
	long int _l0, isranp_v, n;
	static long int last, nmid;
	float s, sprev, term, u;
	static float sum[NMAX-(0)+1], tlast;
	static LOGICAL32 first = TRUE;
	static float elimit = ZERO;
	static float xmsave = 0.e0;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Snums = &rancs2.snums[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 2000-12-01 ISRANP Krogh  Removed unused parameter MONE.
	 *>> 1996-03-30 ISRANP Krogh  Added external statement, rem. F90 comment.
	 *>> 1995-11-16 ISRANP Krogh  Converted from SFTRAN to Fortran 77.
	 *>> 1994-10-19 ISRANP Krogh  Changes to use M77CON
	 *>> 1994-06-24 ISRANP CLL Changed common to use RANC[D/S]1 & RANC[D/S]2.
	 *>> 1994-06-23 CLL Changed name to I[D/S]RANP.
	 *>> 1992-03-16 CLL
	 *>> 1991-11-26 CLL Reorganized common. Using RANCM[A/D/S].
	 *>> 1991-11-22 CLL Added call to RAN0, and SGFLAG in common.
	 *>> 1991-01-15 CLL Reordered common contents for efficiency.
	 *>> 1990-01-23 CLL  Making names in common same in all subprogams.
	 *>> 1987-04-23 IRANP  Lawson  Initial code.
	 *        Returns one pseudorandom integer from the Poisson distribution
	 *     with mean and variance = XMEAN.
	 *        The probability of occurrence of the nonnegative integer k in
	 *     the Poisson distribution with mean XMEAN is */
	/*                P(k) = exp(-XMEAN) * XMEAN**k / k! */
	/*        Let SUM(n) denote the sum for k = 0 to n of P(k).
	 *     Let U be a random sample from a uniform
	 *     distribution on [0.0, 1.0].  The returned value of ISRANP will be
	 *     the smallest integer such that S(ISRANP) .ge. U.
	 *     This variable has mean and variance = XMEAN.
	 *     Reference: Richard H. Snow, Algorithm 342, Comm ACM, V. 11,
	 *     Dec 1968, p. 819.
	 *     Code based on subprogram written for JPL by Stephen L. Ritchie,
	 *     Heliodyne Corp. and Wiley R. Bunton, JPL, 1969.
	 *     Adapted to Fortran 77 for the JPL MATH77 library by C. L. Lawson &
	 *     S. Y. Chiu, JPL, Apr 1987.
	 *     ------------------------------------------------------------------
	 *                    Subprogram argument and result
	 *
	 *     XMEAN [in]  Mean value for the desired distribution.
	 *           Must be positive and not so large that exp(-XMEAN)
	 *           will underflow.  For example, on a machine with underflow
	 *           limit 0.14E-38, XMEAN must not exceed 88.
	 *
	 *     ISRANP [Returned function value]  Will be set to a nonnegative
	 *           integer value if computation is successful.
	 *           Will be set to -1 if XMEAN is out of range.
	 *     ------------------------------------------------------------------
	 *--S replaces "?": I?RANP, ?ERM1, ?ERV1, ?RANUA
	 *--&                 RANC?1, RANC?2, ?PTR, ?NUMS, ?GFLAG
	 *     Generic intrinsic functions referenced: EXP, LOG, MIN, NINT
	 *     Other MATH77 subprogram referenced: RAN0
	 *     Other MATH77 subprograms needed: ERMSG, ERFIN, AMACH
	 *     Common referenced: /RANCS1/, /RANCS2/
	 *     ------------------------------------------------------------------ */
 
 
	/*     ------------------------------------------------------------------ */
	if (first || (xmean != xmsave))
	{
		/*                                      Set for new XMEAN
		 *     R1MACH(1) gives the underflow limit.
		 *     ELIMIT gives a limit for arguments to the exponential function
		 *     such that if XMEAN .le. ELIMIT,  exp(-XMEAN) should not
		 *     underflow.
		 * */
		if (elimit == ZERO)
			elimit = -logf( FLT_MIN*TWO );
		if (xmean <= ZERO || xmean > elimit)
		{
			serm1( "ISRANP", 1, 0, "Require 0 .lt. XMEAN .le. ELIMIT"
			 , "XMEAN", xmean, ',' );
			serv1( "ELIMIT", elimit, '.' );
			isranp_v = -1;
		}
		else
		{
			last = 0;
			tlast = expf( -xmean );
			sum[0] = tlast;
			xmsave = xmean;
			nmid = nint( xmean );
		}
		if (first)
		{
			first = FALSE;
			ran0();
		}
	}
	/*               SET U = UNIFORM RANDOM NUMBER */
	rancs1.sptr -= 1;
	if (rancs1.sptr == 0)
	{
		sranua( rancs2.snums, M );
		rancs1.sptr = M;
	}
	u = Snums[rancs1.sptr];
 
	if (u > sum[last])
	{
		/*                          COMPUTE MORE SUMS */
		n = last;
		sprev = sum[last];
		term = tlast;
L_50:
		;
		n += 1;
		term = term*xmean/(float)( n );
		s = sprev + term;
		if (s == sprev)
		{
			isranp_v = n;
			return( isranp_v );
		}
		if (n <= NMAX)
		{
			last = n;
			sum[last] = s;
			tlast = term;
		}
		if (u <= s)
		{
			isranp_v = n;
			return( isranp_v );
		}
		sprev = s;
		goto L_50;
	}
	else
	{
		/*                            SEARCH THROUGH STORED SUMS
		 *     Here we already know that U .le. SUM(LAST).
		 *     It is most likely that U will be near SUM(NMID).
		 * */
		if (nmid < last)
		{
			if (u > sum[nmid])
			{
				for (n = nmid + 1; n <= last; n++)
				{
					if (u <= sum[n])
					{
						isranp_v = n;
						return( isranp_v );
					}
				}
			}
		}
 
		for (n = min( nmid, last ) - 1; n >= 0; n--)
		{
			if (u > sum[n])
			{
				isranp_v = n + 1;
				return( isranp_v );
			}
		}
		isranp_v = 0;
	}
	return( isranp_v );
} /* end of function */