Proyectos de Subversion Moodle

Rev

| Ultima modificación | Ver Log |

Rev Autor Línea Nro. Línea
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
}