# include # include # include # include "besselzero.h" /******************************************************************************/ void jn_zeros ( int n, int nt, double z[] ) /******************************************************************************/ /* Purpose: jn_zeros() computes the first NT zeros of Bessel function Jn(x). Licensing: This code is distributed under the MIT license. Modified: 08 June 2025 Author: John Burkardt Input: integer N, the order of the Bessel functions. integer NT, the number of zeros. Output: real Z(NT), the first NT zeros of Jn(x). */ { double *rj1; double *ry0; double *ry1; rj1 = ( double * ) malloc ( nt * sizeof ( double ) ); ry0 = ( double * ) malloc ( nt * sizeof ( double ) ); ry1 = ( double * ) malloc ( nt * sizeof ( double ) ); jyzo ( n, nt, z, rj1, ry0, ry1 ); free ( rj1 ); free ( ry0 ); free ( ry1 ); return; } /******************************************************************************/ void yn_zeros ( int n, int nt, double z[] ) /******************************************************************************/ /* Purpose: yn_zeros() computes the first NT zeros of Bessel function Yn(x). Licensing: This code is distributed under the MIT license. Modified: 08 June 2025 Author: John Burkardt Input: integer N, the order of the Bessel functions. integer NT, the number of zeros. Output: real Z(NT), the first NT zeros of Yn(x). */ { double *rj0; double *rj1; double *ry1; rj0 = ( double * ) malloc ( nt * sizeof ( double ) ); rj1 = ( double * ) malloc ( nt * sizeof ( double ) ); ry1 = ( double * ) malloc ( nt * sizeof ( double ) ); jyzo ( n, nt, rj0, rj1, z, ry1 ); free ( rj0 ); free ( rj1 ); free ( ry1 ); return; } /******************************************************************************/ void jyndd ( int n, double x, double *bjn, double *djn, double *fjn, double *byn, double *dyn, double *fyn ) /******************************************************************************/ /* jyndd(): Bessel functions Jn(x) and Yn(x), first and second derivatives. Licensing: This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, they give permission to incorporate this routine into a user program provided that the copyright is acknowledged. Modified: 08 June 2025 Author: Original FORTRAN77 version by Shanjie Zhang, Jianming Jin. This version by John Burkardt. Reference: Shanjie Zhang, Jianming Jin, Computation of Special Functions, Wiley, 1996, ISBN: 0-471-11963-6, LC: QA351.C45. Input: integer N, the order. real X, the argument. Output: real BJN, DJN, FJN, BYN, DYN, FYN, the values of Jn(x), Jn'(x), Jn"(x), Yn(x), Yn'(x), Yn"(x). */ { double bj[102]; double bs; double by[102]; double e0; double ec; double f; double f0; double f1; int k; int m; int mt; int nt; double s1; double su; for ( nt = 1; nt <= 900; nt++ ) { mt = ( int ) ( 0.5 * log10 ( 6.28 * nt ) - nt * log10 ( 1.36 * fabs ( x ) / nt ) ); if ( 20 < mt ) { break; } } m = nt; bs = 0.0; f0 = 0.0; f1 = 1.0E-35; su = 0.0; for ( k = m; 0 <= k; k-- ) { f = 2.0 * ( k + 1.0 ) * f1 / x - f0; if ( k <= n + 1 ) { bj[k] = f; } if ( ( k % 2 ) == 0 ) { bs = bs + 2.0 * f; if ( k != 0 ) { su = su + pow ( -1.0, k / 2 ) * f / k; } } f0 = f1; f1 = f; } for ( k = 0; k <= n + 1; k++ ) { bj[k] = bj[k] / ( bs - f ); } *bjn = bj[n]; ec = 0.5772156649015329; e0 = 0.3183098861837907; s1 = 2.0 * e0 * ( log ( x / 2.0 ) + ec ) * bj[0]; f0 = s1 - 8.0 * e0 * su / ( bs - f ); f1 = ( bj[1] * f0 - 2.0 * e0 / x ) / bj[0]; by[0] = f0; by[1] = f1; for ( k = 2; k <= n + 1; k++ ) { f = 2.0 * ( k - 1.0 ) * f1 / x - f0; by[k] = f; f0 = f1; f1 = f; } *byn = by[n]; *djn = - bj[n+1] + n * bj[n] / x; *dyn = - by[n+1] + n * by[n] / x; *fjn = ( n * n / ( x * x ) - 1.0 ) * *bjn - *djn / x; *fyn = ( n * n / ( x * x ) - 1.0 ) * *byn - *dyn / x; return; } /******************************************************************************/ void jyzo ( int n, int nt, double rj0[], double rj1[], double ry0[], double ry1[] ) /******************************************************************************/ /* Purpose: jyzo() computes the zeros of Bessel functions Jn(x), Yn(x) and derivatives. Licensing: This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, they give permission to incorporate this routine into a user program provided that the copyright is acknowledged. Modified: 08 June 2025 Author: Original FORTRAN77 version by Shanjie Zhang, Jianming Jin. This version by John Burkardt. Reference: Shanjie Zhang, Jianming Jin, Computation of Special Functions, Wiley, 1996, ISBN: 0-471-11963-6, LC: QA351.C45. Input: integer N, the order of the Bessel functions. integer NT, the number of zeros. Output: real RJ0(NT), RJ1(NT), RY0(NT), RY1(NT), the zeros of Jn(x), Jn'(x), Yn(x), Yn'(x). */ { double bjn; double byn; double djn; double dyn; double fjn; double fyn; int l; double n_r8; double x; double x0; n_r8 = ( double ) ( n ); if ( n <= 20 ) { x = 2.82141 + 1.15859 * n_r8; } else { x = n + 1.85576 * pow ( n_r8, 0.33333 ) + 1.03315 / pow ( n_r8, 0.33333 ); } l = 0; while ( true ) { x0 = x; jyndd ( n, x, &bjn, &djn, &fjn, &byn, &dyn, &fyn ); x = x - bjn / djn; if ( 1.0E-09 < fabs ( x - x0 ) ) { continue; } l = l + 1; rj0[l-1] = x; x = x + 3.1416 + ( 0.0972 + 0.0679 * n_r8 - 0.000354 * n_r8 * n_r8 ) / l; if ( nt <= l ) { break; } } if ( n <= 20 ) { x = 0.961587 + 1.07703 * n_r8; } else { x = n_r8 + 0.80861 * pow ( n_r8, 0.33333 ) + 0.07249 / pow ( n_r8, 0.33333 ); } if ( n == 0 ) { x = 3.8317; } l = 0; while ( true ) { x0 = x; jyndd ( n, x, &bjn, &djn, &fjn, &byn, &dyn, &fyn ); x = x - djn / fjn; if ( 1.0E-09 < fabs ( x - x0 ) ) { continue; } l = l + 1; rj1[l-1] = x; x = x + 3.1416 + ( 0.4955 + 0.0915 * n_r8 - 0.000435 * n_r8 * n_r8 ) / l; if ( nt <= l ) { break; } } if ( n <= 20 ) { x = 1.19477 + 1.08933 * n_r8; } else { x = n_r8 + 0.93158 * pow ( n_r8, 0.33333 ) + 0.26035 / pow ( n_r8, 0.33333 ); } l = 0; while ( true ) { x0 = x; jyndd ( n, x, &bjn, &djn, &fjn, &byn, &dyn, &fyn ); x = x - byn / dyn; if ( 1.0E-09 < fabs ( x - x0 ) ) { continue; } l = l + 1; ry0[l-1] = x; x = x + 3.1416 + ( 0.312 + 0.0852 * n_r8 - 0.000403 * n_r8 * n_r8 ) / l; if ( nt <= l ) { break; } } if ( n <= 20 ) { x = 2.67257 + 1.16099 * n_r8; } else { x = n_r8 + 1.8211 * pow ( n_r8, 0.33333 ) + 0.94001 / pow ( n_r8, 0.33333 ); } l = 0; while ( true ) { x0 = x; jyndd ( n, x, &bjn, &djn, &fjn, &byn, &dyn, &fyn ); x = x - dyn / fyn; if ( 1.0E-09 < fabs ( x - x0 ) ) { continue; } l = l + 1; ry1[l-1] = x; x = x + 3.1416 + ( 0.197 + 0.0643 * n_r8 -0.000286 * n_r8 * n_r8 ) / l; if ( nt <= l ) { break; } } return; }