#ifndef _UTILS1_H_
#define _UTILS1_H_

#define PI (2.0*asin(1.0))
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define MIN(a,b) (a < b ? a : b)
#define MAX(a,b) (a > b ? a : b)
#define NR_END 1
#define FREE_ARG char*

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

/*--------------------------------------------------------------------------*/
/*    UNIFORM DEVIATES                                                      */
/*--------------------------------------------------------------------------*/

double ran2(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;

	if (*idum <= 0 || !iy) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if (*idum < 0) *idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = *idum;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}


#define M1 	259200L
#define IA1 	7141L
#define IC1 	54773L
#define RM1	(1.0/M1)
#define M2	134456L
#define IA2	8121L
#define IC2	28411L
#define RM2 	(1.0/M2)
#define M3	243000L
#define IA3	4561L
#define IC3	51349L

double ran3(int *idum)
{
 static long ix1, ix2, ix3;
 static double r[98];
 double temp;
 static int iff=0;
 int j;

 if(*idum < 0 || iff == 0)
 {
  iff = 1;
  ix1 = (IC1 - (*idum)) % M1;
  ix1 = (IA1 * ix1 + IC1) % M1;
  ix2 = ix1 % M2;
  ix1 = (IA1 * ix1 + IC1) % M1;
  ix3 = ix1 % M3;
  for (j=1; j<=97; j++)
  {
   ix1 = (IA1 * ix1 + IC1) % M1;
   ix2 = (IA2 * ix2 + IC2) % M2;
   r[j] = (ix1 + ix2 * RM2) * RM1;
  }
  *idum = 1;
 }

 ix1 = (IA1 * ix1 + IC1) % M1;
 ix2 = (IA2 * ix2 + IC2) % M2;
 ix3 = (IA3 * ix3 + IC3) % M3;

 j = (int)(1 + ((97 * ix3)/M3));
 if (j > 97 || j < 1)
 {
  printf("ERROR: This cannot happen in ran3()\n"); exit(0);
 }

 temp = r[j];
 r[j] = (ix1 + ix2 * RM2) * RM1;

 return temp;

}



double gaussian(double mean,double sd,double r1,double r2)
{
 double next_int;

 next_int = mean + (sd * cos(2.0*PI*r1) * sqrt(-2*log(r2)));
 return next_int;
}

#endif // #ifndef _UTILS1_H_