123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193 |
- /* Bessel J function
- 1st arg x
- 2nd arg n
- Recurrence relation
- besselj(x,n) = (2/x) (n-1) besselj(x,n-1) - besselj(x,n-2)
- besselj(x,1/2) = sqrt(2/pi/x) sin(x)
- besselj(x,-1/2) = sqrt(2/pi/x) cos(x)
- For negative n, reorder the recurrence relation as
- besselj(x,n-2) = (2/x) (n-1) besselj(x,n-1) - besselj(x,n)
- Substitute n+2 for n to obtain
- besselj(x,n) = (2/x) (n+1) besselj(x,n+1) - besselj(x,n+2)
- Examples
- besselj(x,3/2) = (1/x) besselj(x,1/2) - besselj(x,-1/2)
- besselj(x,-3/2) = -(1/x) besselj(x,-1/2) - besselj(x,1/2)
- */
- #include "stdafx.h"
- #include "defs.h"
- #include "math.h"
- void
- eval_besselj(void)
- {
- push(cadr(p1));
- eval();
- push(caddr(p1));
- eval();
- besselj();
- }
- void
- besselj(void)
- {
- save();
- yybesselj();
- restore();
- }
- #define X p1
- #define N p2
- #define SGN p3
- void
- yybesselj(void)
- {
- double d;
- int n;
- N = pop();
- X = pop();
- push(N);
- n = pop_integer();
- // numerical result
- if (isdouble(X) && n != (int) 0x80000000) {
- d = jn(n, X->u.d);
- push_double(d);
- return;
- }
- // bessej(0,0) = 1
- if (iszero(X) && iszero(N)) {
- push_integer(1);
- return;
- }
- // besselj(0,n) = 0
- if (iszero(X) && n != (int) 0x80000000) {
- push_integer(0);
- return;
- }
- // half arguments
- if (N->k == NUM && MEQUAL(N->u.q.b, 2)) {
- // n = 1/2
- if (MEQUAL(N->u.q.a, 1)) {
- push_integer(2);
- push_symbol(PI);
- divide();
- push(X);
- divide();
- push_rational(1, 2);
- power();
- push(X);
- sine();
- multiply();
- return;
- }
- // n = -1/2
- if (MEQUAL(N->u.q.a, -1)) {
- push_integer(2);
- push_symbol(PI);
- divide();
- push(X);
- divide();
- push_rational(1, 2);
- power();
- push(X);
- cosine();
- multiply();
- return;
- }
- // besselj(x,n) = (2/x) (n-sgn(n)) besselj(x,n-sgn(n)) - besselj(x,n-2*sgn(n))
- push_integer(MSIGN(N->u.q.a));
- SGN = pop();
- push_integer(2);
- push(X);
- divide();
- push(N);
- push(SGN);
- subtract();
- multiply();
- push(X);
- push(N);
- push(SGN);
- subtract();
- besselj();
- multiply();
- push(X);
- push(N);
- push_integer(2);
- push(SGN);
- multiply();
- subtract();
- besselj();
- subtract();
- return;
- }
- #if 0 // test cases needed
- if (isnegativeterm(X)) {
- push(X);
- negate();
- push(N);
- power();
- push(X);
- push(N);
- negate();
- power();
- multiply();
- push_symbol(BESSELJ);
- push(X);
- negate();
- push(N);
- list(3);
- multiply();
- return;
- }
- if (isnegativeterm(N)) {
- push_integer(-1);
- push(N);
- power();
- push_symbol(BESSELJ);
- push(X);
- push(N);
- negate();
- list(3);
- multiply();
- return;
- }
- #endif
- push(symbol(BESSELJ));
- push(X);
- push(N);
- list(3);
- }
|