#include <math.h>

static double *yloc, f0loc, scale, lprec1;

static double fun(double *y1, int *flag)
    {
	double y1sc, z, s, v, lerch;
	int iter;
	extern double *yloc, f0loc, scale, lprec1;
		
	y1sc = (*y1) * scale;
	z = pow(f0loc,yloc[2]);
	s = 1;
	v = y1sc/yloc[2];
	*flag = lerchphi(&z, &s, &v, &lprec1, &lerch, &iter);
	if (*flag) return 1.0;
	else return (yloc[4] - pow(f0loc,y1sc)/yloc[2]/yloc[3]*lerch);
    }
	
int zqnt(double *y, double *f0, double *lprec, double *ztol, double *sol)
    {
      	double yu, yl, ya, yb, scale1, factor, funu, funl, funup, funlp, yascaled, ybscaled;
	extern double *yloc, f0loc, scale, lprec1;
	int flag = 0;
	double zeroin_();
		
	lprec1 = *lprec;
	f0loc = *f0;
	yloc = y;
								
	yu = 1.0;
	yl = 1.0;
	factor = 10.0;
	scale = 1.0;
	scale1 = 0.1;
	flag = 0;
	funl = fun(&yl,&flag);
	if (flag)
	    {
		*sol = 1.0;
		return flag;
	    }
	flag = 0;
	funu = fun(&yu,&flag);
	if (flag)
	    {
		*sol = 1.0;
		return flag;
	    }	
	funlp = funl;
	funup = funu;
										
	while (funu*funup > 0.0 && funl*funlp > 0.0)
	    {
		funlp = funl;
		funup = funu;	
		scale1 *= factor;
		yu *= factor;
		yl /= factor;
		flag = 0;
		funl = fun(&yl,&flag);
		if (flag)
		    {
			*sol = 1.0;
			return flag;
		    }
		flag = 0;
		funu = fun(&yu,&flag);
		if (flag)
		    {
			*sol = 1.0;
			return flag;
		    }
	    }
				
	if (funu*funup <= 0.0) 
	    {
		ya = yu/factor;
		yb = yu;
		scale = scale1;
	    }
	if (funl*funlp <= 0.0) 
	    {
		ya = yl;
		yb = yl*factor;
		scale = 1/factor/scale1;
	    }
										
	yascaled = ya/scale;
	ybscaled = yb/scale;
	*sol = zeroin_(&yascaled, &ybscaled, &fun, ztol, flag)*scale;
	if (flag) 
	    {
		*sol = 1.0;
		return flag;
	    }
	
	return flag;
		
    }