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