/*

This file contains all Numerical Recipes that I use.

*/


/********* Begin of Numerical Recipes routines **********/


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


/*
 *  numerical recipes random gaussian number generator
 */

float gasdev(idum)
     int *idum;
{
  static int iset=0;
  static float gset;
  float fac,r,v1,v2;
  float ran1();
  
  if  (iset == 0) {
    do {
      v1=2.0*ran1(idum)-1.0;
      v2=2.0*ran1(idum)-1.0;
      r=v1*v1+v2*v2;
    } while (r >= 1.0);
    fac=sqrt(-2.0*log(r)/r);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  } else {
    iset=0;
    return gset;
  }
}

/*
 *	gaussian - return a gaussian random number of
 *		   variance 1 and mean 0
 */

float gaussian ()
{
	int seed = 0;

	return ( gasdev(&seed) );
}


/*
 *  numerical recipes random number generator
 */

float ran1(idum)
     int *idum;
{
  static long ix1,ix2,ix3;
  static float r[98];
  float temp;
  static int iff=0;
  int j;
  void nrerror();
  
  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=1 + ((97*ix3)/M3);
  if (j > 97 || j < 1) nrerror("RAN1: This cannot happen.");
  temp=r[j];
  r[j]=(ix1+ix2*RM2)*RM1;
  return temp;
}


/*
 *  numerical recipes error routine
 */

void nrerror(error_text)
     char error_text[];
{
  /*	void exit();
   */
  fprintf(stderr,"Numerical Recipes run-time error...\n");
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"...now exiting to system...\n");
  exit(1);
}

#undef M1
#undef IA1
#undef IC1
#undef RM1
#undef M2
#undef IA2
#undef IC2
#undef RM2
#undef M3
#undef IA3
#undef IC3

/********* End of Numerical Recipes routines **********/