(* :Name: NumericalMath`BesselZeros` *)

(* :Context: NumericalMath`BesselZeros` *)

(* :Title: Zeros of Bessel Functions *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
This package provides functions that find the zeros of Bessel functions.
*)

(* :Copyright: Copyright 1992-1999,  Wolfram Research, Inc.
*)

(* :History:
	Jerry B. Keiper added extrapolation guesses, September 1994.
	Original version by Jerry B. Keiper, February 1992.
*)

(* :Mathematica Version: 2.1 *)

(* :Source: Abramowitz and Stegun, Handbook of Mathematical Functions,
	equations 9.5.12, 9.5.13, 9.5.27-33. *)

(* :Limitations:
*)

BeginPackage["NumericalMath`BesselZeros`"]

Unprotect[BesselJZeros, BesselYZeros, BesselJPrimeZeros, BesselYPrimeZeros,
	BesselJYJYZeros, BesselJPrimeYPrimeJPrimeYPrimeZeros,
	BesselJPrimeYJYPrimeZeros]; 

BesselJZeros::usage =
"BesselJZeros[nu, n] gives a list of the first n zeros of the order nu
BesselJ function. BesselJZeros[nu, {m, n}] gives a list of the mth through
the nth zeros."

BesselYZeros::usage =
"BesselYZeros[nu, n] gives a list of the first n zeros of the order nu
BesselY function. BesselYZeros[nu, {m, n}] gives a list of the mth through
the nth zeros."

BesselJPrimeZeros::usage =
"BesselJPrimeZeros[nu, n] gives a list of the first n zeros of the
derivative of the order nu BesselJ function. BesselJPrimeZeros[nu, {m, n}]
gives a list of the mth through the nth zeros."

BesselYPrimeZeros::usage =
"BesselYPrimeZeros[nu, n] gives a list of the first n zeros of the
derivative of the order nu BesselY function. BesselYPrimeZeros[nu, {m, n}]
gives a list of the mth through the nth zeros."

BesselJYJYZeros::usage =
"BesselJYJYZeros[nu, l, n] gives a list of the first n zeros of
BesselJ[nu, z] BesselY[nu, l z] - BesselJ[nu, l z] BesselY[nu, z].
BesselJYJYZeros[nu, l, {m, n}] gives a list of the mth through the
nth zeros."

BesselJPrimeYPrimeJPrimeYPrimeZeros::usage =
"BesselJPrimeYPrimeJPrimeYPrimeZeros[nu, l, n] gives a list of the
first n zeros of
BesselJ'[nu, z] BesselY'[nu, l z] - BesselJ'[nu, l z] BesselY'[nu, z].
BesselJPrimeYPrimeJPrimeYPrimeZeros[nu, l, {m, n}] gives a list of
the mth through the nth zeros."

BesselJPrimeYJYPrimeZeros::usage =
"BesselJPrimeYJYPrimeZeros[nu, l, n] gives a list of the first n zeros of
BesselJ'[nu, z] BesselY[nu, l z] - BesselJ[nu, l z] BesselY'[nu, z].
BesselJPrimeYJYPrimeZeros[nu, l, {m, n}] gives a list of the mth
through the nth zeros."

Options[BesselJZeros] =
Options[BesselYZeros] =
Options[BesselJPrimeZeros] =
Options[BesselYPrimeZeros] =
Options[BesselJYJYZeros] =
Options[BesselJPrimeYPrimeJPrimeYPrimeZeros] =
Options[BesselJPrimeYJYPrimeZeros] =
	{WorkingPrecision -> $MachinePrecision, AccuracyGoal -> Automatic}

Begin["NumericalMath`BesselZeros`Private`"]

x$fr;  (* variable to be used by FindRoot *)

vnl = ((TrueQ[NumberQ[#] && (# >= 0)])&)    (* validity test for nu and l *)

mu[nu_] := 4 nu^2;

(* using rational interpolation, predict the predecessor of
	the first element in x *)
pred[x_List] :=
    Module[{b = InterpolatingPolynomial[Drop[x,1]-Drop[x,-1],x$fr] /. x$fr->0},
	If[x[[1]] <= 1.2 b,
		{x$fr, 0.000001, 0.3 x[[1]]},
		{x$fr, x[[1]] - 1.2 b, x[[1]] - 0.8 b}]
	]
	
BesselJZeros[nu_?vnl, {m_Integer, n_Integer}, opts___] :=
    Module[{ans}, ans /; (ans = bz[nu, m, n, BesselJZeros, opts]) =!= $Failed];

BesselJZeros[nu_?vnl, n_Integer, opts___] :=
    Module[{ans}, ans /; (ans = bz[nu, 1, n, BesselJZeros, opts]) =!= $Failed];

BesselYZeros[nu_?vnl, {m_Integer, n_Integer}, opts___] :=
    Module[{ans}, ans /; (ans = bz[nu, m, n, BesselYZeros, opts]) =!= $Failed];

BesselYZeros[nu_?vnl, n_Integer, opts___] :=
    Module[{ans}, ans /; (ans = bz[nu, 1, n, BesselYZeros, opts]) =!= $Failed];

BesselJPrimeZeros[nu_?vnl, {m_Integer, n_Integer}, opts___] :=
    Module[{ans},
	ans /; (ans = bz[nu, m, n, BesselJPrimeZeros, opts]) =!= $Failed];

BesselJPrimeZeros[nu_?vnl, n_Integer, opts___] :=
    Module[{ans},
	ans /; (ans = bz[nu, 1, n, BesselJPrimeZeros, opts]) =!= $Failed];

BesselYPrimeZeros[nu_?vnl, {m_Integer, n_Integer}, opts___] :=
    Module[{ans},
	ans /; (ans = bz[nu, m, n, BesselYPrimeZeros, opts]) =!= $Failed];

BesselYPrimeZeros[nu_?vnl, n_Integer, opts___] :=
    Module[{ans},
	ans /; (ans = bz[nu, 1, n, BesselYPrimeZeros, opts]) =!= $Failed];

bz[nu_, m_, n_, f_, opts___] :=
    Module[{i, wp, ag, ff, his, xp, nn, d},
	wp = WorkingPrecision /. {opts} /. Options[f];
	ag = AccuracyGoal /. {opts} /. Options[f];
	If[ag === Automatic, ag = wp - 6];
	If[!TrueQ[0 < wp && 0 < ag < wp], Return[$Failed]];
	ff = bzfunc[f, nu];
	d = 1;
	nn = n+1;
	While[!NumberQ[guess[f, nn, nu]], nn += Ceiling[d *= 1.2]];
	his = Table[xp = guess[f, i, nu];
		x$fr /. FindRoot[ff[x$fr], {x$fr, xp-0.02, xp+0.02},
			WorkingPrecision -> wp, AccuracyGoal -> ag,
			MaxIterations -> 35],
		{i, nn, nn+3}];
	If[!VectorQ[his, NumberQ], Return[$Failed]];
	While[!OrderedQ[his],
		nn += Ceiling[d *= 1.2];
		his = Table[xp = guess[f, i, nu];
			x$fr /. FindRoot[ff[x$fr], {x$fr, xp-0.02, xp+0.02},
				WorkingPrecision -> wp, AccuracyGoal -> ag,
				MaxIterations -> 35],
			{i, nn, nn+3}];
		If[!VectorQ[his, NumberQ], Return[$Failed]];
		];
	xp = Drop[his,1]-Drop[his,-1];
	While[Max[Abs[(Drop[xp,1]-Drop[xp,-1])/xp[[1]]]] > 0.25,
		nn += Ceiling[d *= 1.2];
		his = Table[xp = guess[f, i, nu];
			x$fr /. FindRoot[ff[x$fr], {x$fr, xp-0.02, xp+0.02},
				WorkingPrecision -> wp, AccuracyGoal -> ag,
				MaxIterations -> 35],
			{i, nn, nn+3}];
		If[!VectorQ[his, NumberQ], Return[$Failed]];
		xp = Drop[his,1]-Drop[his,-1]
		];
	Do[ xp = x$fr /. FindRoot[ff[x$fr], Evaluate[pred[his]],
				WorkingPrecision -> $MachinePrecision,
				AccuracyGoal -> 6, MaxIterations -> 35];
	    If[!NumberQ[xp], Return[$Failed]];
	    his = RotateRight[his];
	    his[[1]] = xp, {nn-n-1}];
	Reverse[Table[
	    xp = x$fr /. FindRoot[ff[x$fr], Evaluate[pred[his]],
				WorkingPrecision -> wp, AccuracyGoal -> ag,
				MaxIterations -> 35];
	    If[!NumberQ[xp], Return[$Failed]];
	    his = RotateRight[his];
	    his[[1]] = xp;
	    xp, {n-m+1}]]
	]

BesselJYJYZeros[nu_?vnl, l_?vnl, {m_Integer, n_Integer}, opts___] :=
    Module[{ans},
	ans /; (ans = bcz[nu, l, m, n, BesselJYJYZeros, opts]) =!= $Failed];

BesselJYJYZeros[nu_?vnl, l_?vnl, n_Integer, opts___] :=
    Module[{ans},
	ans /; (ans = bcz[nu, l, 1, n, BesselJYJYZeros, opts]) =!= $Failed];

BesselJPrimeYPrimeJPrimeYPrimeZeros[nu_?vnl, l_?vnl,
					{m_Integer, n_Integer}, opts___] :=
    Module[{ans}, ans /; (ans = bcz[nu, l, m, n,
	BesselJPrimeYPrimeJPrimeYPrimeZeros, opts]) =!= $Failed];

BesselJPrimeYPrimeJPrimeYPrimeZeros[nu_?vnl, l_?vnl, n_Integer, opts___] :=
    Module[{ans}, ans /; (ans = bcz[nu, l, 1, n,
	BesselJPrimeYPrimeJPrimeYPrimeZeros, opts]) =!= $Failed];

BesselJPrimeYJYPrimeZeros[nu_?vnl, l_?vnl, {m_Integer, n_Integer}, opts___] :=
    Module[{ans}, ans /;
	(ans = bcz[nu, l, m, n, BesselJPrimeYJYPrimeZeros, opts]) =!= $Failed];

BesselJPrimeYJYPrimeZeros[nu_?vnl, l_?vnl, n_Integer, opts___] :=
    Module[{ans}, ans /;
	(ans = bcz[nu, l, 1, n, BesselJPrimeYJYPrimeZeros, opts]) =!= $Failed];

bcz[nu_, l_, m_, n_, f_, opts___] :=
    Module[{i, wp, ag, ff, his, xp, nn, d},
	wp = WorkingPrecision /. {opts} /. Options[f];
	ag = AccuracyGoal /. {opts} /. Options[f];
	If[ag === Automatic, ag = wp - 6];
	If[!TrueQ[0 < wp && 0 < ag < wp && l != 1], Return[$Failed]];
	ff = bczfunc[f, nu, l];
	d = 1;
	nn = n+1;
	While[!NumberQ[cguess[f, nn, nu, l]], nn += Ceiling[d *= 1.2]];
	his = Table[xp = cguess[f, i, nu, l];
		x$fr /. FindRoot[ff[x$fr], {x$fr, xp-0.02, xp+0.02},
			WorkingPrecision -> wp, AccuracyGoal -> ag,
			MaxIterations -> 35],
		{i, nn, nn+3}];
	If[!VectorQ[his, NumberQ], Return[$Failed]];
	While[!OrderedQ[his],
		nn += Ceiling[d *= 1.2];
		his = Table[xp = cguess[f, i, nu, l];
			x$fr /. FindRoot[ff[x$fr], {x$fr, xp-0.02, xp+0.02},
				WorkingPrecision -> wp, AccuracyGoal -> ag,
				MaxIterations -> 35],
			{i, nn, nn+3}];
		If[!VectorQ[his, NumberQ], Return[$Failed]];
		];
	xp = Drop[his,1]-Drop[his,-1];
	While[Max[Abs[(Drop[xp,1]-Drop[xp,-1])/xp[[1]]]] > 0.25,
		nn += Ceiling[d *= 1.2];
		his = Table[xp = cguess[f, i, nu, l];
			x$fr /. FindRoot[ff[x$fr], {x$fr, xp-0.02, xp+0.02},
				WorkingPrecision -> wp, AccuracyGoal -> ag,
				MaxIterations -> 35],
			{i, nn, nn+3}];
		If[!VectorQ[his, NumberQ], Return[$Failed]];
		xp = Drop[his,1]-Drop[his,-1]
		];
	Do[ xp = x$fr /. FindRoot[ff[x$fr], Evaluate[pred[his]],
				WorkingPrecision -> $MachinePrecision,
				AccuracyGoal -> 6, MaxIterations -> 35];
	    If[!NumberQ[xp], Return[$Failed]];
	    his = RotateRight[his];
	    his[[1]] = xp, {nn-n-1}];
	Reverse[Table[
	    xp = x$fr /. FindRoot[ff[x$fr], Evaluate[pred[his]],
				WorkingPrecision -> wp, AccuracyGoal -> ag,
				MaxIterations -> 35];
	    If[!NumberQ[xp], Return[$Failed]];
	    his = RotateRight[his];
	    his[[1]] = xp;
	    xp, {n-m+1}]]
	]

(* -------------------- J, Y, JPrime YPrime --------------------------- *)

beta[BesselJZeros, s_, nu_] := N[Pi] (s + nu/2 - 0.25);
beta[BesselYZeros, s_, nu_] := N[Pi] (s + nu/2 - 0.75);
beta[BesselJPrimeZeros, s_, nu_] := N[Pi] (s + nu/2 - 0.75);
beta[BesselYPrimeZeros, s_, nu_] := N[Pi] (s + nu/2 - 0.25);

guess[BesselJZeros, s_, nu_] := guess1[BesselJZeros, s, nu];
guess[BesselYZeros, s_, nu_] := guess1[BesselYZeros, s, nu];
guess[BesselJPrimeZeros, s_, nu_] := guess2[BesselJPrimeZeros, s, nu];
guess[BesselYPrimeZeros, s_, nu_] := guess2[BesselYPrimeZeros, s, nu];

guess1[f_, s_, nu_] :=		(* A&S 9.5.12 *)
    Module[{b = beta[f, s, nu], m = mu[nu], z},
	z = 64(m-1)(6949m^3-153855m^2+1585743m-6277237)/(105 (8b)^7);
	If[Abs[z] > 1/4, Return[$Failed]];
	z += 32(m-1)(83m^2-982m+3779)/(15 (8b)^5);
	z += (m-1)/(8b) + 4(m-1)(7m-31)/(3 (8b)^3);
	b-z
	]

guess2[f_, s_, nu_] :=		(* A&S 9.5.13 *)
    Module[{b = beta[f, s, nu], m = mu[nu], z},
	z = 64(6949m^4+296492m^3-1248002m^2+7414380m-5853627)/(105 (8b)^7);
	If[Abs[z] > 1/4, Return[$Failed]];
	z += 32(83m^3+2075m^2-3039m+3537)/(15 (8b)^5);
	z += (m+3)/(8b) + 4(7m^2+82m-9)/(3 (8b)^3);
	b-z
	]

bzfunc[BesselJZeros, nu_] := BesselJ[nu, #]&
bzfunc[BesselYZeros, nu_] := BesselY[nu, #]&
bzfunc[BesselJPrimeZeros, nu_] := (Evaluate[Derivative[0,1][BesselJ][nu, #]])&
bzfunc[BesselYPrimeZeros, nu_] := (Evaluate[Derivative[0,1][BesselY][nu, #]])&

(* -------------------- cross products --------------------------- *)

cguess[BesselJYJYZeros, s_, nu_, l_] := guess3[s, nu, l];
cguess[BesselJPrimeYPrimeJPrimeYPrimeZeros, s_, nu_, l_] := guess4[s, nu, l];
cguess[BesselJPrimeYJYPrimeZeros, s_, nu_, l_] := guess5[s, nu, l];

guess3[s_, nu_, ll_] :=		(* A&S 9.5.28 *)
    Module[{b, m = mu[nu], z, p, q, r, l},
	l = If[ll > 1, ll, 1/ll];
	b = N[Pi] s/(l-1);
	p = (m-1)/(8l);
	q = (m-1)(m-25)(l^3-1)/(6 (4l)^3 (l-1));
	r = (m-1)(m^2-114m+1073)(l^5-1)/(5 (4l)^5 (l-1));
	p = 1/(l-1);
	z = (r - 4 p q + 2 p^3)/b^5;
	If[Abs[z] > p/4, Return[$Failed]];
	z = b + p/b + (q - p^2)/b^3 + z;
	If[ll < 1, z *= l];
	z
	]

guess4[s_, nu_, ll_] :=		(* A&S 9.5.31 *)
    Module[{b, m = mu[nu], z, p, q, r, l},
	l = If[ll > 1, ll, 1/ll];
	b = N[Pi] (s-1)/(l-1);
	p = (m+3)/(8l);
	q = (m^2+46m-63)(l^3-1)/(6 (4l)^3 (l-1));
	r = (m^3+185m^2-2053m+1899)(l^5-1)/(5 (4l)^5 (l-1));
	p = 1/(l-1);
	z = (r - 4 p q + 2 p^3)/b^5;
	If[Abs[z] > p/4, Return[$Failed]];
	z = b + p/b + (q - p^2)/b^3 + z;
	If[ll < 1, z *= l];
	z
	]

guess5[s_, nu_, ll_] :=		(* A&S 9.5.33 *)
    Module[{b, m = mu[nu], z, p, q, r, l},
	l = If[ll > 1, ll, 1/ll];
	b = N[Pi] (s-1/2)/(l-1);
	p = ((m+3)l-(m-1))/(8l(l-1));
	q = ((m^2+46m-63)l^3-(m-1)(m-25))/(6 (4l)^3 (l-1));
	r = ((m^3+185m^2-2053m+1899)l^5-(m-1)(m^2-114m+1073))/(5 (4l)^5 (l-1));
	p = 1/(l-1);
	z = (r - 4 p q + 2 p^3)/b^5;
	If[Abs[z] > p/4, Return[$Failed]];
	z = b + p/b + (q - p^2)/b^3 + z;
	If[ll < 1, z *= l];
	z
	]

bczfunc[BesselJYJYZeros, nu_, l_] :=
	(BesselJ[nu, #] BesselY[nu, l #] - BesselJ[nu, l #] BesselY[nu, #])&

bczfunc[BesselJPrimeYPrimeJPrimeYPrimeZeros, nu_, l_] := (Evaluate[
	(Derivative[0,1][BesselJ][nu, #] Derivative[0,1][BesselY][nu, l #] -
	Derivative[0,1][BesselJ][nu, l #] Derivative[0,1][BesselY][nu, #])])&

bczfunc[BesselJPrimeYJYPrimeZeros, nu_, l_] := (Evaluate[
	(Derivative[0,1][BesselJ][nu, #] BesselY[nu, l #] -
	BesselJ[nu, l #] Derivative[0,1][BesselY][nu, #])])&

End[ ] (* NumericalMath`BesselZeros`Private` *)

Protect[BesselJZeros, BesselYZeros, BesselJPrimeZeros, BesselYPrimeZeros,
	BesselJYJYZeros, BesselJPrimeYPrimeJPrimeYPrimeZeros,
	BesselJPrimeYJYPrimeZeros];

EndPackage[] (* NumericalMath`BesselZeros` *)

(* Tests:
BesselJ[1, #]& /@ BesselJZeros[1, 3] 
BesselY[2, #]& /@ BesselYZeros[2, 4]
(Evaluate[Derivative[0,1][BesselJ][3, #]])& /@ BesselJPrimeZeros[3, 4]
(Evaluate[Derivative[0,1][BesselY][3, #]])& /@ BesselYPrimeZeros[3, 4]
(BesselJ[1, #] BesselY[1, 2 #] - BesselJ[1, 2 #] BesselY[1, #])& /@
	BesselJYJYZeros[1, 2, 3]
Evaluate[Derivative[0,1][BesselJ][1,#] Derivative[0,1][BesselY][1,2#] -
	Derivative[0,1][BesselJ][1,2#] Derivative[0,1][BesselY][1,#]]& /@
	BesselJPrimeYPrimeJPrimeYPrimeZeros[1, 2, 3]
Evaluate[Derivative[0,1][BesselJ][2,#] BesselY[2, 11/10 #] -
	BesselJ[2, 11/10 #] Derivative[0,1][BesselY][2,#]]& /@
	BesselJPrimeYJYPrimeZeros[2, 11/10, 3,
		WorkingPrecision -> 30, AccuracyGoal -> 20]
Evaluate[Derivative[0,1][BesselJ][2,#] BesselY[2, 1/10 #] -
	BesselJ[2, 1/10 #] Derivative[0,1][BesselY][2,#]]& /@
	BesselJPrimeYJYPrimeZeros[2, 1/10, 3]
(BesselJ[1, #] BesselY[1, 1/2 #] - BesselJ[1, 1/2 #] BesselY[1, #])& /@
	BesselJYJYZeros[1, 1/2, 3]
(BesselJ[1, #] BesselY[1, 9/10 #] - BesselJ[1, 9/10 #] BesselY[1, #])& /@
	BesselJYJYZeros[1, 9/10, 3]
Evaluate[Derivative[0,1][BesselJ][1,#] Derivative[0,1][BesselY][1,#/2] -
	Derivative[0,1][BesselJ][1,#/2] Derivative[0,1][BesselY][1,#]]& /@
	BesselJPrimeYPrimeJPrimeYPrimeZeros[1, 1/2, 3]

BesselJ[1, #]& /@ BesselJZeros[1, {3, 5}] 
BesselY[2, #]& /@ BesselYZeros[2, {4, 6}]
(Evaluate[Derivative[0,1][BesselJ][3, #]])& /@ BesselJPrimeZeros[3, {2, 4}]
(Evaluate[Derivative[0,1][BesselY][3, #]])& /@ BesselYPrimeZeros[3, {2, 4}]
(BesselJ[1, #] BesselY[1, 2 #] - BesselJ[1, 2 #] BesselY[1, #])& /@
	BesselJYJYZeros[1, 2, {3, 5}]
Evaluate[Derivative[0,1][BesselJ][1,#] Derivative[0,1][BesselY][1,2#] -
	Derivative[0,1][BesselJ][1,2#] Derivative[0,1][BesselY][1,#]]& /@
	BesselJPrimeYPrimeJPrimeYPrimeZeros[1, 2, {2, 3}]
Evaluate[Derivative[0,1][BesselJ][2,#] BesselY[2, 11/10 #] -
	BesselJ[2, 11/10 #] Derivative[0,1][BesselY][2,#]]& /@
	BesselJPrimeYJYPrimeZeros[2, 11/10, {707, 709}]
*)


