#include <stdio.h>
#include <math.h>


/* 	OUTLINE:

 We want to find any zeros in some interval [Xo,X1]. 
 We call the following slatec function from C. We supply in with an initial 
 guess, and criteria for the error, and integer variable that allows the
 function to respond to its exiting status (see below for rough meanings).
 
 Slatec gets back to us, hopefully, by putting the answer into Xo. In this
 case, the code in FLAG is reported as 1, or 2. (For case 1, the interval
 [Xo,X1] where the zero is pinned down to, has been shrunk to specified
 tolerance by dfzero.) For Case 2, the function at Xo is exactly zero now;
 dfzero did not care to change X1, however, so that the interval [Xo,X1] 
 may not not neccessarily have shrunk to specification. 
 Slatec may also get back to us by returning a diverging point, signaled by 
 FLAG=3. FLAG=4 signals that the interval has collapsed but that there are 
 possibly no zeros within the interval; no change of sign has been detected 
 in this interval; this still leaves the possibility of zeroes of even order,
 which is up to the user to test.
 FLAG=5 may be returned when change when too many iterations have been
 performed, without success. The user may want to split up the interval into
 smaller pieces.
 
           FZERO(double (&funct)(double*),
 		 double &Xo,
 		 double &X1,
 		 double &Xguess,
 		 double &RE --relative error, of how close funct(Xo) is to zero,
 		 double &AE --absolute error, how big the diff btw Xo, X1, 
 		      	      when done with the procedure,
 		 int &FLAG
 	      )

FLAG has the following meaning:
1 : 	good deal; the value of the function is close enough to zero for x
	in the interval [Xo,X1], and |X1-Xo| has become as small or smaller 
	as we have requested. 

2 :     the value of f at Xo is exactly zero. (X1 ignored, not to be used.)
        
3 :     The function changes sign in [X0,X1], but increases in magnitude as we 
	shrink the interval. The function probably blows up here.


4 :	The interval [X0,X1 has collapsed to a small number; however no change
	in sign was found to occur within it. This may or may not be a zero of
	even multiplicity.

5 :	No zeros found within allowed program cycles. Program quit.

*/


/* DEFINE a DATATYPE FOR (pointer to) double double FORTRAN type FUNCTION */
typedef double ddfortran(double*);
typedef double (*pddfortran)(double*);

void dfzero(ddfortran*,double*,double*,double*, double*,double*,int*);

ddfortran ddfortranSin;   /* example of fortran double double type function: */
double ddfortranSin(double* x) {   /*   sin() in fortran format   	     */
    double y;
    y=(*x);
    return(sin(y));
}



int dfZero(ddfortran* pfun, double* a, double* b, double re, double ae) { 
  int flag;
  double noguess, A,B;
  noguess=*a; A=*a; B=*b;
  dfzero(pfun, a, b, &noguess, &re, &ae, &flag);
  if (flag==5) {
     puts("Warning: dfzero has crashed while evaluating the interval ");
     printf("(%g,%g).\n Any zeros in that interval will be ignored.\n",A,B);
  }
  return(flag);
} 
 

 
int TestProg0() {    /* A TEST PROGRAM FOR DFZERO                       */
float f; double a,b,guess,re,ae,A,B; int flag;
ddfortran* pFun;
pFun=&ddfortranSin;
 re=.001; ae=.001;
 puts("What is your initial guess? range min? max?\n");
 scanf("%f",&f); guess=1.*f;
 scanf("%f",&f); A=a=1.*f;
 scanf("%f",&f); B=b=1.*f; 
 printf("calling dfzero(&Sin,[%g],[%g],[%g],[%g],&flag)\n",
 					  a,  b, guess, re, ae );
 dfzero(pFun, &a, &b, &guess, &re, &ae, &flag); /* CALL DFZERO: */
 printf("===>  [Xo,X1]==>[%g,%g]\n flag=%d\n\n",a,b,flag);
 
 printf("calling dfZero(&Sin,[%g],[%g],%g,%g)\n", a,b,re,ae);
 dfZero(pFun, &A, &B, re, ae); 
 printf("===>  [Xo,X1]==>[%g,%g]\n\n",A,B);
 return 0;
}
 
int main() { TestProg0(); return 0;}

/* compile me:  cc file.c -lm -lxlf90 -lslatec */