/* @(#)besz.c	1.12 (8/9/93) */
/* BESZ 
   fast zero finder for bessel functions J_nu(x)
   command line switches:
   -n number of zeroes 
   -d dimension of working array (d >> 4 n)
   -i index nu (double)
   -v verbose (gives deviation from asymptotic estimate)
   -p precision (1.e-8 ... 1.e-18)
   -F enforce d > 5n+offset
   -s dump setup

   usage: besz [-p <prec>][-v][-i <index>][-d <dim>][-n <num>]
   algorithm: inverse recursion + asymptotic estimate for speeding up

   E.Onofri (c) 1993 (onofri@parma.infn.it)

   */


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

#define M_PI  3.1415926535897931160E0  
#define sqr(x) ((x)*(x))
#define cube(x) ((x)*(x)*(x))
#define PREC 1.E-16
#define NOFFSET 50
#define CLEAN 0
#define VERBOSE 1
#define FALSE 0
#define TRUE  1


char    mode = CLEAN;
char    force = 0;
int	npts =  500;
int	neig =  20;
int     Nd ;
int     offs = NOFFSET;
char    dump = FALSE;
double  tol,maxerror = 0.0;

double bes_nu=0.0, prec;

void examine_command_line(argc,argv)
int argc;
char **argv;
{
  extern char *optarg;
  extern int  optind;

  int c;
  while ((c = getopt(argc, argv, "svFi:n:d:p:o:")) != EOF) {

    switch(c) {

    case 'v':
      mode = VERBOSE;
      break;

    case 'F':
      force=TRUE;
      break;

    case 'i':
      bes_nu = atof(optarg);
      break;
      
    case 'd':
      npts = atoi(optarg);
      break;

    case 'n':
      neig = atoi(optarg);
      break;
      
    case 'o':
      offs=atoi(optarg);
      break;

    case 's':
      dump = TRUE;
      break;


    case 'p':
      prec  = atof(optarg)/100;
      if(prec<1.0E-18) {
	fprintf(stderr,"Accuracy cannot be that high!\n");
	fprintf(stderr,"assuming prec=%le\n",PREC);
	prec=PREC/100;
      }
      break;
      
    default:
      fprintf(stderr,"Usage: besz [-options]\nOptions:\n");
      fprintf(stderr,"       -i<double>  bessel index (float)    \n");
      fprintf(stderr,"       -d<int>     dimension  (int)      \n");
      fprintf(stderr,"       -n<int>     number of zeroes (int) \n");
      fprintf(stderr,"       -p<double>  precision (%le)   \n");
      fprintf(stderr,"       -F          enforce dim>5 n + offs   \n");
      fprintf(stderr,"       -v          verbose output   \n");
      fprintf(stderr,"       -s          dump setup (stderr)\n");

      exit(1);
      break;
    }
  }
}



double	bisec(k,Max)
int k;
double Max;

{

  register int i, index;
  register double p0, p1, p2; 
  register double x, xmin=0.0, xmax=2.0*Max;
  register double delta=xmax-xmin;

  while( delta > prec) {

    x = xmin+xmax;

    index=1;
    
    p0=1.0;
    p1=x*bes_nu;   
    
    for(i = 1;i< Nd;++i) {  

      p2 = x*p1*((double)i+bes_nu) - p0;

      if(p1*p2<0)
	index++;
      
      p0 = p1;
      p1 = p2;
    }

    if(index <= k)
      xmax = 0.5*x;
    else
      xmin = 0.5*x;
    delta*=0.5;
  }
   
  return x*0.5;
    
}
  
void main(argc,argv)
int argc;
char **argv;

{
  char format[60], str[30];
  double eigenv, error;
  int i=1,j,k;
  double tmp, *asy;
  double Max = 10.0;
  prec=PREC;

  examine_command_line(argc,argv);


  if(force)
    npts=5*neig+NOFFSET;
  else
    neig=(neig*5<npts+NOFFSET ? neig : (npts-NOFFSET)/5);

  Nd  = NOFFSET;
  bes_nu+=1.0;
  tol=prec*10;
  
  asy=(double *)calloc(neig+1,sizeof(double));


  for(k=1;k<=neig;++k) {
    double mu,nu;
    double beta;
    nu=bes_nu-1.0;
    mu=4.0*sqr(nu);
    beta=(k+0.5*nu-0.25)*M_PI;
    asy[k]= beta - (mu-1.0)/8./beta - 4.0*(mu-1.0)*(7.0*mu-31.0)/3.0/cube(8.0*beta);
    asy[k] -=  32.*(mu-1.0)*(83.0*sqr(mu)-982*mu+3779.0)/15./cube(8.0*beta)/sqr(8.0*beta);
    asy[k] -= 64.*(mu-1.0)*(6949.*cube(mu)-153855.0*sqr(mu)+1585743*mu-6277237)/105./sqr(cube(8.0*beta))/(8.*beta); 

  }

  if((int)log10(1./tol/asy[k-1])<2) {
    fprintf(stderr,"Warning: you are asking for too many roots \n");
    fprintf(stderr,"         output will be truncated          \n");
    fprintf(stderr,"         Hit any key to continue \n");
    getchar();
  }
    
  
  for(i=1;i<=neig;++i) {
    Max=eigenv= bisec(i,Max);
    Nd+=5;
    error = fabs(eigenv*asy[i]-1.0);

    sprintf(str,"%d.%d",(int)log10(1./tol), (int)log10(1./tol*eigenv)-1); 

    if(mode==VERBOSE) {
      strcpy(format,"z_{%g,%d}:\t %");
      strcat(format,str);
      strcat(format,"f : %g\n\0");
      printf(format,bes_nu-1.0, i,1./eigenv,error);
    }
    else {
      strcpy(format,"%d\t %");
      strcat(format,str);
      strcat(format,"f \n\0");
      printf(format,i,1./eigenv);

    }
    if(fabs(error)<tol)
      break;
  }

  for(++i;i<=neig;++i) {
   sprintf(str,"%d.%d",(int)log10(1./tol), (int)log10(1./tol/asy[i])-1);
    if(mode==VERBOSE) {
      strcpy(format,"z_{%g,%d}:\t %");
      strcat(format,str);
      strcat(format,"f : %g\n\0");
      printf(format,bes_nu-1.0, i,asy[i],0.0);
    }
    else {
      strcpy(format,"%d\t %");
      strcat(format,str);
      strcat(format,"f \n\0");
      printf(format,i,asy[i]);
    }

  }
  if(dump==TRUE) {
    fprintf(stderr,"Present setup:\n");
    fprintf(stderr,"besz %s",(mode==CLEAN?"":"-v "));
    fprintf(stderr,"%s",(force!=TRUE?"":"-F "));
    fprintf(stderr,"%s",(dump!=TRUE?"":"-s "));
    fprintf(stderr,"-n %d ",neig);
    fprintf(stderr,"-d %d ",npts);
    fprintf(stderr,"-i %g ",bes_nu-1.0);
    fprintf(stderr,"-o %d ",offs);
    fprintf(stderr,"-p %g \n\n",tol);

  }

}