besselj.cpp

   1 /* Bessel J function
   2
   3         1st arg         x
   4
   5         2nd arg         n
   6
   7 Recurrence relation
   8
   9         besselj(x,n) = (2/x) (n-1) besselj(x,n-1) - besselj(x,n-2)
  10
  11         besselj(x,1/2) = sqrt(2/pi/x) sin(x)
  12
  13         besselj(x,-1/2) = sqrt(2/pi/x) cos(x)
  14
  15 For negative n, reorder the recurrence relation as
  16
  17         besselj(x,n-2) = (2/x) (n-1) besselj(x,n-1) - besselj(x,n)
  18
  19 Substitute n+2 for n to obtain
  20
  21         besselj(x,n) = (2/x) (n+1) besselj(x,n+1) - besselj(x,n+2)
  22
  23 Examples
  24
  25         besselj(x,3/2) = (1/x) besselj(x,1/2) - besselj(x,-1/2)
  26
  27         besselj(x,-3/2) = -(1/x) besselj(x,-1/2) - besselj(x,1/2)
  28 */
  29
  30 #include "stdafx.h"
  31 #include "defs.h"
  32 #include "math.h"
  33 void
  34 eval_besselj(void)
  35 {
  36         push(cadr(p1));
  37         eval();
  38         push(caddr(p1));
  39         eval();
  40         besselj();
  41 }
  42
  43 void
  44 besselj(void)
  45 {
  46         save();
  47         yybesselj();
  48         restore();
  49 }
  50
  51 #define X p1
  52 #define N p2
  53 #define SGN p3
  54
  55 void
  56 yybesselj(void)
  57 {
  58         double d;
  59         int n;
  60
  61         N = pop();
  62         X = pop();
  63
  64         push(N);
  65         n = pop_integer();
  66
  67         // numerical result
  68
  69         if (isdouble(X) && n != (int) 0x80000000) {
  70                 d = jn(n, X->u.d);
  71                 push_double(d);
  72                 return;
  73         }
  74
  75         // bessej(0,0) = 1
  76
  77         if (iszero(X) && iszero(N)) {
  78                 push_integer(1);
  79                 return;
  80         }
  81
  82         // besselj(0,n) = 0
  83
  84         if (iszero(X) && n != (int) 0x80000000) {
  85                 push_integer(0);
  86                 return;
  87         }
  88
  89         // half arguments
  90
  91         if (N->k == NUM && MEQUAL(N->u.q.b, 2)) {
  92
  93                 // n = 1/2
  94
  95                 if (MEQUAL(N->u.q.a, 1)) {
  96                         push_integer(2);
  97                         push_symbol(PI);
  98                         divide();
  99                         push(X);
 100                         divide();
 101                         push_rational(1, 2);
 102                         power();
 103                         push(X);
 104                         sine();
 105                         multiply();
 106                         return;
 107                 }
 108
 109                 // n = -1/2
 110
 111                 if (MEQUAL(N->u.q.a, -1)) {
 112                         push_integer(2);
 113                         push_symbol(PI);
 114                         divide();
 115                         push(X);
 116                         divide();
 117                         push_rational(1, 2);
 118                         power();
 119                         push(X);
 120                         cosine();
 121                         multiply();
 122                         return;
 123                 }
 124
 125                 // besselj(x,n) = (2/x) (n-sgn(n)) besselj(x,n-sgn(n)) - besselj(x,n-2*sgn(n))
 126
 127                 push_integer(MSIGN(N->u.q.a));
 128                 SGN = pop();
 129
 130                 push_integer(2);
 131                 push(X);
 132                 divide();
 133                 push(N);
 134                 push(SGN);
 135                 subtract();
 136                 multiply();
 137                 push(X);
 138                 push(N);
 139                 push(SGN);
 140                 subtract();
 141                 besselj();
 142                 multiply();
 143                 push(X);
 144                 push(N);
 145                 push_integer(2);
 146                 push(SGN);
 147                 multiply();
 148                 subtract();
 149                 besselj();
 150                 subtract();
 151
 152                 return;
 153         }
 154
 155 #if 0 // test cases needed
 156         if (isnegativeterm(X)) {
 157                 push(X);
 158                 negate();
 159                 push(N);
 160                 power();
 161                 push(X);
 162                 push(N);
 163                 negate();
 164                 power();
 165                 multiply();
 166                 push_symbol(BESSELJ);
 167                 push(X);
 168                 negate();
 169                 push(N);
 170                 list(3);
 171                 multiply();
 172                 return;
 173         }
 174
 175         if (isnegativeterm(N)) {
 176                 push_integer(-1);
 177                 push(N);
 178                 power();
 179                 push_symbol(BESSELJ);
 180                 push(X);
 181                 push(N);
 182                 negate();
 183                 list(3);
 184                 multiply();
 185                 return;
 186         }
 187 #endif
 188
 189         push(symbol(BESSELJ));
 190         push(X);
 191         push(N);
 192         list(3);
 193 }