BesselJ.php 5.8 KB

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