void PIC::inject_beam(double boundary, double vel, double angle, double ndensity, double Tempxy, double Tempz, double dt, long *d1)
{
 double Ninj=vel*ndensity*PI*pow(boundary,2)*dt/macroN; 
 double gamma=1.0/sqrt(1.0-pow(vel/clight,2));
 double utherm=clight*sqrt(pow(charge_p*Tempxy/(mass_p*clight*clight)+1.0,2)-1.0);
 double r_amp, phase;
 double u, v, w;
 for(int j=0; j<Ninj; j++) 
   { 
	Particle part_inj;
	r_amp=boundary*sqrt(ran1(d1));
	phase=2.0*PI*ran1(d1);	
	part_inj.x=r_amp*cos(phase);
	part_inj.y=r_amp*sin(phase);
	part_inj.z=-0.5*length;
	
	do { u = 2.0 * ran1(d1) - 1.0 ;  
	     v = 2.0 * ran1(d1) - 1.0 ;              
	     w = u*u + v*v ; }
	     while (w >= 1.0) ;
	
	part_inj.ux = part_inj.x/r_amp*vel*angle + utherm * u * sqrt(-2.0*log(w)/w) ;
	part_inj.uy = part_inj.y/r_amp*vel*angle + utherm * v * sqrt(-2.0*log(w)/w) ;
	part_inj.uz = gamma*vel;
    part.push_back(part_inj);
   }	
}

float poidev(float xm, long int *idum)
{
 static float sq,alxm,g,oldm=(-1.0);
 float em,t,y;
 
 if (xm < 12.0) {
    if (xm != oldm) {
       oldm=xm;
       g=exp(-xm);
    }
    em = -1;
    t=1.0;
    do {
       ++em;
       t *= ran1(idum);
    } while (t > g);
 } else {
    if (xm != oldm) {
       oldm=xm;
       sq=sqrt(2.0*xm);
       alxm=log(xm);
       g=xm*alxm-gammln(xm+1.0);
    }
    do {
       do {
          y=tan(PI*ran1(idum));
	  em=sq*y+xm;
       } while (em < 0.0);
       em=floor(em);
       t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
    } while (ran1(idum) > t);
 }

 return em;
}

/* ********************************************************/
double gasdev(long int *idum) {
/*
 *  Gaussian randon number generator. Uses ran1.
***********************************************************/
    double ran1(long int *idum);
    static int iset = 0;
    static double gset;
    double fac, rsq, v1, v2;

    if (*idum < 0) iset = 0;
    if (iset == 0) {
        do {
            v1 = 2.0 * ran1(idum) - 1.0;
            v2 = 2.0 * ran1(idum) - 1.0;
            rsq = v1 * v1 + v2 * v2;
        } while (rsq >= 1.0 || rsq == 0.0);
        fac = sqrt(-2.0 * log(rsq) / rsq);
        gset = v1 * fac;
        iset = 1;
        return v2 * fac;
    } else {
        iset = 0;
        return gset;
    }
}

void gen_Ak(double *Ak, int kx, int ky, int kz, double Pk, int dots)
{
    int nkx, nky, nkz;
    double phi, rnd1, rnd2, r;
    
    
    rnd1 = ran1(&iiseed);
    phi = 2 * PI * rnd1;                /* phase */

    rnd2 = ran1(&iiseed);
    r = sqrt(-2 * log(rnd2)*Pk);         /* amplitude */


    Ak[2*(kx*dots*dots + ky*dots + kz)] = r*cos(phi); /* real part */
    Ak[2*(kx*dots*dots + ky*dots + kz)+1] = r*sin(phi); /* imag. part */
    

    if((kx != 0) && (ky != 0) && (kz != 0))
      {
	  nkx = dots - kx;
	  nky = dots - ky;
	  nkz = dots - kz;

	  Ak[2*(nkx*dots*dots + nky*dots + nkz)] = r*cos(phi); /* real part */
	  Ak[2*(nkx*dots*dots + nky*dots + nkz)+1] = -r*sin(phi); /* imag. part */
      }
}

void PIC::set_const_density_xy(double boundary, double T0, long *d1)
{ 
 double u,v,w;	
 double ut=clight*sqrt(pow(fabs(charge_p)*T0/(mass_p*clight*clight)+1.0,2)-1.0);
 for (int j=0 ; j < part.size(); j++)
  {
   	do {
        u = boundary*(2.0*ran1(d1)-1.0); 
     	v = boundary*(2.0*ran1(d1)-1.0);
	   } while (sqrt(pow(u,2)+pow(v,2)) > boundary);
    		
	part[j].x=u;  	
	part[j].y=v;
	part[j].z=length * ran1(d1);
	
	do { u = 2.0 * ran1(d1) - 1.0 ;  
	     v = 2.0 * ran1(d1) - 1.0 ;              
	     w = u*u + v*v ; }
	     while (w >= 1.0) ;

		part[j].ux = ut * u * sqrt(-2.0*log(w)/w) ;
	    part[j].uy = ut * v * sqrt(-2.0*log(w)/w) ;
		part[j].uz = 0.0; 
	 
  }	
}

void Pic::waterbag_xy(double emittance_x, double emittance_y, double alpha_x, double alpha_y,
		double beta_x, double beta_y, double D0, double Ds0, double x0,
		double xs0, double y0, double ys0, int size, long *d){
  emittance_x *= 6.;
  emittance_y *= 6.;
  long j;
  double x, y, xs, ys;
  double xmax = sqrt(emittance_x*beta_x);
  double xsmax = sqrt(emittance_x * (1.0+alpha_x*alpha_x) / beta_x);
  double ymax = sqrt(emittance_y*beta_y);
  double ysmax = sqrt(emittance_y * (1.0+alpha_y*alpha_y) / beta_y);

  for(j=size; j<pics.size(); ++j){
    do{
      x = xmax*(2.0*ran1(d)-1.0);
      y = ymax*(2.0*ran1(d)-1.0);
      xs = xsmax*(2.0*ran1(d)-1.0);
      ys = ysmax*(2.0*ran1(d)-1.0);
  }while((x*x/beta_x + beta_x * pow(xs+alpha_x/beta_x*x, 2)) / emittance_x
	   + (y*y/beta_y + beta_y * pow(ys+alpha_y/beta_y*y, 2)) / emittance_y > 1.0);
    pics[j].x = x + D0*pics[j].dp + x0;  // stuff for head tail modes removed; SP
    pics[j].xs = xs + Ds0*pics[j].dp + xs0;  // stuff for head tail modes removed; SP
    pics[j].y = y+y0;                        // offset in vertrical space SA
    pics[j].ys = ys+ys0;
  }
  init_old_coord();
}

float gasdev(long *idum)
/* Returns a normally distributed variate with zero mean and
a unit variance, using ran1(idum) as the source of uniform deviates */
{
	float ran1(long *idum);
	static int iset=0;
	static float gset;
	float fac, rsq,v1,v2;

	if (*idum < 0) iset = 0;
	if (iset == 0) {
		do {
			v1 = 2.0*ran1(idum)-1.0;
			v2 = 2.0*ran1(idum)-1.0;
			rsq = v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac = sqrt(-2.0*log(rsq)/rsq);
		/*Now amke the Box-Muller transformation to get 
                  two normal deviates. Return one and save the other
		  for next time */
		gset = v1*fac;
		iset = 1;
		return (v2*fac);
	} else {
		iset = 0;
		return (gset);
	}
}

/*-------------------*
 *  Function gasdev  |
 *-------------------*
 | a slight modification of crecipes version
*/
double gasdev(double m, double v)                            // m=v=u (mean poisson >30)
{
  static int iset=0;
  static float gset=0.0;
  float fac=0.0, r=0.0, v1=0.0, v2=0.0;
  double ran1();
  if(iset==0) 
    {
      do              // Compute probability r=v1^2+v2^2
        {
          v1=2.0*ran1()-1.0;
          v2=2.0*ran1()-1.0;
          r=v1*v1+v2*v2;
        } while(r>=1.0);
      fac=sqrt(-2.0*log(r)/r);
      gset=v1*fac;
      iset=1;
      return(m+sqrt(v)*v2*fac);
    } 
  else
    {
      iset=0;
      return(m+sqrt(v)*gset);
    }
}

void _gaussdev(float *xmv, long n)
{
  /* Returns a normally distributed deviate with zero mean and unit variance,
     using ran1() as the source of uniform deviates. */

  /*  float ran1(long *idum); */
  static int iset=0;
  static float gset;
  float fac,rsq,v1,v2;
  long i;

  for (i=0;i<n;i++) {
    if (iset == 0) {
      do {
  v1=2.0*ran1()-1.0;
  v2=2.0*ran1()-1.0;
  rsq=v1*v1+v2*v2;
      } while (rsq >= 1.0 || rsq == 0.0);
      fac=sqrt(-2.0*log(rsq)/rsq);
      gset=v1*fac;
      iset=1;
      xmv[i] = v2*fac;
    } else {
      iset=0;
      xmv[i] = gset;
    }
  }
}

/*---------*
 *  Pick2  |
 *---------*
 | Pick two numbers (i and j) from n
 | Called in street.c by pick2_chrom
*/
int pick2(int n, int *i, int *j)
{
  double ran1();
  *i=n*ran1();
  while((*j=n*ran1())==*i);
  return(0);
}

float gasdev(long *idum)
/************************************************************
 *	Code of this function is taken from the book			*
 *	"Neumerical Recepies in C", 1992, p289.					*
 *	Returns a normally distributed deviate with zero mean	*
 *	and unit variance, using ran1(idum) as the source of	*
 *	uniform deviates.										*
 ************************************************************/
{
	float ran1(long *idum);
	static int iset=0;
	static float gset;
	float fac, rsq, v1, v2;
	
	if (*idum <0) iset=0;
	if (iset ==0){
		do{
			v1=2.0*ran1(idum)-1.0;
			v2=2.0*ran1(idum)-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq>=1.0 || rsq ==0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	}
	else{
		iset=0;
		return gset;
	}
}

/*
 *----------------------------------------------------------------
 *	Gamma Distribution
 *	Returns a deviate distributed as a gamma distribution
 *	of integer order ia, i.e., a waiting time to the iath
 *	event in a Poisson process of unit mean, using ran1(idum)
 *	as the source of uniform deviates
 *	p. 292 of the book
 *----------------------------------------------------------------
 */
float gamdev(int ia, int *idum)
{
	int j;
	float am,e,s,v1,v2,x,y;
    
	
	if (ia < 1){
		printf("Error in routine GAMDEV!");
		exit ( 0 );
	}
	if (ia < 6) {
		x=1.0;
		for (j=1;j<=ia;j++) x *= ran1(idum);
		x = -log(x);
	} else {
		do {
			do {
				do {
					v1=2.0*ran1(idum)-1.0;
					v2=2.0*ran1(idum)-1.0;
				} while (v1*v1+v2*v2 > 1.0);
				y=v2/v1;
				am=ia-1;
				s=sqrt(2.0*am+1.0);
				x=s*y+am;
			} while (x <= 0.0);
			e=(1.0+y*y)*exp(am*log(x/am)-s*y);
		} while (ran1(idum) > e);
	}
	return x;
}

// normal random variate generator using box-muller
double Random::gran(double s, double m)
{        /* mean m, standard deviation s */
  double x1, x2, w, y1;
  static double y2;
  static int use_last = 0;

  if (use_last)        /* use value from previous call */
    {
      y1 = y2;
      use_last = 0;
    }
  else
    {
      do {
	x1 = 2.0 * ran1() - 1.0;
	x2 = 2.0 * ran1() - 1.0;
	w = x1 * x1 + x2 * x2;
      } while ( w >= 1.0 );
      w = sqrt( (-2.0 * log( w ) ) / w );
      y1 = x1 * w;
      y2 = x2 * w;
      use_last = 1;
    }
  return( m + y1 * s );
}

int16_t awgn(awgn_state_t *s)
{
    double fac;
    double r;
    double v1;
    double v2;
    double amp;

    if (s->iset == 0)
    {
        do
        {
            v1 = 2.0*ran1(s) - 1.0;
            v2 = 2.0*ran1(s) - 1.0;
            r = v1*v1 + v2*v2;
        }
        while (r >= 1.0);
        fac = sqrt(-2.0*log(r)/r);
        s->gset = v1*fac;
        s->iset = 1;
        amp = v2*fac*s->rms;
    }
    else
    {
        s->iset = 0;
        amp = s->gset*s->rms;
    }
    return fsaturate(amp);
}

/**
 * Tries to flip spin i, j. If accepted, the energy and magnetization is updated
 * 
 * @param i
 * @param j
 */
void IsingLattice2D::flipSpin(int i, int j, bool force) {
    double dE =     2*spin_matrix(i,j)*(
                    spin_matrix(i, periodic(j+1)) +
                    spin_matrix(i ,periodic(j-1)) +
                    spin_matrix(periodic(i+1), j) +
                    spin_matrix(periodic(i-1), j));
    if (dE <= 0 || force) {
        M += -2*spin_matrix(i, j);
        E += dE;
        spin_matrix(i,j) *= -1;
        accepted ++;
    }
    else {
        // Check for acceptance;
        if (dE == 4) {
            if (ran1(&idum) <= p4) {
                // accept move
                M += -2*spin_matrix(i, j);
                E += dE;
                spin_matrix(i, j) *= -1;
                accepted ++;
            } else rejected ++;
        }
        if (dE == 8) {
            if (ran1(&idum) <= p8) {
                // Accept move
                M += -2*spin_matrix(i, j);
                E += dE;
                spin_matrix(i, j) *= -1;
                accepted ++;
            } else rejected ++;
        }
    }
}

void IsingLattice2D::Metropolis() {
    for (int n = 0; n<N*N; n++) {
        int i = (int) (ran1(&idum)*N);
        int j = (int) (ran1(&idum)*N);
        flipSpin(i, j, false);
    }
}

/*!
 * Numerical Recipes Gaussian random number generator
 */
double gasdev(long *idum)
{
  static int iset              = 0;
  static double gset;
  double fac, rsq, v1, v2;
  int i;
  double sum;
  if (0) {
  if (*idum<0)            iset = 0;
  if (iset==0) {
    do {
      v1                       = 2.0*ran1(idum)-1.0;
      v2                       = 2.0*ran1(idum)-1.0;
      rsq                      = v1*v1 + v2*v2;
    } while ((rsq>=1.0) || (rsq==0.0));
    fac                        = sqrt(-2.0*log(rsq)/rsq);
    gset                       = v1*fac;
    iset                       = 1;
    return v2*fac;
  } else {
    iset                       = 0;
    return gset;
  }
  } else {
    sum = 0;
    for (i=0; i<12; ++i) {
      sum += ran1(idum)-0.5;
    }
    return sum;
  }


}

/* Next function is from numerical recipes in C. */
float randn_internal(long *idum, int reset) {
  static int iset=0;
  static float gset;
  float fac, rsq, v1, v2;

  if (reset) {
    iset = 0;
  }

  if (iset==0) {
    do {
      v1 = 2.0*ran1(idum, reset)-1.0;
      v2 = 2.0*ran1(idum, reset)-1.0;
      rsq = v1*v1+v2*v2;
    } while(rsq >= 1.0 || rsq == 0.0);

    fac = sqrt(-2.0*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return v2*fac;
  } else {
    iset = 0;
    return gset;
  }
}

int main(void)
{
	float d,**a,**al,*b,*x;
	unsigned long i,j,*indx;
	long idum=(-1);

	a=matrix(1,7,1,4);
	x=vector(1,7);
	b=vector(1,7);
	al=matrix(1,7,1,2);
	indx=lvector(1,7);
	for (i=1;i<=7;i++) {
		x[i]=ran1(&idum);
		for (j=1;j<=4;j++) {
			a[i][j]=ran1(&idum);
		}
	}
	banmul(a,7,2,1,x,b);
	for (i=1;i<=7;i++) printf("%ld %12.6f %12.6f\n",i,b[i],x[i]);
	bandec(a,7,2,1,al,indx,&d);
	banbks(a,7,2,1,al,indx,b);
	for (i=1;i<=7;i++) printf("%ld %12.6f %12.6f\n",i,b[i],x[i]);
	free_lvector(indx,1,7);
	free_matrix(al,1,7,1,2);
	free_vector(b,1,7);
	free_vector(x,1,7);
	free_matrix(a,1,7,1,4);
	return 0;
}

float gamdev(int ia, long *idum)
{
	float ran1(long *idum);
	void nrerror(char error_text[]);
	int j;
	float am,e,s,v1,v2,x,y;

	if (ia < 1) nrerror("Error in routine gamdev");
	if (ia < 6) {
		x=1.0;
		for (j=1;j<=ia;j++) x *= ran1(idum);
		x = -log(x);
	} else {
		do {
			do {
				do {
					v1=ran1(idum);
					v2=2.0*ran1(idum)-1.0;
				} while (v1*v1+v2*v2 > 1.0);
				y=v2/v1;
				am=ia-1;
				s=sqrt(2.0*am+1.0);
				x=s*y+am;
			} while (x <= 0.0);
			e=(1.0+y*y)*exp(am*log(x/am)-s*y);
		} while (ran1(idum) > e);
	}
	return x;
}