/* the functions int this file allow to calculate 
   the hypergeometric function */

#include <hypgeo.h>
#define EPS 1.0e-6
#define ONE Complex(1.0,0.0)

dcomplex aa,bb,cc,z0,dz;

int kmax,kount;
double *xp,**yp,dxsav;

dcomplex hypgeo(dcomplex a, dcomplex b, dcomplex c, dcomplex z)
{
	void bsstep(double y[], double dydx[], int nv, double *xx, double htry,
		double eps, double yscal[], double *hdid, double *hnext,
		void (*derivs)(double, double [], double []));
	void hypdrv(double s, double yy[], double dyyds[]);
	void hypser(dcomplex a, dcomplex b, dcomplex c, dcomplex z,
		dcomplex *series, dcomplex *deriv);
	void odeint(double ystart[], int nvar, double x1, double x2,
		double eps, double h1, double hmin, int *nok, int *nbad,
		void (*derivs)(double, double [], double []),
		void (*rkqs)(double [],double [],int,double *,double, double,
		double [],double *,double *,void (*)(double,double [],double [])));
	int nbad,nok;
	dcomplex ans,y[3];
	double *yy;

	kmax=0;
	if (z.re*z.re+z.im*z.im <= 0.25) {
		hypser(a,b,c,z,&ans,&y[2]);
		return ans;
	}
	else if (z.re < 0.0) z0=Complex(-0.5,0.0);
	else if (z.re <= 1.0) z0=Complex(0.5,0.0);
	else z0=Complex(0.0,z.im >= 0.0 ? 0.5 : -0.5);
	aa=a;
	bb=b;
	cc=c;
	dz=Csub(z,z0);
	hypser(aa,bb,cc,z0,&y[1],&y[2]);
	yy=dvector(1,4);
	yy[1]=y[1].re;
	yy[2]=y[1].im;
	yy[3]=y[2].re;
	yy[4]=y[2].im;
	odeint(yy,4,0.0,1.0,EPS,0.1,0.0001,&nok,&nbad,hypdrv,bsstep);
	y[1]=Complex(yy[1],yy[2]);
	free_dvector(yy,1,4);
	return y[1];
}


void hypser(dcomplex a, dcomplex b, dcomplex c, dcomplex z, dcomplex *series,
	    dcomplex *deriv)
{
	void nrerror(char error_text[]);
	int n;
	dcomplex aa,bb,cc,fac,temp;

	deriv->re=0.0;
	deriv->im=0.0;
	fac=Complex(1.0,0.0);
	temp=fac;
	aa=a;
	bb=b;
	cc=c;
	for (n=1;n<=1000;n++) {
		fac=Cmul(fac,Cmul(aa,Cmul(bb,Cinv(cc))));
		deriv->re +=fac.re;
		deriv->im +=fac.im;
		fac=Cmul(fac,CRmul(1.0/n,z));
		*series=Cadd(temp,fac);
		if (series->re == temp.re && series->im == temp.im) return;
		temp= *series;
		aa=Cadd(aa,ONE);
		bb=Cadd(bb,ONE);
		cc=Cadd(cc,ONE);

	}
	nrerror("convergence failure in hypser");
}


#undef ONE
#undef EPS

/* (C) Copr. 1986-92 Numerical Recipes Software V7&584. */