function agauss(x,av,sig)
c evaluate gaussian probability integral
c using routine from Bevington : Data reduction and Error 
c Analysis for the Physical Sciences (p48) 
c	x- limit for integral
c	av- distribution mean
c	sig- standard deviation of distribution
c returns (1/sqrt(2.*pi)) integral from -z to z of exp(-.5*x**2)dx
c where z=abs(x-av)/sig
c normalized so z=infinity gives 1
	z=abs(x-av)/sig
	agauss=0.
	if (z) 42,42,21
21	term=.7071067812*z
22	sum=term
	y2=(z**2.)/2.
	denom=1.
c accumulate sum of terms
31	denom=denom+2.
	term=term*(y2*2./denom)
	sum=sum+term
	if (term/sum -1.e-5) 41,41,31
41	agauss= 1.128379167*sum*exp(-1.*y2)
42	return
	end