/* -------------------------------------------------------------------------
! Program to calculate the first zeroes (root abscissas) of the first 
! kind Bessel function of integer order N using the subroutine ROOTJ.
! --------------------------------------------------------------------------
! SAMPLE RUN:
!
! (Calculate the first 10 zeroes of 1st kind Bessel function of order 2).
! 
! Zeroes of Bessel Function of order: 2
! 
! Number of calculated zeroes: 10 
!
! Table of root abcissas (5 items per line)
!    5.135622    8.417244   11.619841   14.795952   17.959819
    21.116997   24.270112   27.420574   30.569204   33.716520
! 
! Table of error codes (5 items per line)
!   0   0   0   0   0
!   0   0   0   0   0
!
! --------------------------------------------------------------------------
! Reference: From Numath Library By Tuan Dang Trong in Fortran 77
!            [BIBLI 18].
!
!                               C++ Release 1.0 By J-P Moreau, Paris
!                                          (www.jpmoreau.fr)
! ------------------------------------------------------------------------ */
#include <stdio.h>
#include <math.h>

#define SIZE  11
    
    double JZ[SIZE];
    int IE[SIZE];
    int i,j, N,NR;

    void SECANT(int N,int NITMX, double TOL, double *ZEROJ, int *IER);
    double BESSJ(int N, double X);
    double BESSJ0(double X);
    double BESSJ1(double X);

    double Sign(double a, double b) {
      if (b>=0.0) return (fabs(a));
      else return (-fabs(a));
    }

//----------------------------------------------------------------------
    void ROOTJ(int N, int NK, double *JZERO, int *IER)  {
/*----------------------------------------------------------------------
!     CALCULATE THE FIRST NK ZEROES OF BESSEL FUNCTION J(N,X)
!
!     INPUTS:
!       N    ORDER OF FUNCTION J (INTEGER >= 0)                  I*4
!       NK   NUMBER OF FIRST ZEROES  (INTEGER > 0)               I*4
!     OUTPUTS:
!       JZERO(NK)  TABLE OF FIRST ZEROES (ABCISSAS)              R*8
!       IER(NK)    TABLE OF ERROR CODES (MUST BE ZEROES)         I*4
!
!     REFERENCE :
!     ABRAMOWITZ M. & STEGUN IRENE A.
!     HANDBOOK OF MATHEMATICAL FUNCTIONS
! -------------------------------------------------------------------- */
    double ZEROJ,B0,B1,B2,B3,B5,B7,T0,T1,T3,T5,T7,FN,FK;
    double C1,C2,C3,C4,C5,F1,F2,F3,TOL,ERRJ, PI;
    int IERROR,K,NITMX;

      TOL=1E-8; NITMX=10; PI=4.0*atan(1.0);
      C1=1.8557571; C2=1.033150; C3=0.00397; C4=0.0908; C5=0.043;
      FN = 1.0*N;

//    FIRST ZERO
      if (N==0) {
        ZEROJ = C1+C2-C3-C4+C5;
        SECANT(N,NITMX,TOL,&ZEROJ,&IERROR);
 	    IER[1]=IERROR;
        JZERO[1]=ZEROJ;
      }
      else {
        F1 = pow(FN,(1.0/3.0));
        F2 = F1*F1*FN;
        F3 = F1*FN*FN;
        ZEROJ = FN+C1*F1+(C2/F1)-(C3/FN)-(C4/F2)+(C5/F3);
        SECANT(N,NITMX,TOL,&ZEROJ,&IERROR);
        IER[1]=IERROR;
        JZERO[1]=ZEROJ;
      }

      T0 = 4.0*FN*FN;
      T1 = T0-1.0;
      T3 = 4.0*T1*(7.0*T0-31.0);
      T5 = 32.0*T1*((83.0*T0-982.0)*T0+3779.0);
      T7 = 64.0*T1*(((6949.0*T0-153855.0)*T0+1585743.0)*T0-6277237.0);

//    OTHER ZEROES
      for (K = 2; K <= NK; K++) {
        JZERO[K] = 0.0;
        FK = 1.0*K;

//      MAC MAHON'S SERIES FOR K>>N

        B0 = (FK+0.5*FN-0.25)*PI;
        B1 = 8.0*B0;
        B2 = B1*B1;
        B3 = 3.0*B1*B2;
        B5 = 5.0*B3*B2;
        B7 = 7.0*B5*B2;

        ZEROJ = B0-(T1/B1)-(T3/B3)-(T5/B5)-(T7/B7);

        ERRJ=fabs(BESSJ(N,ZEROJ));

//      IMPROVE SOLUTION USING PROCEDURE SECANT

        if (ERRJ > TOL) SECANT(N,NITMX,TOL,&ZEROJ,&IERROR);
        JZERO[K]=ZEROJ;
        IER[K]=IERROR;
      }
    }
// ------------------------------------------------------------------------------
  void SECANT(int N,int NITMX, double TOL, double *ZEROJ, int *IER)  {
    //Labels: e5,e10,e15,e20
    
    double P0,P1,Q0,Q1,DP,P;
    double C[3];
    int IT,NEV,NTRY;

      C[1]=0.95; C[2]=0.999;
      NTRY=1;
      *IER=0;

e5:   P0 = C[NTRY]*(*ZEROJ);

      P1 = *ZEROJ;
      NEV = 2;
      Q0 = BESSJ(N,P0);
      Q1 = BESSJ(N,P1);
      for (IT = 1; IT <= NITMX; IT++) {
        if (Q1==Q0) goto e15;
        P = P1-Q1*(P1-P0)/(Q1-Q0);
        DP = P-P1;
        if (IT==1) goto e10;
        if (fabs(DP) < TOL) goto e20;

e10:    NEV = NEV+1;
        P0 = P1;
        Q0 = Q1;
        P1 = P;
        Q1 = BESSJ(N,P1);
      }
e15:  NTRY++;
      if (NTRY <= 2) goto e5;
      *IER = NTRY;
e20:  *ZEROJ = P;
  }

// --------------------------------------------------------------------
    double BESSJ(int N, double X) {

/*    THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION 
      OF ORDER N, INTEGER FOR ANY REAL X. WE USE HERE THE CLASSICAL 
      RECURRENT FORMULA, WHEN  X > N. FOR X < N, THE MILLER'S ALGORITHM
      IS USED TO AVOID OVERFLOWS.
      REFERENCE :
      C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
      MATHEMATICAL TABLES, VOL.5, 1962.  */
      int IACC = 40; 
  	  double BIGNO = 1e10, BIGNI = 1E-10;
      double TOX,BJM,BJ,BJP,SUM,TMP;
      int J,JSUM,M;

      if (N==0) return BESSJ0(X);
      if (N==1) return BESSJ1(X);
      if (X==0.0) return 0.0;

      TOX = 2.0/X;
      if (X > N) {
        BJM = BESSJ0(X);
        BJ  = BESSJ1(X);
        for (J = 1; J<N; J++) {
          BJP = J*TOX*BJ-BJM;
          BJM = BJ;
          BJ  = BJP;
          TMP = BJ;
        }
		return TMP;
      }
      else {
        M = (int) (2*((N+floor(sqrt(IACC*N))) / 2));
        TMP = 0.0;
        JSUM = 0;
        SUM = 0.0;
        BJP = 0.0;
        BJ  = 1.0;
        for (J = M; J>0; J--) {
          BJM = J*TOX*BJ-BJP;
          BJP = BJ;
          BJ  = BJM;
          if (fabs(BJ) > BIGNO) {
            BJ  *= BIGNI;
            BJP *= BIGNI;
            TMP *= BIGNI;
            SUM *= BIGNI;
          }
          if (JSUM != 0)  SUM += BJ;
          JSUM = 1-JSUM;
          if (J == N) TMP = BJP;
          SUM = 2.0*SUM-BJ;
        }
        return (TMP/SUM);
      }
    }
// ----------------------------------------------------------------------
    double BESSJ0(double X)  {
      
/*    THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION 
      OF ORDER 0 FOR ANY REAL X. WE USE HERE THE POLYNOMIAL APPROXIMATION
      BY SERIES OF CHEBYSHEV POLYNOMIALS FOR 0<X<8 AND 0<8/X<1.
      REFERENCES :
      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
      VOL.5, 1962.  */

      double AX,FR,FS,Z,FP,FQ,TMP,XX;
      double Y,P1,P2,P3,P4,P5,R1,R2,R3,R4,R5,R6,Q1,Q2,Q3,Q4,Q5,S1,S2,S3,S4,S5,S6;

      P1=1.0; P2=-0.1098628627E-2; P3=0.2734510407E-4;
      P4=-0.2073370639E-5; P5=0.2093887211E-6;
      Q1=-0.1562499995E-1; Q2=0.1430488765E-3; Q3=-0.6911147651E-5;
      Q4=0.7621095161E-6; Q5=-0.9349451520E-7;
      R1=57568490574.0; R2=-13362590354.0; R3=651619640.7;
      R4=-11214424.18; R5=77392.33017; R6=-184.9052456;
      S1=57568490411.0; S2=1029532985.0; S3=9494680.718;
      S4=59272.64853; S5=267.8532712; S6=1.0;

      if (X==0.0) return 1.0;
      AX = fabs(X);
      if (AX < 8) {
        Y = X*X;
        FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
        FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
        TMP = FR/FS;
      }
      else {
        Z = 8./AX;
        Y = Z*Z;
        XX = AX-0.785398164;
        FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
        FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
        TMP = sqrt(0.636619772/AX)*(FP*cos(XX)-Z*FQ*sin(XX));
      }
      return TMP;
    }
// ----------------------------------------------------------------------
    double BESSJ1(double X)  {
      
/*    THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION 
      OF ORDER 0 FOR ANY REAL X. WE USE HERE THE POLYNOMIAL APPROXIMATION
      BY SERIES OF CHEBYSHEV POLYNOMIALS FOR 0<X<8 AND 0<8/X<1.
      REFERENCES :
      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
      VOL.5, 1962.  */
    
      double AX,FR,FS,Z,FP,FQ,TMP,XX;
      double Y,P1,P2,P3,P4,P5,P6,R1,R2,R3,R4,R5,R6,Q1,Q2,Q3,Q4,Q5,S1,S2,S3,S4,S5,S6;

      P1=1.0; P2=0.183105E-2; P3=-0.3516396496E-4;
      P4=0.2457520174E-5; P5=-0.240337019E-6; P6=0.636619772;
      Q1=0.04687499995; Q2=-0.2002690873E-3; Q3=0.8449199096E-5;
      Q4=-0.88228987E-6; Q5=0.105787412E-6;
      R1=72362614232.0; R2=-7895059235.0; R3=242396853.1;
      R4=-2972611.439; R5=15704.48260; R6=-30.16036606;
      S1=144725228442.0; S2=2300535178.0; S3=18583304.74;
      S4=99447.43394; S5=376.9991397; S6=1.0;

      AX = fabs(X);
      if (AX < 8) {
        Y = X*X;
        FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
        FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
        TMP = X*(FR/FS);
      }
      else {
        Z = 8./AX;
        Y = Z*Z;
        XX = AX-2.35619491;
        FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
        FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
        TMP = sqrt(P6/AX)*(cos(XX)*FP-Z*sin(XX)*FQ)*Sign(S6,X);
      }
	  return TMP;
    }

void main()  {

  N=2;
  NR=10;

  ROOTJ(N,NR,JZ,IE);

  printf("\n Zeroes of Bessel Function of order: %d\n\n", N);
  printf(" Number of calculated zeroes: %d\n\n", NR);
  printf(" Table of root abcissas (5 items per line)\n");
  j=1;
  for (i=1; i<=NR; i++) {
    printf("%12.6f", JZ[i]); j++;
    if (j % 6 == 0) {
      printf("\n"); j=1;
    }
  }
  printf("\n Table of error codes (5 items per line)\n");
  j=1;
  for (i=1; i<=NR; i++) {
    printf("   %d",IE[i]); j++;
    if (j % 6 == 0) {
      printf("\n"); j=1;
    }
  };
  printf("\n");

}

// end of file trootj.cpp