/*
#pragma module g05caf "C1.00"	don't include if doing plus-list compiles
*/
/*
   Bart Z. Lederman	06-Jul-2006 Replacement for the FORTRAN
				    routine.

   I won't guarantee that this does exactly what the old routine did in
   every possible case, but it does return random numbers, and icalcv
   works with it; which is all I expect from it.
*/

/*
==========================================================================
			random number generator #2
==========================================================================
*/

/*
   Replacements for the FORTRAN functions
*/

#ifndef MIN
#define MIN(x,y)        ((x) < (y) ? (x) : (y))
#endif
#ifndef MAX
#define MAX(x,y)        ((x) > (y) ? (x) : (y))
#endif

/*
   These were variables in FORTRAN, can be declared as constants here.
*/

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0 / IM1)
#define IMM1 (IM1 - 1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define EPS (1.E-7)
#define RNMX (1. - EPS)

/*    int idum2iv[NTAB], iy;						*/
/*	  SAVE            iv,iy,idum2					*/
/*	  DATA            idum2/123456789/,iv/NTAB*0/,iy/0/		*/

    int idum2 = 123456789;
    int iy = 0;
    int iv[NTAB] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
		    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
/*
   This one apparently needs to be calculated.
*/
    int NDIV = (1 + IMM1 / NTAB);

double g05caf (int tidum)
{
    int idum, j, k;
/*
   This modified version of RAN2 (2nde edition) use an intermediate
   number tidum (the r=function parameter) because g05caf is often
   called as g05caf(0). The test is also .LT. and not .LE. for the same
   reason.  Therefore, only a negative seed value will initialize the
   random number generator. The program will still crash on a call such
   as g05caf(-2000) since it expect g05caf(int *tidum).
*/
    if (tidum < 0)
    {
	idum  = tidum;
	tidum = 0;
    };

    if (idum < 0)
    {
	idum = MAX (-idum, 1);
	idum2 = idum;

	for (j = NTAB + 8; j <= 1; j--)	/* may need to adjust end point	*/
	{
	    k = idum / IQ1;
	    idum = IA1 * (idum - k * IQ1) - k * IR1;
	    if (idum < 0) idum = idum + IM1;
	    if (j <= NTAB) iv[j] = idum;
	};

	iy = iv[1];
    };

    k = idum / IQ1;

    idum = IA1 * (idum - k * IQ1) - k * IR1;

    if (idum < 0) idum = idum + IM1;

    k = idum2 / IQ2;

    idum2 = IA2 * (idum2 - k * IQ2) - k * IR2;

    if (idum2 < 0) idum2 = idum2 + IM2;

/*    j = 1 + iy / NDIV;    */
    j = 1 + iy;
    j = j / NDIV;

    iy = iv[j] - idum2;
    iv[j] = idum;

    if (iy < 1) iy = iy + IMM1;

    return (MIN (AM * iy, RNMX));
}