[Summary view]

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

   1  <?php
   2  
   3  namespace PhpOffice\PhpSpreadsheet\Calculation\Engineering;
   4  
   5  use PhpOffice\PhpSpreadsheet\Calculation\Exception;
   6  use PhpOffice\PhpSpreadsheet\Calculation\Functions;
   7  
   8  class BesselI
   9  {
  10      /**
  11       * BESSELI.
  12       *
  13       *    Returns the modified Bessel function In(x), which is equivalent to the Bessel function evaluated
  14       *        for purely imaginary arguments
  15       *
  16       *    Excel Function:
  17       *        BESSELI(x,ord)
  18       *
  19       * NOTE: The MS Excel implementation of the BESSELI function is still not accurate.
  20       *       This code provides a more accurate calculation
  21       *
  22       * @param mixed $x A float value at which to evaluate the function.
  23       *                                If x is nonnumeric, BESSELI returns the #VALUE! error value.
  24       * @param mixed $ord The integer order of the Bessel function.
  25       *                                If ord is not an integer, it is truncated.
  26       *                                If $ord is nonnumeric, BESSELI returns the #VALUE! error value.
  27       *                                If $ord < 0, BESSELI returns the #NUM! error value.
  28       *
  29       * @return float|string Result, or a string containing an error
  30       */
  31      public static function BESSELI($x, $ord)
  32      {
  33          $x = Functions::flattenSingleValue($x);
  34          $ord = Functions::flattenSingleValue($ord);
  35  
  36          try {
  37              $x = EngineeringValidations::validateFloat($x);
  38              $ord = EngineeringValidations::validateInt($ord);
  39          } catch (Exception $e) {
  40              return $e->getMessage();
  41          }
  42  
  43          if ($ord < 0) {
  44              return Functions::NAN();
  45          }
  46  
  47          $fResult = self::calculate($x, $ord);
  48  
  49          return (is_nan($fResult)) ? Functions::NAN() : $fResult;
  50      }
  51  
  52      private static function calculate(float $x, int $ord): float
  53      {
  54          // special cases
  55          switch ($ord) {
  56              case 0:
  57                  return self::besselI0($x);
  58              case 1:
  59                  return self::besselI1($x);
  60          }
  61  
  62          return self::besselI2($x, $ord);
  63      }
  64  
  65      private static function besselI0(float $x): float
  66      {
  67          $ax = abs($x);
  68  
  69          if ($ax < 3.75) {
  70              $y = $x / 3.75;
  71              $y = $y * $y;
  72  
  73              return 1.0 + $y * (3.5156229 + $y * (3.0899424 + $y * (1.2067492
  74                                  + $y * (0.2659732 + $y * (0.360768e-1 + $y * 0.45813e-2)))));
  75          }
  76  
  77          $y = 3.75 / $ax;
  78  
  79          return (exp($ax) / sqrt($ax)) * (0.39894228 + $y * (0.1328592e-1 + $y * (0.225319e-2 + $y * (-0.157565e-2
  80                              + $y * (0.916281e-2 + $y * (-0.2057706e-1 + $y * (0.2635537e-1 +
  81                                          $y * (-0.1647633e-1 + $y * 0.392377e-2))))))));
  82      }
  83  
  84      private static function besselI1(float $x): float
  85      {
  86          $ax = abs($x);
  87  
  88          if ($ax < 3.75) {
  89              $y = $x / 3.75;
  90              $y = $y * $y;
  91              $ans = $ax * (0.5 + $y * (0.87890594 + $y * (0.51498869 + $y * (0.15084934 + $y * (0.2658733e-1 +
  92                                      $y * (0.301532e-2 + $y * 0.32411e-3))))));
  93  
  94              return ($x < 0.0) ? -$ans : $ans;
  95          }
  96  
  97          $y = 3.75 / $ax;
  98          $ans = 0.2282967e-1 + $y * (-0.2895312e-1 + $y * (0.1787654e-1 - $y * 0.420059e-2));
  99          $ans = 0.39894228 + $y * (-0.3988024e-1 + $y * (-0.362018e-2 + $y * (0.163801e-2 +
 100                          $y * (-0.1031555e-1 + $y * $ans))));
 101          $ans *= exp($ax) / sqrt($ax);
 102  
 103          return ($x < 0.0) ? -$ans : $ans;
 104      }
 105  
 106      private static function besselI2(float $x, int $ord): float
 107      {
 108          if ($x === 0.0) {
 109              return 0.0;
 110          }
 111  
 112          $tox = 2.0 / abs($x);
 113          $bip = 0;
 114          $ans = 0.0;
 115          $bi = 1.0;
 116  
 117          for ($j = 2 * ($ord + (int) sqrt(40.0 * $ord)); $j > 0; --$j) {
 118              $bim = $bip + $j * $tox * $bi;
 119              $bip = $bi;
 120              $bi = $bim;
 121  
 122              if (abs($bi) > 1.0e+12) {
 123                  $ans *= 1.0e-12;
 124                  $bi *= 1.0e-12;
 125                  $bip *= 1.0e-12;
 126              }
 127  
 128              if ($j === $ord) {
 129                  $ans = $bip;
 130              }
 131          }
 132  
 133          $ans *= self::besselI0($x) / $bi;
 134  
 135          return ($x < 0.0 && (($ord % 2) === 1)) ? -$ans : $ans;
 136      }
 137  }