[Summary view]

Differences Between: [Versions 400 and 401] [Versions 401 and 402] [Versions 401 and 403]

   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($x, $ord)
  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          // special cases
  62          switch ($ord) {
  63              case 0:
  64                  return self::besselJ0($x);
  65              case 1:
  66                  return self::besselJ1($x);
  67          }
  68  
  69          return self::besselJ2($x, $ord);
  70      }
  71  
  72      private static function besselJ0(float $x): float
  73      {
  74          $ax = abs($x);
  75  
  76          if ($ax < 8.0) {
  77              $y = $x * $x;
  78              $ans1 = 57568490574.0 + $y * (-13362590354.0 + $y * (651619640.7 + $y * (-11214424.18 + $y *
  79                              (77392.33017 + $y * (-184.9052456)))));
  80              $ans2 = 57568490411.0 + $y * (1029532985.0 + $y * (9494680.718 + $y * (59272.64853 + $y *
  81                              (267.8532712 + $y * 1.0))));
  82  
  83              return $ans1 / $ans2;
  84          }
  85  
  86          $z = 8.0 / $ax;
  87          $y = $z * $z;
  88          $xx = $ax - 0.785398164;
  89          $ans1 = 1.0 + $y * (-0.1098628627e-2 + $y * (0.2734510407e-4 + $y * (-0.2073370639e-5 + $y * 0.2093887211e-6)));
  90          $ans2 = -0.1562499995e-1 + $y * (0.1430488765e-3 + $y * (-0.6911147651e-5 + $y *
  91                      (0.7621095161e-6 - $y * 0.934935152e-7)));
  92  
  93          return sqrt(0.636619772 / $ax) * (cos($xx) * $ans1 - $z * sin($xx) * $ans2);
  94      }
  95  
  96      private static function besselJ1(float $x): float
  97      {
  98          $ax = abs($x);
  99  
 100          if ($ax < 8.0) {
 101              $y = $x * $x;
 102              $ans1 = $x * (72362614232.0 + $y * (-7895059235.0 + $y * (242396853.1 + $y *
 103                              (-2972611.439 + $y * (15704.48260 + $y * (-30.16036606))))));
 104              $ans2 = 144725228442.0 + $y * (2300535178.0 + $y * (18583304.74 + $y * (99447.43394 + $y *
 105                              (376.9991397 + $y * 1.0))));
 106  
 107              return $ans1 / $ans2;
 108          }
 109  
 110          $z = 8.0 / $ax;
 111          $y = $z * $z;
 112          $xx = $ax - 2.356194491;
 113  
 114          $ans1 = 1.0 + $y * (0.183105e-2 + $y * (-0.3516396496e-4 + $y * (0.2457520174e-5 + $y * (-0.240337019e-6))));
 115          $ans2 = 0.04687499995 + $y * (-0.2002690873e-3 + $y * (0.8449199096e-5 + $y *
 116                      (-0.88228987e-6 + $y * 0.105787412e-6)));
 117          $ans = sqrt(0.636619772 / $ax) * (cos($xx) * $ans1 - $z * sin($xx) * $ans2);
 118  
 119          return ($x < 0.0) ? -$ans : $ans;
 120      }
 121  
 122      private static function besselJ2(float $x, int $ord): float
 123      {
 124          $ax = abs($x);
 125          if ($ax === 0.0) {
 126              return 0.0;
 127          }
 128  
 129          if ($ax > $ord) {
 130              return self::besselj2a($ax, $ord, $x);
 131          }
 132  
 133          return self::besselj2b($ax, $ord, $x);
 134      }
 135  
 136      private static function besselj2a(float $ax, int $ord, float $x)
 137      {
 138          $tox = 2.0 / $ax;
 139          $bjm = self::besselJ0($ax);
 140          $bj = self::besselJ1($ax);
 141          for ($j = 1; $j < $ord; ++$j) {
 142              $bjp = $j * $tox * $bj - $bjm;
 143              $bjm = $bj;
 144              $bj = $bjp;
 145          }
 146          $ans = $bj;
 147  
 148          return ($x < 0.0 && ($ord % 2) == 1) ? -$ans : $ans;
 149      }
 150  
 151      private static function besselj2b(float $ax, int $ord, float $x)
 152      {
 153          $tox = 2.0 / $ax;
 154          $jsum = false;
 155          $bjp = $ans = $sum = 0.0;
 156          $bj = 1.0;
 157          for ($j = 2 * ($ord + (int) sqrt(40.0 * $ord)); $j > 0; --$j) {
 158              $bjm = $j * $tox * $bj - $bjp;
 159              $bjp = $bj;
 160              $bj = $bjm;
 161              if (abs($bj) > 1.0e+10) {
 162                  $bj *= 1.0e-10;
 163                  $bjp *= 1.0e-10;
 164                  $ans *= 1.0e-10;
 165                  $sum *= 1.0e-10;
 166              }
 167              if ($jsum === true) {
 168                  $sum += $bj;
 169              }
 170              $jsum = !$jsum;
 171              if ($j === $ord) {
 172                  $ans = $bjp;
 173              }
 174          }
 175          $sum = 2.0 * $sum - $bj;
 176          $ans /= $sum;
 177  
 178          return ($x < 0.0 && ($ord % 2) === 1) ? -$ans : $ans;
 179      }
 180  }