/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
package org.opensourcephysics.numerics.specialfunctions;
import java.util.HashMap;
import java.util.Map;
import org.opensourcephysics.numerics.Function;
/*
* Bessel function of order n by J E Hasbun for the OSP library.
* Ref : http : //www.alglib.net/specialfunctions
*
* Using the Algorithm of Stephen L.Moshier Returns the Bessel function
* of order n, where n is a (possibly negative) integer. The ratio of jn(x)
* to j0(x) is computed by backward recurrence.First the ratio jn / jn -
* is found by a continued fraction expansion.Then the recurrence
* relating successive orders is applied until j0 or j1 is reached.
* If n = 0 or 1 the routine for j0 or j1 is called directly.
* ACCURACY - Absolute error :
* arithmetic range #trials peak rms
* IEEE 0, 3050004.4e-167.9e-17
* Not suitable for large n or x. Use jv() (fractional order) instead.
*
* OSP getFunction(n) method added by W. Christian in order to use the OSP Function interface.
*
* @author Javier E Hasbun 2008.
* @author Wolfgang Christian
*/
public class Bessel {
static final Map<Integer, Function> functionMap = new HashMap<Integer, Function>();
static final Map<Integer, Function> derivativeMap = new HashMap<Integer, Function>();
/**
* Gets the Bessel function with the given order.
*/
public static synchronized Function getFunction(int n) {
if(n<0) {
throw new IllegalArgumentException(Messages.getString("Bessel.0.neg_order")); //$NON-NLS-1$
}
Function f = functionMap.get(n);
if(f!=null) {
return f;
}
f = new BesselFunction(n);
functionMap.put(n, f); // polynomial was not in the list so add it.
return f;
}
/**
* Gets the derivative of the Bessel function with the given order.
*/
public static synchronized Function getDerivative(int n) {
if(n<0) {
throw new IllegalArgumentException(Messages.getString("Bessel.1.neg_order")); //$NON-NLS-1$
}
Function f = derivativeMap.get(n);
if(f!=null) {
return f;
}
f = new BesselDerivative(n);
derivativeMap.put(n, f); // polynomial was not in the list so add it.
return f;
}
/**
* Computes the Bessel function of order n at x.
* @param n
* @param x
* @return
*/
public static double besseln(int n, double x) {
int sg, k;
double y, tmp, pk, xk, pkm1, r, pkm2;
if(n<0) {
n = -n;
if(n%2==0) {
sg = 1;
} else {
sg = -1;
}
} else {
sg = 1;
}
if(x<0) {
if(n%2!=0) {
sg = -sg;
}
x = -x;
}
if(n==0) {
y = sg*bessel0(x);
return y;
}
if(n==1) {
y = sg*bessel1(x);
return y;
}
if(n==2) {
if(x==0) {
y = 0;
} else {
y = sg*(2.0*bessel1(x)/x-bessel0(x));
}
return y;
}
if(x<1.e-12) {
y = 0;
return y;
}
k = 53;
pk = 2*(n+k);
tmp = pk;
xk = x*x;
while(k!=0) {
pk = pk-2.0;
tmp = pk-xk/tmp;
k = k-1;
}
tmp = x/tmp;
pk = 1.0;
pkm1 = 1.0/tmp;
k = n-1;
r = 2*k;
while(k!=0) {
pkm2 = (pkm1*r-pk*x)/x;
pk = pkm1;
pkm1 = pkm2;
r = r-2.0;
k = k-1;
}
if(Math.abs(pk)>Math.abs(pkm1)) {
tmp = bessel1(x)/pk;
} else {
tmp = bessel0(x)/pkm1;
}
return sg*tmp;
}
/**
* Computes the derivative of the Bessel function of order n at x.
* @param n
* @param x
* @return
*/
public static double besselnDerivative(int n, double x) {
int m, qm, qp, nm, np;
double bjn, bjnm, bjnp;
m = n;
if(n==0) {
bjn = besseln(1, x);
return -bjn;
}
qm = 1;
qp = 1;
nm = m-1;
np = m+1;
if(m<0) {
if(nm<0) {
nm = -nm;
qm = -1;
}
if(np<0) {
np = -np;
qp = -1;
}
}
bjnm = besseln(nm, x);
bjnp = besseln(np, x);
bjnm = Math.pow(qm, nm)*bjnm;
bjnp = Math.pow(qp, np)*bjnp;
return(bjnm-bjnp)/2.;
}
/**
* Computes nt zeroes of the n-th order Bessel function
* @param n
* @param nt
* @return
*/
public static double[] besselnZeros(int n, int nt) {
int l;
double x, x0, bjn, djn;
double[] rj0 = new double[nt];
if(n<0) {
n = -n;
}
if(n<=20) {
x = 2.82141+1.15859*n;
} else {
x = n+1.85576*Math.pow(n, 0.33333)+1.03315/Math.pow(n, 0.33333);
}
l = 0;
while(true) {
while(true) {
x0 = x;
bjn = besseln(n, x);
djn = besselnDerivative(n, x);
x = x-bjn/djn; // Newton-Raphson step
if(!(Math.abs(x-x0)>1.0e-6)) {
break;
}
}
rj0[l] = x;
l = l+1;
x = x+Math.PI+(0.0972+0.0679*n-0.000354*Math.pow(n, 2))/l;
if(!(l<nt)) {
break;
}
}
return rj0;
}
/*
* Comutes the Bessel function of order 0 at x.
* Ref : http : //www.alglib.net/specialfunctions
* Using the Algorithm of Stephen L.Moshier
* The domain is divided into the intervals[0, 5] and (5, infinity).
* In the first interval the following rational approximation is used:
* 2 2
* (w - r ) (w - r ) P (w) / Q (w)
* 1 2 3 8
* where w = x^2 and the two r 'sare zeros of the function. In the second
* interval, the Hankel asymptotic expansion is employed with two
* rational functions of degree 6 / 6 and 7 / 7.
*/
public static double bessel0(double x) {
double nn, pzero, qzero, xsq, p1, q1, y;
double[] zz = new double[2];
double[] p = {26857.86856980014981415848441, -40504123.71833132706360663322, 25071582855.36881945555156435, -8085222034853.793871199468171, 1434354939140344.111664316553, -136762035308817138.6865416609, 6382059341072356562.289432465, -117915762910761053603.8440800, 493378725179413356181.6813446};
double[] q = {1.0, 1363.063652328970604442810507, 1114636.098462985378182402543, 669998767.2982239671814028660, 312304311494.1213172572469442, 112775673967979.8507056031594, 30246356167094626.98627330784, 5428918384092285160.200195092, 493378725179413356211.3278438};
if(x<0) {
x = -x;
}
if(x>8.0) {
zz = besselasympt0(x);
pzero = zz[0];
qzero = zz[1];
nn = x-Math.PI/4;
y = Math.sqrt(2/Math.PI/x)*(pzero*Math.cos(nn)-qzero*Math.sin(nn));
return y;
}
xsq = x*x;
p1 = p[0];
for(int i = 1; i<9; i++) {
p1 = p[i]+p1*xsq;
}
q1 = q[0];
for(int i = 1; i<9; i++) {
q1 = q[i]+q1*xsq;
}
return p1/q1;
}
/*
* Computes the Bessel function of order 1.
* Ref : http : //www.alglib.net/specialfunctions
* Using the Algorithm of Stephen L.Moshier
* The domain is divided into the intervals[0, 8] and (8, infinity).
* In the first interval a 24term Chebyshev expansion is used.In the
* second, the asymptotic trigonometric representation is employed
* using two rational functions of degree 5 / 5.
*/
public static double bessel1(double x) {
double s, pzero, qzero, nn, p1, q1, y, xsq;
double[] zz = new double[2];
double[] p = {2701.122710892323414856790990, -4695753.530642995859767162166, 3413234182.301700539091292655, -1322983480332.126453125473247, 290879526383477.5409737601689, -35888175699101060.50743641413, 2316433580634002297.931815435, -66721065689249162980.20941484, 581199354001606143928.050809};
double[] q = {1.0, 1606.931573481487801970916749, 1501793.594998585505921097578, 1013863514.358673989967045588, 524371026216.7649715406728642, 208166122130760.7351240184229, 60920613989175217.46105196863, 11857707121903209998.37113348, 1162398708003212287858.529400};
s = Math.signum(x);
if(x<0) {
x = -x;
}
if(x>8.0) {
zz = besselasympt1(x);
pzero = zz[0];
qzero = zz[1];
nn = x-3*Math.PI/4;
y = Math.sqrt(2/Math.PI/x)*(pzero*Math.cos(nn)-qzero*Math.sin(nn));
if(s<0) {
y = -y;
}
return y;
}
xsq = x*x;
p1 = p[0];
for(int i = 1; i<9; i++) {
p1 = p[i]+p1*xsq;
}
q1 = q[0];
for(int i = 1; i<9; i++) {
q1 = q[i]+q1*xsq;
}
return s*x*p1/q1;
}
private static double[] besselasympt0(double x) {
double xsq, p2, q2, p3, q3, pzero, qzero;
double[] zz = new double[2];
double[] p = {0.0, 2485.271928957404011288128951, 153982.6532623911470917825993, 2016135.283049983642487182349, 8413041.456550439208464315611, 12332384.76817638145232406055, 5393485.083869438325262122897};
double[] q = {1.0, 2615.700736920839685159081813, 156001.7276940030940592769933, 2025066.801570134013891035236, 8426449.050629797331554404810, 12338310.22786324960844856182, 5393485.083869438325560444960};
double[] pp = {0.0, -4.887199395841261531199129300, -226.2630641933704113967255053, -2365.956170779108192723612816, -8239.066313485606568803548860, -10381.41698748464093880530341, -3984.617357595222463506790588};
double[] qq = {1.0, 408.7714673983499223402830260, 15704.89191515395519392882766, 156021.3206679291652539287109, 533291.3634216897168722255057, 666745.4239319826986004038103, 255015.5108860942382983170882};
xsq = 64.0/x/x;
p2 = p[0];
for(int i = 1; i<7; i++) {
p2 = p[i]+p2*xsq;
}
q2 = q[0];
for(int i = 1; i<7; i++) {
q2 = q[i]+q2*xsq;
}
p3 = pp[0];
for(int i = 1; i<7; i++) {
p3 = pp[i]+p3*xsq;
}
q3 = qq[0];
for(int i = 1; i<7; i++) {
q3 = qq[i]+q3*xsq;
}
pzero = p2/q2;
qzero = 8*p3/q3/x;
zz[0] = pzero;
zz[1] = qzero;
return zz;
}
private static double[] besselasympt1(double x) {
double xsq, p2, q2, p3, q3, pzero, qzero;
double[] zz = new double[2];
double[] p = {-1611.616644324610116477412898, -109824.0554345934672737413139, -1523529.351181137383255105722, -6603373.248364939109255245434, -9942246.505077641195658377899, -4435757.816794127857114720794};
double[] q = {1.0, -1455.009440190496182453565068, -107263.8599110382011903063867, -1511809.506634160881644546358, -6585339.479723087072826915069, -9934124.389934585658967556309, -4435757.816794127856828016962};
double[] pp = {35.26513384663603218592175580, 1706.375429020768002061283546, 18494.26287322386679652009819, 66178.83658127083517939992166, 85145.16067533570196555001171, 33220.91340985722351859704442};
double[] qq = {1.0, 863.8367769604990967475517183, 37890.22974577220264142952256, 400294.4358226697511708610813, 1419460.669603720892855755253, 1819458.042243997298924553839, 708712.8194102874357377502472};
xsq = 64.0/x/x;
p2 = p[0];
for(int i = 1; i<6; i++) {
p2 = p[i]+p2*xsq;
}
q2 = q[0];
for(int i = 1; i<7; i++) {
q2 = q[i]+q2*xsq;
}
p3 = pp[0];
for(int i = 1; i<6; i++) {
p3 = pp[i]+p3*xsq;
}
q3 = qq[0];
for(int i = 1; i<7; i++) {
q3 = qq[i]+q3*xsq;
}
pzero = p2/q2;
qzero = 8*p3/q3/x;
zz[0] = pzero;
zz[1] = qzero;
return zz;
}
/**
* Computes the Bessel function of order n.
* @author Wolfgang Christian
*/
static class BesselFunction implements Function {
int n;
BesselFunction(int n) {
this.n = n;
}
/**
* Evaluates the Bessel function.
*/
public double evaluate(final double x) {
if(n==0) {
return Bessel.bessel0(x);
}
if(n==1) {
return Bessel.bessel1(x);
}
return Bessel.besseln(n, x);
}
}
/**
* Computes the derivative of the Bessel function or order n.
* @author Wolfgang Christian
*/
static class BesselDerivative implements Function {
int n;
BesselDerivative(int n) {
this.n = n;
}
/**
* Evaluates the derivative of the Bessel function.
*/
public double evaluate(final double x) {
return Bessel.besselnDerivative(n, x);
}
}
}
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must also be released
* under the GNU GPL license.
*
* This software is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/