1441 |
ariadna |
1 |
<?php
|
|
|
2 |
|
|
|
3 |
namespace PhpOffice\PhpSpreadsheet\Calculation\Engineering;
|
|
|
4 |
|
|
|
5 |
use PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled;
|
|
|
6 |
use PhpOffice\PhpSpreadsheet\Calculation\Exception;
|
|
|
7 |
use PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
|
|
|
8 |
|
|
|
9 |
class BesselJ
|
|
|
10 |
{
|
|
|
11 |
use ArrayEnabled;
|
|
|
12 |
|
|
|
13 |
/**
|
|
|
14 |
* BESSELJ.
|
|
|
15 |
*
|
|
|
16 |
* Returns the Bessel function
|
|
|
17 |
*
|
|
|
18 |
* Excel Function:
|
|
|
19 |
* BESSELJ(x,ord)
|
|
|
20 |
*
|
|
|
21 |
* NOTE: The MS Excel implementation of the BESSELJ function is still not accurate, particularly for higher order
|
|
|
22 |
* values with x < -8 and x > 8. This code provides a more accurate calculation
|
|
|
23 |
*
|
|
|
24 |
* @param mixed $x A float value at which to evaluate the function.
|
|
|
25 |
* If x is nonnumeric, BESSELJ returns the #VALUE! error value.
|
|
|
26 |
* Or can be an array of values
|
|
|
27 |
* @param mixed $ord The integer order of the Bessel function.
|
|
|
28 |
* If ord is not an integer, it is truncated.
|
|
|
29 |
* If $ord is nonnumeric, BESSELJ returns the #VALUE! error value.
|
|
|
30 |
* If $ord < 0, BESSELJ returns the #NUM! error value.
|
|
|
31 |
* Or can be an array of values
|
|
|
32 |
*
|
|
|
33 |
* @return array|float|string Result, or a string containing an error
|
|
|
34 |
* If an array of numbers is passed as an argument, then the returned result will also be an array
|
|
|
35 |
* with the same dimensions
|
|
|
36 |
*/
|
|
|
37 |
public static function BESSELJ(mixed $x, mixed $ord): array|string|float
|
|
|
38 |
{
|
|
|
39 |
if (is_array($x) || is_array($ord)) {
|
|
|
40 |
return self::evaluateArrayArguments([self::class, __FUNCTION__], $x, $ord);
|
|
|
41 |
}
|
|
|
42 |
|
|
|
43 |
try {
|
|
|
44 |
$x = EngineeringValidations::validateFloat($x);
|
|
|
45 |
$ord = EngineeringValidations::validateInt($ord);
|
|
|
46 |
} catch (Exception $e) {
|
|
|
47 |
return $e->getMessage();
|
|
|
48 |
}
|
|
|
49 |
|
|
|
50 |
if ($ord < 0) {
|
|
|
51 |
return ExcelError::NAN();
|
|
|
52 |
}
|
|
|
53 |
|
|
|
54 |
$fResult = self::calculate($x, $ord);
|
|
|
55 |
|
|
|
56 |
return (is_nan($fResult)) ? ExcelError::NAN() : $fResult;
|
|
|
57 |
}
|
|
|
58 |
|
|
|
59 |
private static function calculate(float $x, int $ord): float
|
|
|
60 |
{
|
|
|
61 |
return match ($ord) {
|
|
|
62 |
|
|
|
63 |
1 => self::besselJ1($x),
|
|
|
64 |
default => self::besselJ2($x, $ord),
|
|
|
65 |
};
|
|
|
66 |
}
|
|
|
67 |
|
|
|
68 |
private static function besselJ0(float $x): float
|
|
|
69 |
{
|
|
|
70 |
$ax = abs($x);
|
|
|
71 |
|
|
|
72 |
if ($ax < 8.0) {
|
|
|
73 |
$y = $x * $x;
|
|
|
74 |
$ans1 = 57568490574.0 + $y * (-13362590354.0 + $y * (651619640.7 + $y * (-11214424.18 + $y
|
|
|
75 |
* (77392.33017 + $y * (-184.9052456)))));
|
|
|
76 |
$ans2 = 57568490411.0 + $y * (1029532985.0 + $y * (9494680.718 + $y * (59272.64853 + $y
|
|
|
77 |
* (267.8532712 + $y * 1.0))));
|
|
|
78 |
|
|
|
79 |
return $ans1 / $ans2;
|
|
|
80 |
}
|
|
|
81 |
|
|
|
82 |
$z = 8.0 / $ax;
|
|
|
83 |
$y = $z * $z;
|
|
|
84 |
$xx = $ax - 0.785398164;
|
|
|
85 |
$ans1 = 1.0 + $y * (-0.1098628627e-2 + $y * (0.2734510407e-4 + $y * (-0.2073370639e-5 + $y * 0.2093887211e-6)));
|
|
|
86 |
$ans2 = -0.1562499995e-1 + $y * (0.1430488765e-3 + $y * (-0.6911147651e-5 + $y
|
|
|
87 |
* (0.7621095161e-6 - $y * 0.934935152e-7)));
|
|
|
88 |
|
|
|
89 |
return sqrt(0.636619772 / $ax) * (cos($xx) * $ans1 - $z * sin($xx) * $ans2);
|
|
|
90 |
}
|
|
|
91 |
|
|
|
92 |
private static function besselJ1(float $x): float
|
|
|
93 |
{
|
|
|
94 |
$ax = abs($x);
|
|
|
95 |
|
|
|
96 |
if ($ax < 8.0) {
|
|
|
97 |
$y = $x * $x;
|
|
|
98 |
$ans1 = $x * (72362614232.0 + $y * (-7895059235.0 + $y * (242396853.1 + $y
|
|
|
99 |
* (-2972611.439 + $y * (15704.48260 + $y * (-30.16036606))))));
|
|
|
100 |
$ans2 = 144725228442.0 + $y * (2300535178.0 + $y * (18583304.74 + $y * (99447.43394 + $y
|
|
|
101 |
* (376.9991397 + $y * 1.0))));
|
|
|
102 |
|
|
|
103 |
return $ans1 / $ans2;
|
|
|
104 |
}
|
|
|
105 |
|
|
|
106 |
$z = 8.0 / $ax;
|
|
|
107 |
$y = $z * $z;
|
|
|
108 |
$xx = $ax - 2.356194491;
|
|
|
109 |
|
|
|
110 |
$ans1 = 1.0 + $y * (0.183105e-2 + $y * (-0.3516396496e-4 + $y * (0.2457520174e-5 + $y * (-0.240337019e-6))));
|
|
|
111 |
$ans2 = 0.04687499995 + $y * (-0.2002690873e-3 + $y * (0.8449199096e-5 + $y
|
|
|
112 |
* (-0.88228987e-6 + $y * 0.105787412e-6)));
|
|
|
113 |
$ans = sqrt(0.636619772 / $ax) * (cos($xx) * $ans1 - $z * sin($xx) * $ans2);
|
|
|
114 |
|
|
|
115 |
return ($x < 0.0) ? -$ans : $ans;
|
|
|
116 |
}
|
|
|
117 |
|
|
|
118 |
private static function besselJ2(float $x, int $ord): float
|
|
|
119 |
{
|
|
|
120 |
$ax = abs($x);
|
|
|
121 |
if ($ax === 0.0) {
|
|
|
122 |
return 0.0;
|
|
|
123 |
}
|
|
|
124 |
|
|
|
125 |
if ($ax > $ord) {
|
|
|
126 |
return self::besselj2a($ax, $ord, $x);
|
|
|
127 |
}
|
|
|
128 |
|
|
|
129 |
return self::besselj2b($ax, $ord, $x);
|
|
|
130 |
}
|
|
|
131 |
|
|
|
132 |
private static function besselj2a(float $ax, int $ord, float $x): float
|
|
|
133 |
{
|
|
|
134 |
$tox = 2.0 / $ax;
|
|
|
135 |
$bjm = self::besselJ0($ax);
|
|
|
136 |
$bj = self::besselJ1($ax);
|
|
|
137 |
for ($j = 1; $j < $ord; ++$j) {
|
|
|
138 |
$bjp = $j * $tox * $bj - $bjm;
|
|
|
139 |
$bjm = $bj;
|
|
|
140 |
$bj = $bjp;
|
|
|
141 |
}
|
|
|
142 |
$ans = $bj;
|
|
|
143 |
|
|
|
144 |
return ($x < 0.0 && ($ord % 2) == 1) ? -$ans : $ans;
|
|
|
145 |
}
|
|
|
146 |
|
|
|
147 |
private static function besselj2b(float $ax, int $ord, float $x): float
|
|
|
148 |
{
|
|
|
149 |
$tox = 2.0 / $ax;
|
|
|
150 |
$jsum = false;
|
|
|
151 |
$bjp = $ans = $sum = 0.0;
|
|
|
152 |
$bj = 1.0;
|
|
|
153 |
for ($j = 2 * ($ord + (int) sqrt(40.0 * $ord)); $j > 0; --$j) {
|
|
|
154 |
$bjm = $j * $tox * $bj - $bjp;
|
|
|
155 |
$bjp = $bj;
|
|
|
156 |
$bj = $bjm;
|
|
|
157 |
if (abs($bj) > 1.0e+10) {
|
|
|
158 |
$bj *= 1.0e-10;
|
|
|
159 |
$bjp *= 1.0e-10;
|
|
|
160 |
$ans *= 1.0e-10;
|
|
|
161 |
$sum *= 1.0e-10;
|
|
|
162 |
}
|
|
|
163 |
if ($jsum === true) {
|
|
|
164 |
$sum += $bj;
|
|
|
165 |
}
|
|
|
166 |
$jsum = $jsum === false;
|
|
|
167 |
if ($j === $ord) {
|
|
|
168 |
$ans = $bjp;
|
|
|
169 |
}
|
|
|
170 |
}
|
|
|
171 |
$sum = 2.0 * $sum - $bj;
|
|
|
172 |
$ans /= $sum;
|
|
|
173 |
|
|
|
174 |
return ($x < 0.0 && ($ord % 2) === 1) ? -$ans : $ans;
|
|
|
175 |
}
|
|
|
176 |
}
|