/* @(#)erf.c	17.1.1.1 (ES0-DMD) 01/25/02 17:49:05 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* % (ESO-DAG) 01/25/02 17:49:05
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENTIFICATION: subroutine erfunc.c
.KEYWORDS:       error function
.PURPOSE:        return the error function erfunc(x);
.USE:            erfunc(x)
.NOTE:           The code is taken for Numerical Recipes in C 
                 The routines concerned are: erf, gammp, gammq, gcf, gser
------------------------------------------------------------------------------ */

#include <math.h>

#define ITMAX 100
#define EPS   3.0e-7

float erfunc(x)
float x;
{
	float gammp();

	return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}


float gammp(a,x)
float a,x;
{
	float gamser,gammcf,gln;
	void gser(),gcf();

	if (x < 0.0 || a <= 0.0) 
           {
           SCTPUT("*** ERROR: Invalid arguments in routine GAMMP");
           return;
    	   }
	if (x < (a+1.0)) {
		gser(&gamser,a,x,&gln);
		return gamser;
	} else {
		gcf(&gammcf,a,x,&gln);
		return 1.0-gammcf;
	}
}

float gammq(a,x)
float a,x;
{
	float gamser,gammcf,gln;
	void gcf(),gser();

	if (x < 0.0 || a <= 0.0)
           {
           SCTPUT("*** ERROR: Invalid arguments in routine GAMMQ");
           return;
    	   }
	if (x < (a+1.0)) {
		gser(&gamser,a,x,&gln);
		return 1.0-gamser;
	} else {
		gcf(&gammcf,a,x,&gln);
		return gammcf;
	}
}



void gser(gamser,a,x,gln)
float a,x,*gamser,*gln;
{
	int n;
	float sum,del,ap;
	float gammln();

	*gln = gammln(a);
	if (x <= 0.0) {
		if (x < 0.0) 
                   { 
                   SCTPUT("*** ERROR: Invalid arguments in routine GSER");
                   return;
    	           }
		*gamser = 0.0;
		return;
	} else {
		ap = a;
		del = sum = 1.0/a;
		for (n=1;n<=ITMAX;n++) {
			ap  += 1.0;
			del *= x/ap;
			sum += del;
			if (fabs(del) < fabs(sum)*EPS) {
				*gamser=sum*exp(-x+a*log(x)-(*gln));
				return;
			}
		}
                SCTPUT(
                   "*** ERROR: a too large, ITMAX too small in routine GSER");
		return;
	}
}


void gcf(gammcf,a,x,gln)
float a,x,*gammcf,*gln;
{
	int n;
	float gammln();
	float gold=0.0,g,fac=1.0,b1=1.0;
	float b0=0.0,anf,ana,an,a1,a0=1.0;

	*gln = gammln(a);
	a1   = x;
	for (n=1;n<=ITMAX;n++) {
		an  = (float) n;
		ana = an-a;
		a0  = (a1+a0*ana)*fac;
		b0  = (b1+b0*ana)*fac;
		anf = an*fac;
		a1  = x*a0+anf*a1;
		b1  = x*b0+anf*b1;
		if (a1) {
			fac = 1.0/a1;
			g   = b1*fac;
			if (fabs((g-gold)/g) < EPS) {
				*gammcf=exp(-x+a*log(x)-(*gln))*g;
				return;
			}
			gold = g;
		}
	}
   SCTPUT("*** ERROR: a too large, ITMAX too small in routine GCF");
   return;
}


float gammln(xx)
float xx;
{
	double x,tmp,ser;
	static double cof[6]={76.18009173,-86.50532033,24.01409822,
		-1.231739516,0.120858003e-2,-0.536382e-5};
	int j;

	x    = xx-1.0;
	tmp  = x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser  = 1.0;
	for (j=0;j<=5;j++) {
		x += 1.0;
		ser += cof[j]/x;
	}
	return (float) (-tmp + log(2.50662827465 * ser));
}


#undef ITMAX
#undef EPS