UNPKG

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