/
bessels.ts
164 lines (143 loc) · 5.02 KB
/
bessels.ts
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
//================================================================
/**
* @packageDocumentation
* @module std
*/
//================================================================
import { MathUtil } from "../../internal/numeric/MathUtil";
import { InvalidArgument } from "../../exception/InvalidArgument";
import { tgamma } from "./gamma";
const INFINITY = 100; // (1 / 30!) is nearby 0.
/*================================================================
ORIGINAL FUNCTIONS
- CYLINDRICAL
- SPHERICAL
==================================================================
FIRST KIND
--------------------------------------------------------------- */
/**
* Bessel function of the 1st kind.
*
* @reference https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J.CE.B1
*/
export function cyl_bessel_j(n: number, x: number): number
{
// VALIDATION
if (x < 0 && Math.floor(n) !== n)
throw new InvalidArgument(`Error on std.cyl_bessel_j(): n must be integer when x is negative -> (n = ${n}, x = ${x}).`);
else if (x === 0 && n !== 0)
throw new InvalidArgument(`Error on std.cyl_bessel_j(): n must be zero when x is zero -> (n = ${n}, x = ${x}).`);
// COMPUTATION
if (n === Math.floor(n))
return _J_int(n, x);
else
return _J_positive(n, x);
}
/**
* Bessel function of the 2nd kind.
*
* @reference https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y.CE.B1
*/
export function cyl_neumann(v: number, x: number): number
{
if (x <= 0)
throw new InvalidArgument(`Error on std.cyl_neumann(): x must be greater than zero -> (x = ${x}).`);
let numerator: number = cyl_bessel_j(v, x) * Math.cos(v*Math.PI) - cyl_bessel_j(-v, x)
let denominator: number = Math.sin(v * Math.PI);
return numerator / denominator;
}
function _J_int(n: number, x: number): number
{
if (n < 0)
return Math.pow(-1, n) * _J_positive(-n, x);
else
return _J_positive(n, x);
}
function _J_positive(v: number, x: number): number
{
let sigma: number = MathUtil.sigma(function (k: number): number
{
let ret: number = Math.pow(-1, k) * Math.pow(x/2, v + 2*k);
ret /= MathUtil.factorial(k) * tgamma(v + k + 1);
return ret;
}, 0, INFINITY);
return sigma;
}
/* ---------------------------------------------------------------
SPHERICAL
--------------------------------------------------------------- */
/**
* Spherical Bessel function of the 1st kind.
*
* @reference https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn.2C_yn
*/
export function sph_bessel(n: number, x: number): number
{
return Math.sqrt(Math.PI / (2*x)) * cyl_bessel_j(n+.5, x);
}
/**
* Spherical Bessel function of the 2nd kind.
*
* @reference https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn.2C_yn
*/
export function sph_neumann(n: number, x: number): number
{
let ret: number = Math.sqrt(Math.PI / (2*x));
ret *= cyl_neumann(n + .5, x);
return ret;
}
/*================================================================
REQGULAR MODIFIED
- FIRST KIND
- SECOND KIND
==================================================================
FIRST KIND
--------------------------------------------------------------- */
/**
* Modified cylindrical Bessel function of the 1st kind.
*
* @reference https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I.CE.B1_.2C_K.CE.B1
*/
export function cyl_bessel_i(n: number, x: number): number
{
// VALIDATION
if (x < 0 && Math.floor(n) !== n)
throw new InvalidArgument(`Error on std.cyl_bessel_i(): n must be integer when x is negative -> (n = ${n}, x = ${x}).`);
else if (x === 0 && n !== 0)
throw new InvalidArgument(`Error on std.cyl_bessel_i(): n must be zero when x is zero -> (n = ${n}, x = ${x}).`);
// COMPUTATION
if (n === .5)
return Math.sqrt(2.0 / (Math.PI*x)) * Math.sinh(x);
else
return _Bessel_i(n, x);
}
function _Bessel_i(v: number, x: number): number
{
return MathUtil.sigma(function (k: number): number
{
let numerator: number = Math.pow(x / 2, v + 2*k);
let denominator: number = MathUtil.factorial(k) * tgamma(v + k + 1);
return numerator / denominator;
}, 0, INFINITY);
}
/* ---------------------------------------------------------------
SECOND KIND
--------------------------------------------------------------- */
/**
* Modified cylindrical Bessel function of the 2nd kind.
*
* @reference https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I.CE.B1_.2C_K.CE.B1
*/
export function cyl_bessel_k(n: number, x: number): number
{
if (x <= 0)
throw new InvalidArgument(`Error on std.cyl_bessel_k(): requires x > 0 -> (x = ${x}).`);
return _Bessel_k(n, x);
}
function _Bessel_k(v: number, z: number): number
{
let ret: number = Math.PI / 2;
ret *= cyl_bessel_i(-v, z) - cyl_bessel_i(v, z);
ret /= Math.sin(v*Math.PI);
return ret;
}