/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:56 */
/*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 "drang.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	M	97
#define	ONE	1.0e0
#define	TWO	2.0e0
		/* end of PARAMETER translations */
 
	/* COMMON translations */
struct t_rancd2 {
	double dnums[M];
	}	rancd2;
struct t_rancd1 {
	long int dptr;
	LOGICAL32 dgflag;
	}	rancd1;
	/* end of COMMON translations */
double /*FUNCTION*/ drang()
{
	double drang_v, s, u3, xx, yy;
	static double r, x, y;
	static LOGICAL32 first = TRUE;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Dnums = &rancd2.dnums[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.
	 *>> 1996-04-16 DRANG WVS SQRT(abs(TWO*log(U3))) avoids sqrt(-0.0)
	 *>> 1994-10-20 DRANG Krogh  Changes to use M77CON
	 *>> 1994-06-24 DRANG CLL Changed common to use RANC[D/S]1 & RANC[D/S]2.
	 *>> 1992-03-16 DRANG CLL
	 *>> 1991-11-26 DRANG CLL Reorganized common. Using RANCM[A/D/S].
	 *>> 1991-11-22 DRANG CLL Added call to RAN0, and DGFLAG in common.
	 *>> 1991-01-15 DRANG CLL Reordered common contents for efficiency.
	 *>> 1990-01-23 DRANG CLL Making names in common same in all subprogams.
	 *>> 1987-06-09 DRANG CLLawson  Initial code.
	 *
	 *     Returns one pseudorandom number from the Gausian (Normal)
	 *     distribution with zero mean and unit standard deviation.
	 *     Method taken from Algorithm 334, Comm. ACM, July 1968, p. 498.
	 *     Implemented at JPL in Univac Fortran V in 1969 by Wiley R. Bunton
	 *     of JPL and Stephen L. Richie of Heliodyne Corp.
	 *
	 *     Adapted to Fortran 77 for the MATH 77 library by C. L. Lawson and
	 *     S. Y. Chiu, JPL, April 1987, 6/9/87.
	 *     ------------------------------------------------------------------
	 *--D replaces "?": ?RANG, ?RANUA, RANC?1, RANC?2, ?PTR, ?NUMS, ?GFLAG
	 *     RANCD1 and RANCD2 are common blocks.
	 *     Uses intrinsic functions, log and sqrt.
	 *     Calls DRANUA to obtain an array of uniform random numbers.
	 *     Calls RAN0 to initialize DPTR and DGFLAG.
	 *     ------------------------------------------------------------------
	 *                        Important common variables
	 *
	 *     DPTR [integer] and DGFLAG [logical]
	 *
	 *          Will be set to DPTR = 1 and DGFLAG = .false.
	 *          when RAN0 is called from this subr if this
	 *          is the first call to RAN0.  Also set to these values when
	 *          RAN1 or RANPUT is called by a user.
	 *
	 *          DGFLAG will be set true on return from this subr when the
	 *          algorithm has values that it can save to reduce the amount
	 *          of computation needed the next time this function is
	 *          referenced.  Will be set false in the contrary case.
	 *
	 *          DPTR is the index of the last location used in the
	 *          common array DNUMS().  This index is counted down.
	 *
	 *     DNUMS() [floating point]  Buffer of previously computed uniform
	 *          random numbers.
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (first)
	{
		first = FALSE;
		ran0();
	}
 
	if (!rancd1.dgflag || rancd1.dptr == 1)
	{
 
		/*     Use the Von Neuman rejection method for choosing a random point
		 *     (X,Y) in the unit circle, X**2 + Y**2 .le. 1.0.
		 *     Then the angle Theta = arctan(Y/X) is random, uniform in
		 *     (-Pi/2, Pi/2), and Phi = 2*Theta is random, uniform in (-Pi, Pi).
		 *     Define S = X**2 + Y**2, then
		 *     sin(Theta) = Y/sqrt(S),    cos(Theta) = X/sqrt(S),
		 *     sin(Phi) = 2*X*Y/S,    and cos(Phi) = (X**2 - Y**2)/S.
		 * */
L_10:
		;
		/*                              Set X = random, uniform in [0., 1.] */
		rancd1.dptr -= 1;
		if (rancd1.dptr == 0)
		{
			dranua( rancd2.dnums, M );
			rancd1.dptr = M;
		}
		x = Dnums[rancd1.dptr];
		/*                              Set Y = random, uniform in [-1., 1.] */
		rancd1.dptr -= 1;
		if (rancd1.dptr == 0)
		{
			dranua( rancd2.dnums, M );
			rancd1.dptr = M;
		}
		y = TWO*Dnums[rancd1.dptr] - ONE;
 
		xx = x*x;
		yy = y*y;
		s = xx + yy;
		if (s > ONE)
			goto L_10;
 
		/*     Choose R randomly from Chi-squared distribution and
		 *     normalize with S.
		 *
		 *                              Set U3 = random, uniform in [0., 1.] */
		rancd1.dptr -= 1;
		if (rancd1.dptr == 0)
		{
			dranua( rancd2.dnums, M );
			rancd1.dptr = M;
		}
		u3 = Dnums[rancd1.dptr];
		/*        Changed -TWO*log(U3) to abs(TWO*log(U3)) because Lahey LF90
		 *        2.00 on a pentium produced -0.0 for -TWO*log(1.0), then got a
		 *        floating point exception on sqrt(-0.0). */
		r = sqrt( fabs( TWO*(log( u3 )) ) )/s;
 
		/*                                Compute result as  R*Sin(PHI)
		 * */
		drang_v = (xx - yy)*r;
		rancd1.dgflag = TRUE;
		return( drang_v );
	}
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *        Come here when DGFLAG is true and DPTR is not 1.
	 *
	 *                                Compute result as  R*Cos(PHI) */
	drang_v = TWO*x*y*r;
	rancd1.dgflag = FALSE;
	return( drang_v );
} /* end of function */