tablepress / libraries / vendor / PhpSpreadsheet / Calculation / Engineering / BesselI.php

BesselI.php in TablePress – Tables in WordPress made easy 3.4, at libraries/vendor/PhpSpreadsheet/Calculation/Engineering/BesselI.php

145 lines 3.9 KB

Raw Download Zip

1	<?php
2
3	namespace TablePress\PhpOffice\PhpSpreadsheet\Calculation\Engineering;
4
5	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled;
6	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Exception;
7	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
8
9	class BesselI
10	{
11	use ArrayEnabled;
12
13	/**
14	* BESSELI.
15	*
16	* Returns the modified Bessel function In(x), which is equivalent to the Bessel function evaluated
17	* for purely imaginary arguments
18	*
19	* Excel Function:
20	* BESSELI(x,ord)
21	*
22	* NOTE: The MS Excel implementation of the BESSELI function is still not accurate.
23	* This code provides a more accurate calculation
24	*
25	* @param mixed $x A float value at which to evaluate the function.
26	* If x is nonnumeric, BESSELI returns the #VALUE! error value.
27	* Or can be an array of values
28	* @param mixed $ord The integer order of the Bessel function.
29	* If ord is not an integer, it is truncated.
30	* If $ord is nonnumeric, BESSELI returns the #VALUE! error value.
31	* If $ord < 0, BESSELI returns the #NUM! error value.
32	* Or can be an array of values
33	*
34	* @return array<mixed>\|float\|string Result, or a string containing an error
35	* If an array of numbers is passed as an argument, then the returned result will also be an array
36	* with the same dimensions
37	*/
38	public static function BESSELI($x, $ord)
39	{
40	if (is_array($x) \|\| is_array($ord)) {
41	return self::evaluateArrayArguments([self::class, __FUNCTION__], $x, $ord);
42	}
43
44	try {
45	$x = EngineeringValidations::validateFloat($x);
46	$ord = EngineeringValidations::validateInt($ord);
47	} catch (Exception $e) {
48	return $e->getMessage();
49	}
50
51	if ($ord < 0) {
52	return ExcelError::NAN();
53	}
54
55	$fResult = self::calculate($x, $ord);
56
57	return (is_nan($fResult)) ? ExcelError::NAN() : $fResult;
58	}
59
60	private static function calculate(float $x, int $ord): float
61	{
62	switch ($ord) {
63	case 0:
64	return self::besselI0($x);
65	case 1:
66	return self::besselI1($x);
67	default:
68	return self::besselI2($x, $ord);
69	}
70	}
71
72	private static function besselI0(float $x): float
73	{
74	$ax = abs($x);
75
76	if ($ax < 3.75) {
77	$y = $x / 3.75;
78	$y = $y * $y;
79
80	return 1.0 + $y * (3.5156229 + $y * (3.0899424 + $y * (1.2067492
81	+ $y * (0.2659732 + $y * (0.360768e-1 + $y * 0.45813e-2)))));
82	}
83
84	$y = 3.75 / $ax;
85
86	return (exp($ax) / sqrt($ax)) * (0.39894228 + $y * (0.1328592e-1 + $y * (0.225319e-2 + $y * (-0.157565e-2
87	+ $y * (0.916281e-2 + $y * (-0.2057706e-1 + $y * (0.2635537e-1
88	+ $y * (-0.1647633e-1 + $y * 0.392377e-2))))))));
89	}
90
91	private static function besselI1(float $x): float
92	{
93	$ax = abs($x);
94
95	if ($ax < 3.75) {
96	$y = $x / 3.75;
97	$y = $y * $y;
98	$ans = $ax * (0.5 + $y * (0.87890594 + $y * (0.51498869 + $y * (0.15084934 + $y * (0.2658733e-1
99	+ $y * (0.301532e-2 + $y * 0.32411e-3))))));
100
101	return ($x < 0.0) ? -$ans : $ans;
102	}
103
104	$y = 3.75 / $ax;
105	$ans = 0.2282967e-1 + $y * (-0.2895312e-1 + $y * (0.1787654e-1 - $y * 0.420059e-2));
106	$ans = 0.39894228 + $y * (-0.3988024e-1 + $y * (-0.362018e-2 + $y * (0.163801e-2
107	+ $y * (-0.1031555e-1 + $y * $ans))));
108	$ans *= exp($ax) / sqrt($ax);
109
110	return ($x < 0.0) ? -$ans : $ans;
111	}
112
113	private static function besselI2(float $x, int $ord): float
114	{
115	if ($x === 0.0) {
116	return 0.0;
117	}
118
119	$tox = 2.0 / abs($x);
120	$bip = 0;
121	$ans = 0.0;
122	$bi = 1.0;
123
124	for ($j = 2 * ($ord + (int) sqrt(40.0 * $ord)); $j > 0; --$j) {
125	$bim = $bip + $j * $tox * $bi;
126	$bip = $bi;
127	$bi = $bim;
128
129	if (abs($bi) > 1.0e+12) {
130	$ans *= 1.0e-12;
131	$bi *= 1.0e-12;
132	$bip *= 1.0e-12;
133	}
134
135	if ($j === $ord) {
136	$ans = $bip;
137	}
138	}
139
140	$ans *= self::besselI0($x) / $bi;
141
142	return ($x < 0.0 && (($ord % 2) === 1)) ? -$ans : $ans;
143	}
144	}
145