/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */
/*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 "srang.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_rancs2 {
	float snums[M];
	}	rancs2;
struct t_rancs1 {
	long int sptr;
	LOGICAL32 sgflag;
	}	rancs1;
	/* end of COMMON translations */
float /*FUNCTION*/ srang()
{
	float s, srang_v, u3, xx, yy;
	static float r, x, y;
	static LOGICAL32 first = TRUE;
		/* 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.
	 *>> 1996-04-16 SRANG WVS SQRT(abs(TWO*log(U3))) avoids sqrt(-0.0)
	 *>> 1994-10-20 SRANG Krogh  Changes to use M77CON
	 *>> 1994-06-24 SRANG CLL Changed common to use RANC[D/S]1 & RANC[D/S]2.
	 *>> 1992-03-16 SRANG CLL
	 *>> 1991-11-26 SRANG CLL Reorganized common. Using RANCM[A/D/S].
	 *>> 1991-11-22 SRANG CLL Added call to RAN0, and SGFLAG in common.
	 *>> 1991-01-15 SRANG CLL Reordered common contents for efficiency.
	 *>> 1990-01-23 SRANG CLL Making names in common same in all subprogams.
	 *>> 1987-06-09 SRANG 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.
	 *     ------------------------------------------------------------------
	 *--S replaces "?": ?RANG, ?RANUA, RANC?1, RANC?2, ?PTR, ?NUMS, ?GFLAG
	 *     RANCS1 and RANCS2 are common blocks.
	 *     Uses intrinsic functions, log and sqrt.
	 *     Calls SRANUA to obtain an array of uniform random numbers.
	 *     Calls RAN0 to initialize SPTR and SGFLAG.
	 *     ------------------------------------------------------------------
	 *                        Important common variables
	 *
	 *     SPTR [integer] and SGFLAG [logical]
	 *
	 *          Will be set to SPTR = 1 and SGFLAG = .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.
	 *
	 *          SGFLAG 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.
	 *
	 *          SPTR is the index of the last location used in the
	 *          common array SNUMS().  This index is counted down.
	 *
	 *     SNUMS() [floating point]  Buffer of previously computed uniform
	 *          random numbers.
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (first)
	{
		first = FALSE;
		ran0();
	}
 
	if (!rancs1.sgflag || rancs1.sptr == 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.] */
		rancs1.sptr -= 1;
		if (rancs1.sptr == 0)
		{
			sranua( rancs2.snums, M );
			rancs1.sptr = M;
		}
		x = Snums[rancs1.sptr];
		/*                              Set Y = random, uniform in [-1., 1.] */
		rancs1.sptr -= 1;
		if (rancs1.sptr == 0)
		{
			sranua( rancs2.snums, M );
			rancs1.sptr = M;
		}
		y = TWO*Snums[rancs1.sptr] - 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.] */
		rancs1.sptr -= 1;
		if (rancs1.sptr == 0)
		{
			sranua( rancs2.snums, M );
			rancs1.sptr = M;
		}
		u3 = Snums[rancs1.sptr];
		/*        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 = sqrtf( fabsf( TWO*(logf( u3 )) ) )/s;
 
		/*                                Compute result as  R*Sin(PHI)
		 * */
		srang_v = (xx - yy)*r;
		rancs1.sgflag = TRUE;
		return( srang_v );
	}
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *        Come here when SGFLAG is true and SPTR is not 1.
	 *
	 *                                Compute result as  R*Cos(PHI) */
	srang_v = TWO*x*y*r;
	rancs1.sgflag = FALSE;
	return( srang_v );
} /* end of function */