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