#include <math.h>

float rangauss(long *idum, float mean, float sigma)

{
  float ran1(long *idum);
  float sigrange, a, x, b, p, rgauss;

  b=1.;
  p=0.;
  sigrange=5.;

  do {

    a=ran1(idum);
    x=2.*(a-0.5)*sigma*sigrange + mean;
    b=ran1(idum);
    p=exp(-(x-mean)*(x-mean)/(2.*sigma*sigma));

  } while (b>p);

  rgauss=x;
  return rgauss;

}