?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Replace Python program auto-generated for documentation categories
[maxima/cygwin.git] / src / bessel.lisp
blob45e9439aa632a6b92a34d5255b697eccb42447e0
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
4 ;;; ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8
9 (in-package :maxima)
10
11 (declare-top (special $hypergeometric_representation))
12
13 ;; When non-NIL, the Bessel functions of half-integral order are
14 ;; expanded in terms of elementary functions.
15
16 (defmvar $besselexpand nil)
17
18 ;; When T Bessel functions with an integer order are reduced to order 0 and 1
19 (defmvar $bessel_reduce nil)
20
21 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
22
23 ;;; Helper functions for this file
24
25 ;; If E is a maxima ratio with a denominator of DEN, return the ratio
26 ;; as a Lisp rational. Otherwise NIL.
27 (defun max-numeric-ratio-p (e den)
28 (if (and (listp e)
29 (eq 'rat (caar e))
30 (= den (third e))
31 (integerp (second e)))
32 (/ (second e) (third e))
33 nil))
34
35 (defun bessel-numerical-eval-p (order arg)
36 ;; Return non-NIL if we should numerically evaluate a bessel
37 ;; function. Basically, both args have to be numbers. If both args
38 ;; are integers, we don't evaluate unless $numer is true.
39 (or (and (numberp order) (complex-number-p arg)
40 (or (floatp order)
41 (floatp ($realpart arg))
42 (floatp ($imagpart arg))))
43 (and $numer (numberp order)
44 (complex-number-p arg))))
45
46 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
47 ;;;
48 ;;; Implementation of the Bessel J function
49 ;;;
50 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
51
52 (defmfun $bessel_j (v z)
53 (simplify (list '(%bessel_j) v z)))
54
55 #||
56 (defprop $bessel_j %bessel_j alias)
57 (defprop $bessel_j %bessel_j verb)
58 (defprop %bessel_j $bessel_j reversealias)
59 (defprop %bessel_j $bessel_j noun)
60
61 ;; Bessel J is a simplifying function.
62
63 (defprop %bessel_j simp-bessel-j operators)
64 ||#
65 ;; Bessel J distributes over lists, matrices, and equations
66
67 (defprop %bessel_j (mlist $matrix mequal) distribute_over)
68
69 ;; Derivatives of the Bessel function.
70 (defprop %bessel_j
71 ((n x)
72 ;; Derivative wrt to order n. A&S 9.1.64. Do we really want to
73 ;; do this? It's quite messy.
74 ;;
75 ;; J[n](x)*log(x/2)
76 ;; - (x/2)^n*sum((-1)^k*psi(n+k+1)/gamma(n+k+1)*(x^2/4)^k/k!,k,0,inf)
77 ((mplus)
78 ((mtimes)
79 ((%bessel_j) n x)
80 ((%log) ((mtimes) ((rat) 1 2) x)))
81 ((mtimes) -1
82 ((mexpt) ((mtimes) x ((rat) 1 2)) n)
83 ((%sum)
84 ((mtimes) ((mexpt) -1 $%k)
85 ((mexpt) ((mfactorial) $%k) -1)
86 ((mqapply) (($psi array) 0) ((mplus) 1 $%k n))
87 ((mexpt) ((%gamma) ((mplus) 1 $%k n)) -1)
88 ((mexpt) ((mtimes) x x ((rat) 1 4)) $%k))
89 $%k 0 $inf)))
90
91 ;; Derivative wrt to arg x.
92 ;; A&S 9.1.27; changed from 9.1.30 so that taylor works on Bessel functions
93 ((mtimes)
94 ((mplus)
95 ((%bessel_j) ((mplus) -1 n) x)
96 ((mtimes) -1 ((%bessel_j) ((mplus) 1 n) x)))
97 ((rat) 1 2)))
98 grad)
99
100 ;; Integral of the Bessel function wrt z
101 (defun bessel-j-integral-2 (v z)
102 (case v
103 (0
104 ;; integrate(bessel_j(0,z),z)
105 ;; = (1/2)*z*(%pi*bessel_j(1,z)*struve_h(0,z)
106 ;; +bessel_j(0,z)*(2-%pi*struve_h(1,z)))
107 `((mtimes) ((rat) 1 2) ,z
108 ((mplus)
109 ((mtimes) $%pi
110 ((%bessel_j) 1 ,z)
111 ((%struve_h) 0 ,z))
112 ((mtimes)
113 ((%bessel_j) 0 ,z)
114 ((mplus) 2 ((mtimes) -1 $%pi ((%struve_h) 1 ,z)))))))
115 (1
116 ;; integrate(bessel_j(1,z),z) = -bessel_j(0,z)
117 `((mtimes) -1 ((%bessel_j) 0 ,z)))
118 (otherwise
119 ;; http://functions.wolfram.com/03.01.21.0002.01
120 ;; integrate(bessel_j(v,z),z)
121 ;; = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2)
122 ;; * hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],-z^2/4)
123 ;; = 2^(-v)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],-z^2/4)
124 ;; / gamma(v+2)
125 `((mtimes)
126 (($hypergeometric)
127 ((mlist)
128 ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) ,v)))
129 ((mlist)
130 ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) ,v))
131 ((mplus) 1 ,v))
132 ((mtimes) ((rat) -1 4) ((mexpt) ,z 2)))
133 ((mexpt) 2 ((mtimes) -1 ,v))
134 ((mexpt) ((%gamma) ((mplus) 2 ,v)) -1)
135 ((mexpt) z ((mplus) 1 ,v))))))
136
137 (putprop '%bessel_j `((v z) nil ,#'bessel-j-integral-2) 'integral)
138
139 ;; Support a simplim%bessel_j function to handle specific limits
140
141 (defprop %bessel_j simplim%bessel_j simplim%function)
142
143 (defun simplim%bessel_j (expr var val)
144 ;; Look for the limit of the arguments.
145 (let ((v (limit (cadr expr) var val 'think))
146 (z (limit (caddr expr) var val 'think)))
147 (cond
148 ;; Handle an argument 0 at this place.
149 ((or (zerop1 z)
150 (eq z '$zeroa)
151 (eq z '$zerob))
152 (let ((sgn ($sign ($realpart v))))
153 (cond ((and (eq sgn '$neg)
154 (not (maxima-integerp v)))
155 ;; bessel_j(v,0), Re(v)<0 and v not an integer
156 '$infinity)
157 ((and (eq sgn '$zero)
158 (not (zerop1 v)))
159 ;; bessel_j(v,0), Re(v)=0 and v #0
160 '$und)
161 ;; Call the simplifier of the function.
162 (t (simplify (list '(%bessel_j) v z))))))
163 ((or (eq z '$inf)
164 (eq z '$minf))
165 ;; bessel_j(v,inf) or bessel_j(v,minf)
166 0)
167 (t
168 ;; All other cases are handled by the simplifier of the function.
169 (simplify (list '(%bessel_j) v z))))))
170
171 #+nil
172 (defun simp-bessel-j (expr ignored z)
173 (declare (ignore ignored))
174 (twoargcheck expr)
175 (let ((order (simpcheck (cadr expr) z))
176 (arg (simpcheck (caddr expr) z))
177 (rat-order nil))
178 (cond
179 ((zerop1 arg)
180 ;; We handle the different case for zero arg carefully.
181 (let ((sgn ($sign ($realpart order))))
182 (cond ((and (eq sgn '$zero)
183 (zerop1 ($imagpart order)))
184 ;; bessel_j(0,0) = 1
185 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1))
186 ((or (floatp order) (floatp arg)) 1.0)
187 (t 1)))
188 ((or (eq sgn '$pos)
189 (maxima-integerp order))
190 ;; bessel_j(v,0) and Re(v)>0 or v an integer
191 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0))
192 ((or (floatp order) (floatp arg)) 0.0)
193 (t 0)))
194 ((and (eq sgn '$neg)
195 (not (maxima-integerp order)))
196 ;; bessel_j(v,0) and Re(v)<0 and v not an integer
197 (simp-domain-error
198 (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.")
199 order arg))
200 ((and (eq sgn '$zero)
201 (not (zerop1 ($imagpart order))))
202 ;; bessel_j(v,0) and Re(v)=0 and v # 0
203 (simp-domain-error
204 (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.")
205 order arg))
206 (t
207 ;; No information about the sign of the order
208 (eqtest (list '(%bessel_j) order arg) expr)))))
209
210 ((complex-float-numerical-eval-p order arg)
211 ;; We have numeric order and arg and $numer is true, or we have either
212 ;; the order or arg being floating-point, so let's evaluate it
213 ;; numerically.
214 ;; The numerical routine bessel-j returns a CL number, so we have
215 ;; to add the conversion to a Maxima-complex-number.
216 (cond ((= 0 ($imagpart order))
217 ;; order is real, arg is real or complex
218 (complexify (bessel-j ($float order)
219 (complex ($float ($realpart arg))
220 ($float ($imagpart arg))))))
221 (t
222 ;; order is complex, arg is real or complex
223 (let (($numer t)
224 ($float t)
225 (order ($float order))
226 (arg ($float arg)))
227 ($float
228 ($rectform
229 (bessel-j-hypergeometric order arg)))))))
230
231 ((and (integerp order) (minusp order))
232 ;; Some special cases when the order is an integer.
233 ;; A&S 9.1.5: J[-n](x) = (-1)^n*J[n](x)
234 (if (evenp order)
235 (take '(%bessel_j) (- order) arg)
236 (mul -1 (take '(%bessel_j) (- order) arg))))
237
238 ((and $besselexpand
239 (setq rat-order (max-numeric-ratio-p order 2)))
240 ;; When order is a fraction with a denominator of 2, we
241 ;; can express the result in terms of elementary functions.
242 (bessel-j-half-order rat-order arg))
243
244 ((and $bessel_reduce
245 (and (integerp order)
246 (plusp order)
247 (> order 1)))
248 ;; Reduce a bessel function of order > 2 to order 1 and 0.
249 ;; A&S 9.1.27: bessel_j(v,z) = 2*(v-1)/z*bessel_j(v-1,z)-bessel_j(v-2,z)
250 (sub (mul 2
251 (- order 1)
252 (inv arg)
253 (take '(%bessel_j) (- order 1) arg))
254 (take '(%bessel_j) (- order 2) arg)))
255
256 ((and $%iargs (multiplep arg '$%i))
257 ;; bessel_j(v, %i*x) = (%i*x)^v/(x^v) * bessel_i(v, x)
258 ;; (From http://functions.wolfram.com/03.01.27.0002.01)
259 (let ((x (coeff arg '$%i 1)))
260 (mul (power (mul '$%i x) order)
261 (inv (power x order))
262 (take '(%bessel_i) order x))))
263
264 ($hypergeometric_representation
265 (bessel-j-hypergeometric order arg))
266
267 (t
268 (eqtest (list '(%bessel_j) order arg) expr)))))
269
270 (def-simplifier bessel_j (order arg)
271 (let ((rat-order nil))
272 (cond
273 ((zerop1 arg)
274 ;; We handle the different case for zero arg carefully.
275 (let ((sgn ($sign ($realpart order))))
276 (cond ((and (eq sgn '$zero)
277 (zerop1 ($imagpart order)))
278 ;; bessel_j(0,0) = 1
279 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1))
280 ((or (floatp order) (floatp arg)) 1.0)
281 (t 1)))
282 ((or (eq sgn '$pos)
283 (maxima-integerp order))
284 ;; bessel_j(v,0) and Re(v)>0 or v an integer
285 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0))
286 ((or (floatp order) (floatp arg)) 0.0)
287 (t 0)))
288 ((and (eq sgn '$neg)
289 (not (maxima-integerp order)))
290 ;; bessel_j(v,0) and Re(v)<0 and v not an integer
291 (simp-domain-error
292 (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.")
293 order arg))
294 ((and (eq sgn '$zero)
295 (not (zerop1 ($imagpart order))))
296 ;; bessel_j(v,0) and Re(v)=0 and v # 0
297 (simp-domain-error
298 (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.")
299 order arg))
300 (t
301 ;; No information about the sign of the order
302 (give-up)))))
303
304 ((complex-float-numerical-eval-p order arg)
305 ;; We have numeric order and arg and $numer is true, or we have either
306 ;; the order or arg being floating-point, so let's evaluate it
307 ;; numerically.
308 ;; The numerical routine bessel-j returns a CL number, so we have
309 ;; to add the conversion to a Maxima-complex-number.
310 (cond ((= 0 ($imagpart order))
311 ;; order is real, arg is real or complex
312 (complexify (bessel-j ($float order)
313 (complex ($float ($realpart arg))
314 ($float ($imagpart arg))))))
315 (t
316 ;; order is complex, arg is real or complex
317 (let (($numer t)
318 ($float t)
319 (order ($float order))
320 (arg ($float arg)))
321 ($float
322 ($rectform
323 (bessel-j-hypergeometric order arg)))))))
324
325 ((and (integerp order) (minusp order))
326 ;; Some special cases when the order is an integer.
327 ;; A&S 9.1.5: J[-n](x) = (-1)^n*J[n](x)
328 (if (evenp order)
329 (take '(%bessel_j) (- order) arg)
330 (mul -1 (take '(%bessel_j) (- order) arg))))
331
332 ((and $besselexpand
333 (setq rat-order (max-numeric-ratio-p order 2)))
334 ;; When order is a fraction with a denominator of 2, we
335 ;; can express the result in terms of elementary functions.
336 (bessel-j-half-order rat-order arg))
337
338 ((and $bessel_reduce
339 (and (integerp order)
340 (plusp order)
341 (> order 1)))
342 ;; Reduce a bessel function of order > 2 to order 1 and 0.
343 ;; A&S 9.1.27: bessel_j(v,z) = 2*(v-1)/z*bessel_j(v-1,z)-bessel_j(v-2,z)
344 (sub (mul 2
345 (- order 1)
346 (inv arg)
347 (take '(%bessel_j) (- order 1) arg))
348 (take '(%bessel_j) (- order 2) arg)))
349
350 ((and $%iargs (multiplep arg '$%i))
351 ;; bessel_j(v, %i*x) = (%i*x)^v/(x^v) * bessel_i(v, x)
352 ;; (From http://functions.wolfram.com/03.01.27.0002.01)
353 (let ((x (coeff arg '$%i 1)))
354 (mul (power (mul '$%i x) order)
355 (inv (power x order))
356 (take '(%bessel_i) order x))))
357
358 ($hypergeometric_representation
359 (bessel-j-hypergeometric order arg))
360
361 (t
362 (give-up)))))
363
364 ;; Returns the hypergeometric representation of bessel_j
365 (defun bessel-j-hypergeometric (order arg)
366 ;; A&S 9.1.69 / http://functions.wolfram.com/03.01.26.0002.01
367 ;; hypergeometric([],[v+1],-z^2/4)*(z/2)^v/gamma(v+1)
368 (mul (inv (take '(%gamma) (add order 1)))
369 (power (div arg 2) order)
370 (take '($hypergeometric)
371 (list '(mlist))
372 (list '(mlist) (add order 1))
373 (neg (div (mul arg arg) 4)))))
374
375 ;; Compute value of Bessel function of the first kind of order ORDER.
376 (defun bessel-j (order arg)
377 (cond
378 ((zerop (imagpart arg))
379 ;; We have numeric args and the arg is purely real.
380 ;; Call the real-valued Bessel function when possible.
381 (let ((arg (realpart arg)))
382 (cond
383 ((= order 0)
384 (slatec:dbesj0 (float arg)))
385 ((= order 1)
386 (slatec:dbesj1 (float arg)))
387 ((minusp order)
388 (cond ((zerop (nth-value 1 (truncate order)))
389 ;; The order is a negative integer. We use A&S 9.1.5:
390 ;; J[-n](z) = (-1)^n*J[n](z)
391 ;; and not the Hankel functions.
392 (if (evenp (floor order))
393 (bessel-j (- order) arg)
394 (- (bessel-j (- order) arg))))
395 (t
396 ;; Bessel function of negative order. We use the Hankel
397 ;; functions to compute this:
398 ;; J[v](z) = (H1[v](z) + H2[v](z)) / 2.
399 ;; This works for negative and positive arg and handles
400 ;; special cases correctly.
401 (let ((result (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg)))))
402 (cond ((= (nth-value 1 (floor order)) 1/2)
403 ;; ORDER is a half-integral-value or a float
404 ;; representation, thereof.
405 (if (minusp arg)
406 ;; arg is negative, the result is purely imaginary
407 (complex 0 (imagpart result))
408 ;; arg is positive, the result is purely real
409 (realpart result)))
410 ;; in all other cases general complex result
411 (t result))))))
412 (t
413 ;; We have a real arg and order > 0 and order not 0 or 1
414 ;; for this case we can call the function dbesj
415 (multiple-value-bind (n alpha) (floor (float order))
416 (let ((jvals (make-array (1+ n) :element-type 'flonum)))
417 (slatec:dbesj (abs (float arg)) alpha (1+ n) jvals 0)
418 (cond ((>= arg 0)
419 (aref jvals n))
420 (t
421 ;; Use analytic continuation formula A&S 9.1.35:
422 ;; J[v](z*exp(m*%pi*%i)) = exp(m*%pi*%i*v)*J[v](z)
423 ;; for an integer m. In particular, for m = 1:
424 ;; J[v](-x) = exp(v*%pi*%i)*J[v](x)
425 ;; and handle special cases
426 (cond ((zerop (nth-value 1 (truncate order)))
427 ;; order is an integer
428 (if (evenp (floor order))
429 (aref jvals n)
430 (- (aref jvals n))))
431 ((= (nth-value 1 (floor order)) 1/2)
432 ;; Order is a half-integral-value and we know that
433 ;; arg < 0, so the result is purely imginary.
434 (if (evenp (floor order))
435 (complex 0 (aref jvals n))
436 (complex 0 (- (aref jvals n)))))
437 ;; In all other cases a general complex result
438 (t
439 (* (cis (* order pi))
440 (aref jvals n))))))))))))
441 (t
442 ;; The arg is complex. Use the complex-valued Bessel function.
443 (cond ((mminusp order)
444 ;; Bessel function of negative order. We use the Hankel function to
445 ;; compute this, because A&S 9.1.3 and 9.1.4 say
446 ;; H1[v](z) = J[v](z) + %i * Y[v](z)
447 ;; and
448 ;; H2[v](z) = J[v](z) - %i * Y[v](z).
449 ;; Thus
450 ;; J[v](z) = (H1[v](z) + H2[v](z)) / 2.
451 ;; Not the most efficient way, but perhaps good enough for maxima.
452 (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg))))
453 (t
454 (multiple-value-bind (n alpha) (floor (float order))
455 (let ((cyr (make-array (1+ n) :element-type 'flonum))
456 (cyi (make-array (1+ n) :element-type 'flonum)))
457 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
458 v-cyr v-cyi v-nz v-ierr)
459 (slatec:zbesj (float (realpart arg))
460 (float (imagpart arg))
461 alpha 1 (1+ n) cyr cyi 0 0)
462 (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz))
463
464 ;; Should check the return status in v-ierr of this routine.
465 (when (plusp v-ierr)
466 (format t "zbesj ierr = ~A~%" v-ierr))
467 (complex (aref cyr n) (aref cyi n))))))))))
468
469 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
470 ;;;
471 ;;; Implementation of the Bessel Y function
472 ;;;
473 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
474
475 (defmfun $bessel_y (v z)
476 (simplify (list '(%bessel_y) v z)))
477
478 #||
479 (defprop $bessel_y %bessel_y alias)
480 (defprop $bessel_y %bessel_y verb)
481 (defprop %bessel_y $bessel_y reversealias)
482 (defprop %bessel_y $bessel_y noun)
483
484 (defprop %bessel_y simp-bessel-y operators)
485 ||#
486
487 ;; Bessel Y distributes over lists, matrices, and equations
488
489 (defprop %bessel_y (mlist $matrix mequal) distribute_over)
490
491 (defprop %bessel_y
492 ((n x)
493 ;; Derivative wrt order n. A&S 9.1.65.
494 ;;
495 ;; cot(n*%pi)*[diff(bessel_j(n,x),n)-%pi*bessel_y(n,x)]
496 ;; - csc(n*%pi)*diff(bessel_j(-n,x),n)-%pi*bessel_j(n,x)
497 ((mplus)
498 ((mtimes) -1 $%pi ((%bessel_j) n x))
499 ((mtimes)
500 -1
501 ((%csc) ((mtimes) $%pi n))
502 ((%derivative) ((%bessel_j) ((mtimes) -1 n) x) n 1))
503 ((mtimes)
504 ((%cot) ((mtimes) $%pi n))
505 ((mplus)
506 ((mtimes) -1 $%pi ((%bessel_y) n x))
507 ((%derivative) ((%bessel_j) n x) n 1))))
508
509 ;; Derivative wrt to arg x. A&S 9.1.27; changed from A&S 9.1.30
510 ;; to be consistent with bessel_j.
511 ((mtimes)
512 ((mplus)
513 ((%bessel_y)((mplus) -1 n) x)
514 ((mtimes) -1 ((%bessel_y) ((mplus) 1 n) x)))
515 ((rat) 1 2)))
516 grad)
517
518 ;; Integral of the Bessel Y function wrt z
519 ;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/21/01/01/
520 (defun bessel-y-integral-2 (n z)
521 (cond
522 ((and ($integerp n) (<= 0 n))
523 (cond
524 (($oddp n)
525 ;; integrate(bessel_y(2*N+1,z),z) , N > 0
526 ;; = -bessel_y(0,z) - 2 * sum(bessel_y(2*k,z),k,1,N)
527 (let* ((k (gensym))
528 (answer `((mplus) ((mtimes) -1 ((%bessel_y) 0 ,z))
529 ((mtimes) -2
530 ((%sum) ((%bessel_y) ((mtimes) 2 ,k) ,z) ,k 1
531 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))))
532 ;; Expand out the sum if n < 10. Otherwise fix up the indices
533 (if (< n 10)
534 (meval `(($ev) ,answer $sum)) ; Is there a better way?
535 (simplify ($niceindices answer)))))
536 (($evenp n)
537 ;; integrate(bessel_y(2*N,z),z) , N > 0
538 ;; = (1/2)*%pi*z*(bessel_y(0,z)*struve_h(-1,z)
539 ;; +bessel_y(1,z)*struve_h(0,z))
540 ;; - 2 * sum(bessel_y(2*k+1,z),k,0,N-1)
541 (let*
542 ((k (gensym))
543 (answer `((mplus)
544 ((mtimes) -2
545 ((%sum)
546 ((%bessel_y) ((mplus) 1 ((mtimes) 2 ,k)) ,z)
547 ,k 0
548 ((mplus)
549 -1
550 ((mtimes) ((rat) 1 2) ,n))))
551 ((mtimes) ((rat) 1 2) $%pi ,z
552 ((mplus)
553 ((mtimes)
554 ((%bessel_y) 0 ,z)
555 ((%struve_h) -1 ,z))
556 ((mtimes)
557 ((%bessel_y) 1 ,z)
558 ((%struve_h) 0 ,z)))))))
559 ;; Expand out the sum if n < 10. Otherwise fix up the indices
560 (if (< n 10)
561 (meval `(($ev) ,answer $sum)) ; Is there a better way?
562 (simplify ($niceindices answer)))))))
563 (t nil)))
564
565 (putprop '%bessel_y `((n z) nil ,#'bessel-y-integral-2) 'integral)
566
567 ;; Support a simplim%bessel_y function to handle specific limits
568
569 (defprop %bessel_y simplim%bessel_y simplim%function)
570
571 (defun simplim%bessel_y (expr var val)
572 ;; Look for the limit of the arguments.
573 (let ((v (limit (cadr expr) var val 'think))
574 (z (limit (caddr expr) var val 'think)))
575 (cond
576 ;; Handle an argument 0 at this place.
577 ((or (zerop1 z)
578 (eq z '$zeroa)
579 (eq z '$zerob))
580 (cond ((zerop1 v)
581 ;; bessel_y(0,0)
582 '$minf)
583 ((integerp v)
584 ;; bessel_y(n,0), n an integer
585 (cond ((evenp v) '$minf)
586 (t (cond ((eq z '$zeroa) '$minf)
587 ((eq z '$zerob) '$inf)
588 (t '$infinity)))))
589 ((not (zerop1 ($realpart v)))
590 ;; bessel_y(v,0), Re(v)#0
591 '$infinity)
592 ((and (zerop1 ($realpart v))
593 (not (zerop1 v)))
594 ;; bessel_y(v,0), Re(v)=0 and v#0
595 '$und)
596 ;; Call the simplifier of the function.
597 (t (simplify (list '(%bessel_y) v z)))))
598 ((or (eq z '$inf)
599 (eq z '$minf))
600 ;; bessel_y(v,inf) or bessel_y(v,minf)
601 0)
602 (t
603 ;; All other cases are handled by the simplifier of the function.
604 (simplify (list '(%bessel_y) v z))))))
605
606 #+nil
607 (defun simp-bessel-y (expr ignored z)
608 (declare (ignore ignored))
609 (twoargcheck expr)
610 (let ((order (simpcheck (cadr expr) z))
611 (arg (simpcheck (caddr expr) z))
612 (rat-order nil))
613 (cond
614 ((zerop1 arg)
615 ;; Domain error for a zero argument.
616 (simp-domain-error
617 (intl:gettext "bessel_y: bessel_y(~:M,~:M) is undefined.") order arg))
618
619 ((complex-float-numerical-eval-p order arg)
620 ;; We have numeric order and arg and $numer is true, or
621 ;; we have either the order or arg being floating-point,
622 ;; so let's evaluate it numerically.
623 (cond ((= 0 ($imagpart order))
624 ;; order is real, arg is real or complex
625 (complexify (bessel-y ($float order)
626 (complex ($float ($realpart arg))
627 ($float ($imagpart arg))))))
628 (t
629 ;; order is complex, arg is real or complex
630 (let (($numer t)
631 ($float t)
632 (order ($float order))
633 (arg ($float arg)))
634 ($float
635 ($rectform
636 (bessel-y-hypergeometric order arg)))))))
637
638 ((and (integerp order) (minusp order))
639 ;; Special case when the order is an integer.
640 ;; A&S 9.1.5: Y[-n](x) = (-1)^n*Y[n](x)
641 (if (evenp order)
642 (take '(%bessel_y) (- order) arg)
643 (mul -1 (take '(%bessel_y) (- order) arg))))
644
645 ((and $besselexpand
646 (setq rat-order (max-numeric-ratio-p order 2)))
647 ;; When order is a fraction with a denominator of 2, we
648 ;; can express the result in terms of elementary functions.
649 ;; From A&S 10.1.1, 10.1.11-12 and 10.1.15:
650 ;;
651 ;; Y[1/2](z) = -J[-1/2](z) is a function of cos(z).
652 ;; Y[-1/2](z) = J[1/2](z) is a function of sin(z).
653 (bessel-y-half-order rat-order arg))
654
655 ((and $bessel_reduce
656 (and (integerp order)
657 (plusp order)
658 (> order 1)))
659 ;; Reduce a bessel function of order > 2 to order 1 and 0.
660 ;; A&S 9.1.27: bessel_y(v,z) = 2*(v-1)/z*bessel_y(v-1,z)-bessel_y(v-2,z)
661 (sub (mul 2
662 (- order 1)
663 (inv arg)
664 (take '(%bessel_y) (- order 1) arg))
665 (take '(%bessel_y) (- order 2) arg)))
666
667 ($hypergeometric_representation
668 (bessel-y-hypergeometric order arg))
669
670 (t
671 (eqtest (list '(%bessel_y) order arg) expr)))))
672
673 (def-simplifier bessel_y (order arg)
674 (let ((rat-order nil))
675 (cond
676 ((zerop1 arg)
677 ;; Domain error for a zero argument.
678 (simp-domain-error
679 (intl:gettext "bessel_y: bessel_y(~:M,~:M) is undefined.") order arg))
680
681 ((complex-float-numerical-eval-p order arg)
682 ;; We have numeric order and arg and $numer is true, or
683 ;; we have either the order or arg being floating-point,
684 ;; so let's evaluate it numerically.
685 (cond ((= 0 ($imagpart order))
686 ;; order is real, arg is real or complex
687 (complexify (bessel-y ($float order)
688 (complex ($float ($realpart arg))
689 ($float ($imagpart arg))))))
690 (t
691 ;; order is complex, arg is real or complex
692 (let (($numer t)
693 ($float t)
694 (order ($float order))
695 (arg ($float arg)))
696 ($float
697 ($rectform
698 (bessel-y-hypergeometric order arg)))))))
699
700 ((and (integerp order) (minusp order))
701 ;; Special case when the order is an integer.
702 ;; A&S 9.1.5: Y[-n](x) = (-1)^n*Y[n](x)
703 (if (evenp order)
704 (take '(%bessel_y) (- order) arg)
705 (mul -1 (take '(%bessel_y) (- order) arg))))
706
707 ((and $besselexpand
708 (setq rat-order (max-numeric-ratio-p order 2)))
709 ;; When order is a fraction with a denominator of 2, we
710 ;; can express the result in terms of elementary functions.
711 ;; From A&S 10.1.1, 10.1.11-12 and 10.1.15:
712 ;;
713 ;; Y[1/2](z) = -J[-1/2](z) is a function of cos(z).
714 ;; Y[-1/2](z) = J[1/2](z) is a function of sin(z).
715 (bessel-y-half-order rat-order arg))
716
717 ((and $bessel_reduce
718 (and (integerp order)
719 (plusp order)
720 (> order 1)))
721 ;; Reduce a bessel function of order > 2 to order 1 and 0.
722 ;; A&S 9.1.27: bessel_y(v,z) = 2*(v-1)/z*bessel_y(v-1,z)-bessel_y(v-2,z)
723 (sub (mul 2
724 (- order 1)
725 (inv arg)
726 (take '(%bessel_y) (- order 1) arg))
727 (take '(%bessel_y) (- order 2) arg)))
728
729 ($hypergeometric_representation
730 (bessel-y-hypergeometric order arg))
731
732 (t
733 (give-up)))))
734
735 ;; Returns the hypergeometric representation of bessel_y
736 (defun bessel-y-hypergeometric (order arg)
737 ;; http://functions.wolfram.com/03.03.26.0002.01
738 ;; -hypergeometric([],[1-v],-z^2/4)*(2/z)^v*gamma(v)/%pi
739 ;; - hypergeometric([],[v+1],-z^2/4)*gamma(-v)*cos(%pi*v)*(z/2)^v/%pi
740 (add (mul -1
741 (power 2 order)
742 (inv (power arg order))
743 (take '(%gamma) order)
744 (inv '$%pi)
745 (take '($hypergeometric)
746 (list '(mlist))
747 (list '(mlist) (sub 1 order))
748 (neg (div (mul arg arg) 4))))
749 (mul -1
750 (inv (power 2 order))
751 (power arg order)
752 (take '(%cos) (mul order '$%pi))
753 (take '(%gamma) (neg order))
754 (inv '$%pi)
755 (take '($hypergeometric)
756 (list '(mlist))
757 (list '(mlist) (add 1 order))
758 (neg (div (mul arg arg) 4))))))
759
760 ;; Bessel function of the second kind, Y[n](z), for real or complex z
761 (defun bessel-y (order arg)
762 (cond
763 ((zerop (imagpart arg))
764 ;; We have numeric args and the first arg is purely
765 ;; real. Call the real-valued Bessel function.
766 ;;
767 ;; For negative values, use the analytic continuation formula
768 ;; A&S 9.1.36:
769 ;;
770 ;; Y[v](z*exp(m*%pi*%i)) = exp(-v*m*%pi*%i) * Y[v](z)
771 ;; + 2*%i*sin(m*v*%pi)*cot(v*%pi)*J[v](z)
772 ;;
773 ;; In particular for m = 1:
774 ;; Y[v](-z) = exp(-v*%pi*%i) * Y[v](z) + 2*%i*cos(v*%pi)*J[v](z)
775 (let ((arg (realpart arg)))
776 (cond
777 ((zerop order)
778 (cond ((>= arg 0)
779 (slatec:dbesy0 (float arg)))
780 (t
781 ;; For v = 0, this simplifies to
782 ;; Y[0](-z) = Y[0](z) + 2*%i*J[0](z)
783 ;; the return value has to be a CL number
784 (+ (slatec:dbesy0 (float (- arg)))
785 (complex 0 (* 2 (slatec:dbesj0 (float (- arg)))))))))
786 ((= order 1)
787 (cond ((>= arg 0)
788 (slatec:dbesy1 (float arg)))
789 (t
790 ;; For v = 1, this simplifies to
791 ;; Y[1](-z) = -Y[1](z) - 2*%i*J[1](v)
792 ;; the return value has to be a CL number
793 (+ (- (slatec:dbesy1 (float (- arg))))
794 (complex 0 (* -2 (slatec:dbesj1 (float (- arg)))))))))
795 ((minusp order)
796 (cond ((zerop (nth-value 1 (truncate order)))
797 ;; Order is a negative integer or float representation.
798 ;; Use A&S 9.1.5: Y[-n](z) = (-1)^n*Y[n](z).
799 (if (evenp (floor order))
800 (bessel-y (- order) arg)
801 (- (bessel-y (- order) arg))))
802 (t
803 ;; Bessel function of negative order. We use the Hankel
804 ;; functions to compute this:
805 ;; Y[v](z) = (H1[v](z) - H2[v](z)) / 2 / %i
806 (let ((result (/ (- (hankel-1 order arg)
807 (hankel-2 order arg))
808 (complex 0 2))))
809 (cond ((= (nth-value 1 (floor order)) 1/2)
810 ;; ORDER is half-integral-value or a float
811 ;; representation thereof.
812 (if (minusp arg)
813 ;; arg is negative, the result is purely imaginary
814 (complex 0 (imagpart result))
815 ;; arg is positive, the result is purely real
816 (realpart result)))
817 ;; in all other cases general complex result
818 (t result))))))
819 (t
820 (multiple-value-bind (n alpha) (floor (float order))
821 (let ((jvals (make-array (1+ n) :element-type 'flonum)))
822 ;; First we do the calculation for an positive argument.
823 (slatec:dbesy (abs (float arg)) alpha (1+ n) jvals)
824
825 ;; Now we look at the sign of the argument
826 (cond ((>= arg 0)
827 (aref jvals n))
828 (t
829 (let* ((dpi (coerce pi 'flonum))
830 (s1 (cis (- (* order dpi))))
831 (s2 (* #c(0 2) (cos (* order dpi)))))
832 (let ((result (+ (* s1 (aref jvals n))
833 (* s2 (bessel-j order (- arg))))))
834 (cond ((zerop (nth-value 1 (truncate order)))
835 ;; ORDER is an integer or a float representation
836 ;; of an integer, and the arg is positive the
837 ;; result is general complex.
838 result)
839 ((= (nth-value 1 (floor order)) 1/2)
840 ;; ORDER is a half-integral-value or an float
841 ;; representation and we have arg < 0. The
842 ;; result is purely imaginary.
843 (complex 0 (imagpart result)))
844 ;; in all other cases general complex result
845 (t result))))))))))))
846 (t
847 (cond ((minusp order)
848 ;; Bessel function of negative order. We use the Hankel function to
849 ;; compute this, because A&S 9.1.3 and 9.1.4 say
850 ;; H1[v](z) = J[v](z) + %i * Y[v](z)
851 ;; and
852 ;; H2[v](z) = J[v](z) - %i * Y[v](z).
853 ;; Thus
854 ;; Y[v](z) = (H1[v](z) - H2[v](z)) / 2 / %i
855 (/ (- (hankel-1 order arg) (hankel-2 order arg))
856 (complex 0 2)))
857 (t
858 (multiple-value-bind (n alpha) (floor (float order))
859 (let ((cyr (make-array (1+ n) :element-type 'flonum))
860 (cyi (make-array (1+ n) :element-type 'flonum))
861 (cwrkr (make-array (1+ n) :element-type 'flonum))
862 (cwrki (make-array (1+ n) :element-type 'flonum)))
863 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi
864 v-nz v-cwrkr v-cwrki v-ierr)
865 (slatec::zbesy (float (realpart arg))
866 (float (imagpart arg))
867 alpha 1 (1+ n) cyr cyi 0 cwrkr cwrki 0)
868 (declare (ignore v-zr v-zi v-fnu v-kode v-n
869 v-cyr v-cyi v-cwrkr v-cwrki v-nz))
870
871 ;; We should check for errors based on the value of v-ierr.
872 (when (plusp v-ierr)
873 (format t "zbesy ierr = ~A~%" v-ierr))
874
875 (complex (aref cyr n) (aref cyi n))))))))))
876
877 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
878 ;;;
879 ;;; Implementation of the Bessel I function
880 ;;;
881 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
882
883 (defmfun $bessel_i (v z)
884 (simplify (list '(%bessel_i) v z)))
885
886 #||
887 (defprop $bessel_i %bessel_i alias)
888 (defprop $bessel_i %bessel_i verb)
889 (defprop %bessel_i $bessel_i reversealias)
890 (defprop %bessel_i $bessel_i noun)
891
892 (defprop %bessel_i simp-bessel-i operators)
893 ||#
894
895 ;; Bessel I distributes over lists, matrices, and equations
896
897 (defprop %bessel_i (mlist $matrix mequal) distribute_over)
898
899 (defprop %bessel_i
900 ((n x)
901 ;; Derivative wrt order n. A&S 9.6.42.
902 ;;
903 ;; I[n](x)*log(x/2)
904 ;; - (x/2)^n*sum(psi(n+k+1)/gamma(n+k+1)*(x^2/4)^k/k!,k,0,inf)
905 ((mplus)
906 ((mtimes)
907 ((%bessel_i) n x)
908 ((%log) ((mtimes) ((rat) 1 2) x)))
909 ((mtimes) -1
910 ((mexpt) ((mtimes) x ((rat) 1 2)) n)
911 ((%sum)
912 ((mtimes)
913 ((mexpt) ((mfactorial) $%k) -1)
914 ((mqapply) (($psi array) 0) ((mplus) 1 $%k n))
915 ((mexpt) ((%gamma) ((mplus) 1 $%k n)) -1)
916 ((mexpt) ((mtimes) x x ((rat) 1 4)) $%k))
917 $%k 0 $inf)))
918 ;; Derivative wrt to x. A&S 9.6.29.
919 ((mtimes)
920 ((mplus) ((%bessel_i) ((mplus) -1 n) x)
921 ((%bessel_i) ((mplus) 1 n) x))
922 ((rat) 1 2)))
923 grad)
924
925 ;; Integral of the Bessel I function wrt z
926 ;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselI/21/01/01/
927 (defun bessel-i-integral-2 (v z)
928 (case v
929 (0
930 ;; integrate(bessel_i(0,z),z)
931 ;; = (1/2)*z*(bessel_i(0,z)*(%pi*struve_l(1,z)+2)
932 ;; -%pi*bessel_i(1,z)*struve_l(0,z))
933 `((mtimes) ((rat) 1 2) ,z
934 ((mplus)
935 ((mtimes) -1 $%pi
936 ((%bessel_i) 1 ,z)
937 ((%struve_l) 0 ,z))
938 ((mtimes)
939 ((%bessel_i) 0 ,z)
940 ((mplus) 2
941 ((mtimes) $%pi ((%struve_l) 1 ,z)))))))
942 (1
943 ;; integrate(bessel_i(1,z),z) = bessel_i(0,z)
944 `((%bessel_i) 0 ,z))
945 (otherwise
946 ;; http://functions.wolfram.com/03.02.21.0002.01
947 ;; integrate(bessel_i(v,z),z)
948 ;; = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2)
949 ;; * hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],z^2/4)
950 ;; = 2^(-v)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],z^2/4)
951 ;; / gamma(v+2)
952 `((mtimes)
953 (($hypergeometric)
954 ((mlist)
955 ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) ,v)))
956 ((mlist)
957 ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) ,v))
958 ((mplus) 1 ,v))
959 ((mtimes) ((rat) 1 4) ((mexpt) ,z 2)))
960 ((mexpt) 2 ((mtimes) -1 ,v))
961 ((mexpt) ((%gamma) ((mplus) 2 ,v)) -1)
962 ((mexpt) z ((mplus) 1 ,v))))))
963
964 (putprop '%bessel_i `((v z) nil ,#'bessel-i-integral-2) 'integral)
965
966 ;; Support a simplim%bessel_i function to handle specific limits
967
968 (defprop %bessel_i simplim%bessel_i simplim%function)
969
970 (defun simplim%bessel_i (expr var val)
971 ;; Look for the limit of the arguments.
972 (let ((v (limit (cadr expr) var val 'think))
973 (z (limit (caddr expr) var val 'think)))
974 (cond
975 ;; Handle an argument 0 at this place.
976 ((or (zerop1 z)
977 (eq z '$zeroa)
978 (eq z '$zerob))
979 (let ((sgn ($sign ($realpart v))))
980 (cond ((and (eq sgn '$neg)
981 (not (maxima-integerp v)))
982 ;; bessel_i(v,0), Re(v)<0 and v not an integer
983 '$infinity)
984 ((and (eq sgn '$zero)
985 (not (zerop1 v)))
986 ;; bessel_i(v,0), Re(v)=0 and v #0
987 '$und)
988 ;; Call the simplifier of the function.
989 (t (simplify (list '(%bessel_i) v z))))))
990 ((eq z '$inf)
991 ;; bessel_i(v,inf)
992 '$inf)
993 ((eq z '$minf)
994 ;; bessel_i(v,minf)
995 '$infinity)
996 (t
997 ;; All other cases are handled by the simplifier of the function.
998 (simplify (list '(%bessel_i) v z))))))
999
1000 #+nil
1001 (defun simp-bessel-i (expr ignored z)
1002 (declare (ignore ignored))
1003 (twoargcheck expr)
1004 (let ((order (simpcheck (cadr expr) z))
1005 (arg (simpcheck (caddr expr) z))
1006 (rat-order nil))
1007 (cond
1008 ((zerop1 arg)
1009 ;; We handle the different case for zero arg carefully.
1010 (let ((sgn ($sign ($realpart order))))
1011 (cond ((and (eq sgn '$zero)
1012 (zerop1 ($imagpart order)))
1013 ;; bessel_i(0,0) = 1
1014 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1))
1015 ((or (floatp order) (floatp arg)) 1.0)
1016 (t 1)))
1017 ((or (eq sgn '$pos)
1018 (maxima-integerp order))
1019 ;; bessel_i(v,0) and Re(v)>0 or v an integer
1020 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0))
1021 ((or (floatp order) (floatp arg)) 0.0)
1022 (t 0)))
1023 ((and (eq sgn '$neg)
1024 (not (maxima-integerp order)))
1025 ;; bessel_i(v,0) and Re(v)<0 and v not an integer
1026 (simp-domain-error
1027 (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.")
1028 order arg))
1029 ((and (eq sgn '$zero)
1030 (not (zerop1 ($imagpart order))))
1031 ;; bessel_i(v,0) and Re(v)=0 and v # 0
1032 (simp-domain-error
1033 (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.")
1034 order arg))
1035 (t
1036 ;; No information about the sign of the order
1037 (eqtest (list '(%bessel_i) order arg) expr)))))
1038
1039 ((complex-float-numerical-eval-p order arg)
1040 (cond ((= 0 ($imagpart order))
1041 ;; order is real, arg is real or complex
1042 (complexify (bessel-i ($float order)
1043 (complex ($float ($realpart arg))
1044 ($float ($imagpart arg))))))
1045 (t
1046 ;; order is complex, arg is real or complex
1047 (let (($numer t)
1048 ($float t)
1049 (order ($float order))
1050 (arg ($float arg)))
1051 ($float
1052 ($rectform
1053 (bessel-i-hypergeometric order arg)))))))
1054
1055 ((and (integerp order) (minusp order))
1056 ;; Some special cases when the order is an integer
1057 ;; A&S 9.6.6: I[-n](x) = I[n](x)
1058 (take '(%bessel_i) (- order) arg))
1059
1060 ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
1061 ;; When order is a fraction with a denominator of 2, we
1062 ;; can express the result in terms of elementary functions.
1063 ;; From A&S 10.2.13 and 10.2.14:
1064 ;;
1065 ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
1066 ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
1067 (bessel-i-half-order rat-order arg))
1068
1069 ((and $bessel_reduce
1070 (and (integerp order)
1071 (plusp order)
1072 (> order 1)))
1073 ;; Reduce a bessel function of order > 2 to order 1 and 0.
1074 ;; A&S 9.6.26: bessel_i(v,z) = -2*(v-1)/z*bessel_i(v-1,z)+bessel_i(v-2,z)
1075 (add (mul -2
1076 (- order 1)
1077 (inv arg)
1078 (take '(%bessel_i) (- order 1) arg))
1079 (take '(%bessel_i) (- order 2) arg)))
1080
1081 ((and $%iargs (multiplep arg '$%i))
1082 ;; bessel_i(v, %i*x) = (%i*x)^v/(x^v) * bessel_j(v, x)
1083 ;; (From http://functions.wolfram.com/03.02.27.0002.01)
1084 (let ((x (coeff arg '$%i 1)))
1085 (mul (power (mul '$%i x) order)
1086 (inv (power x order))
1087 (take '(%bessel_j) order x))))
1088
1089 ($hypergeometric_representation
1090 (bessel-i-hypergeometric order arg))
1091
1092 (t
1093 (eqtest (list '(%bessel_i) order arg) expr)))))
1094
1095 (def-simplifier bessel_i (order arg)
1096 (let ((rat-order nil))
1097 (cond
1098 ((zerop1 arg)
1099 ;; We handle the different case for zero arg carefully.
1100 (let ((sgn ($sign ($realpart order))))
1101 (cond ((and (eq sgn '$zero)
1102 (zerop1 ($imagpart order)))
1103 ;; bessel_i(0,0) = 1
1104 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1))
1105 ((or (floatp order) (floatp arg)) 1.0)
1106 (t 1)))
1107 ((or (eq sgn '$pos)
1108 (maxima-integerp order))
1109 ;; bessel_i(v,0) and Re(v)>0 or v an integer
1110 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0))
1111 ((or (floatp order) (floatp arg)) 0.0)
1112 (t 0)))
1113 ((and (eq sgn '$neg)
1114 (not (maxima-integerp order)))
1115 ;; bessel_i(v,0) and Re(v)<0 and v not an integer
1116 (simp-domain-error
1117 (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.")
1118 order arg))
1119 ((and (eq sgn '$zero)
1120 (not (zerop1 ($imagpart order))))
1121 ;; bessel_i(v,0) and Re(v)=0 and v # 0
1122 (simp-domain-error
1123 (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.")
1124 order arg))
1125 (t
1126 ;; No information about the sign of the order
1127 (give-up)))))
1128
1129 ((complex-float-numerical-eval-p order arg)
1130 (cond ((= 0 ($imagpart order))
1131 ;; order is real, arg is real or complex
1132 (complexify (bessel-i ($float order)
1133 (complex ($float ($realpart arg))
1134 ($float ($imagpart arg))))))
1135 (t
1136 ;; order is complex, arg is real or complex
1137 (let (($numer t)
1138 ($float t)
1139 (order ($float order))
1140 (arg ($float arg)))
1141 ($float
1142 ($rectform
1143 (bessel-i-hypergeometric order arg)))))))
1144
1145 ((and (integerp order) (minusp order))
1146 ;; Some special cases when the order is an integer
1147 ;; A&S 9.6.6: I[-n](x) = I[n](x)
1148 (take '(%bessel_i) (- order) arg))
1149
1150 ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
1151 ;; When order is a fraction with a denominator of 2, we
1152 ;; can express the result in terms of elementary functions.
1153 ;; From A&S 10.2.13 and 10.2.14:
1154 ;;
1155 ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
1156 ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
1157 (bessel-i-half-order rat-order arg))
1158
1159 ((and $bessel_reduce
1160 (and (integerp order)
1161 (plusp order)
1162 (> order 1)))
1163 ;; Reduce a bessel function of order > 2 to order 1 and 0.
1164 ;; A&S 9.6.26: bessel_i(v,z) = -2*(v-1)/z*bessel_i(v-1,z)+bessel_i(v-2,z)
1165 (add (mul -2
1166 (- order 1)
1167 (inv arg)
1168 (take '(%bessel_i) (- order 1) arg))
1169 (take '(%bessel_i) (- order 2) arg)))
1170
1171 ((and $%iargs (multiplep arg '$%i))
1172 ;; bessel_i(v, %i*x) = (%i*x)^v/(x^v) * bessel_j(v, x)
1173 ;; (From http://functions.wolfram.com/03.02.27.0002.01)
1174 (let ((x (coeff arg '$%i 1)))
1175 (mul (power (mul '$%i x) order)
1176 (inv (power x order))
1177 (take '(%bessel_j) order x))))
1178
1179 ($hypergeometric_representation
1180 (bessel-i-hypergeometric order arg))
1181
1182 (t
1183 (give-up)))))
1184
1185 ;; Returns the hypergeometric representation of bessel_i
1186 (defun bessel-i-hypergeometric (order arg)
1187 ;; A&S 9.6.47 / http://functions.wolfram.com/03.02.26.0002.01
1188 ;; hypergeometric([],[v+1],z^2/4)*(z/2)^v/gamma(v+1)
1189 (mul (inv (take '(%gamma) (add order 1)))
1190 (power (div arg 2) order)
1191 (take '($hypergeometric)
1192 (list '(mlist))
1193 (list '(mlist) (add order 1))
1194 (div (mul arg arg) 4))))
1195
1196 ;; Compute value of Modified Bessel function of the first kind of order ORDER
1197 (defun bessel-i (order arg)
1198 (cond
1199 ((zerop (imagpart arg))
1200 ;; We have numeric args and the first arg is purely
1201 ;; real. Call the real-valued Bessel function. Use special
1202 ;; routines for order 0 and 1, when possible
1203 (let ((arg (realpart arg)))
1204 (cond
1205 ((zerop order)
1206 (slatec:dbesi0 (float arg)))
1207 ((= order 1)
1208 (slatec:dbesi1 (float arg)))
1209 ((or (minusp order) (< arg 0))
1210 (multiple-value-bind (order-int order-frac) (floor order)
1211 (cond ((zerop order-frac)
1212 ;; Order is an integer. From A&S 9.6.6 and 9.6.30 we have
1213 ;; I[-n](z) = I[n](z)
1214 ;; and
1215 ;; I[n](-z) = (-1)^n*I[n](z)
1216 (if (< arg 0)
1217 (if (evenp order-int)
1218 (bessel-i (abs order) (abs arg))
1219 (- (bessel-i (abs order) (abs arg))))
1220 (bessel-i (abs order) arg)))
1221 (t
1222 ;; Order or arg is negative and order is not an integer, use
1223 ;; the bessel-j function for calculation. We know from
1224 ;; the definition I[v](x) = x^v/(%i*x)^v * J[v](%i*x).
1225 (let* ((arg (float arg))
1226 (result (* (expt arg order)
1227 (expt (complex 0 arg) (- order))
1228 (bessel-j order (complex 0 arg)))))
1229 ;; Try to clean up result if we know the result is
1230 ;; purely real or purely imaginary.
1231 (cond ((>= arg 0)
1232 ;; Result is purely real for arg >= 0
1233 (realpart result))
1234 ((zerop order-frac)
1235 ;; Order is an integer or a float representation of
1236 ;; an integer, the result is purely real.
1237 (realpart result))
1238 ((= order-frac 1/2)
1239 ;; Order is half-integral-value or a float
1240 ;; representation and arg < 0, the result
1241 ;; is purely imaginary.
1242 (complex 0 (imagpart result)))
1243 (t result)))))))
1244 (t
1245 ;; Now the case order > 0 and arg >= 0
1246 (multiple-value-bind (n alpha) (floor (float order))
1247 (let ((jvals (make-array (1+ n) :element-type 'flonum)))
1248 (slatec:dbesi (float (realpart arg)) alpha 1 (1+ n) jvals 0)
1249 (aref jvals n)))))))
1250
1251 ((and (zerop (realpart arg))
1252 (zerop (rem order 1)))
1253 ;; Handle the case for a pure imaginary arg and integer order.
1254 ;; In this case, the result is purely real or purely imaginary,
1255 ;; so we want to make sure that happens.
1256 ;;
1257 ;; bessel_i(n, %i*x) = (%i*x)^n/x^n * bessel_j(n,x)
1258 ;; = %i^n * bessel_j(n,x)
1259 (let* ((n (floor order))
1260 (result (bessel-j (float n) (imagpart arg))))
1261 (cond ((evenp n)
1262 ;; %i^(2*m) = (-1)^m, where n=2*m
1263 (if (evenp (/ n 2))
1264 result
1265 (- result)))
1266 ((oddp n)
1267 ;; %i^(2*m+1) = %i*(-1)^m, where n = 2*m+1
1268 (if (evenp (floor n 2))
1269 (complex 0 result)
1270 (complex 0 (- result)))))))
1271
1272 (t
1273 ;; The arg is complex. Use the complex-valued Bessel function.
1274 (multiple-value-bind (n alpha) (floor (abs (float order)))
1275 ;; Evaluate the function for positive order and fixup the result later.
1276 (let ((cyr (make-array (1+ n) :element-type 'flonum))
1277 (cyi (make-array (1+ n) :element-type 'flonum)))
1278 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
1279 v-cyr v-cyi v-nz v-ierr)
1280 (slatec::zbesi (float (realpart arg))
1281 (float (imagpart arg))
1282 alpha 1 (1+ n) cyr cyi 0 0)
1283 (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz))
1284 ;; We should check for errors here based on the value of v-ierr.
1285 (when (plusp v-ierr)
1286 (format t "zbesi ierr = ~A~%" v-ierr))
1287
1288 ;; We have evaluated I[|order|](arg), now we look at
1289 ;; the the sign of the order.
1290 (cond ((minusp order)
1291 ;; From A&S 9.6.2:
1292 ;; I[-a](z) = I[a](z) + (2/%pi)*sin(%pi*a)*K[a](z)
1293 (+ (complex (aref cyr n) (aref cyi n))
1294 (let ((dpi (coerce pi 'flonum)))
1295 (* (/ 2.0 dpi)
1296 (sin (* dpi (- order)))
1297 (bessel-k (- order) arg)))))
1298 (t
1299 (complex (aref cyr n) (aref cyi n))))))))))
1300
1301 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1302 ;;;
1303 ;;; Implementation of the Bessel K function
1304 ;;;
1305 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1306
1307 (defmfun $bessel_k (v z)
1308 (simplify (list '(%bessel_k) v z)))
1309
1310 #||
1311 (defprop $bessel_k %bessel_k alias)
1312 (defprop $bessel_k %bessel_k verb)
1313 (defprop %bessel_k $bessel_k reversealias)
1314 (defprop %bessel_k $bessel_k noun)
1315
1316 (defprop %bessel_k simp-bessel-k operators)
1317 ||#
1318
1319 ;; Bessel K distributes over lists, matrices, and equations
1320
1321 (defprop %bessel_k (mlist $matrix mequal) distribute_over)
1322
1323 (defprop %bessel_k
1324 ((n x)
1325 ;; Derivativer wrt order n. A&S 9.6.43.
1326 ;;
1327 ;; %pi/2*csc(n*%pi)*['diff(bessel_i(-n,x),n)-'diff(bessel_i(n,x),n)]
1328 ;; - %pi*cot(n*%pi)*bessel_k(n,x)
1329 ((mplus)
1330 ((mtimes) -1 $%pi
1331 ((%bessel_k) n x)
1332 ((%cot) ((mtimes) $%pi n)))
1333 ((mtimes)
1334 ((rat) 1 2)
1335 $%pi
1336 ((%csc) ((mtimes) $%pi n))
1337 ((mplus)
1338 ((%derivative) ((%bessel_i) ((mtimes) -1 n) x) n 1)
1339 ((mtimes) -1
1340 ((%derivative) ((%bessel_i) n x) n 1)))))
1341 ;; Derivative wrt to x. A&S 9.6.29.
1342 ((mtimes)
1343 -1
1344 ((mplus) ((%bessel_k) ((mplus) -1 n) x)
1345 ((%bessel_k) ((mplus) 1 n) x))
1346 ((rat) 1 2)))
1347 grad)
1348
1349 ;; Integral of the Bessel K function wrt z
1350 ;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/21/01/01/
1351 (defun bessel-k-integral-2 (n z)
1352 (cond
1353 ((and ($integerp n) (<= 0 n))
1354 (cond
1355 (($oddp n)
1356 ;; integrate(bessel_k(2*N+1,z),z) , N > 0
1357 ;; = -(-1)^((n-1)/2)*bessel_k(0,z)
1358 ;; + 2*sum((-1)^(k+(n-1)/2-1)*bessel_k(2*k,z),k,1,(n-1)/2)
1359 (let* ((k (gensym))
1360 (answer `((mplus)
1361 ((mtimes) -1 ((%bessel_k) 0 ,z)
1362 ((mexpt) -1
1363 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n))))
1364 ((mtimes) 2
1365 ((%sum)
1366 ((mtimes) ((%bessel_k) ((mtimes) 2 ,k) ,z)
1367 ((mexpt) -1
1368 ((mplus) -1 ,k
1369 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))
1370 ,k 1 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))))
1371 ;; Expand out the sum if n < 10. Otherwise fix up the indices
1372 (if (< n 10)
1373 (meval `(($ev) ,answer $sum)) ; Is there a better way?
1374 (simplify ($niceindices answer)))))
1375 (($evenp n)
1376 ;; integrate(bessel_k(2*N,z),z) , N > 0
1377 ;; = (1/2)*(-1)^(n/2)*%pi*z*(bessel_k(0,z)*struve_l(-1,z)
1378 ;; +bessel_k(1,z)*struve_l(0,z))
1379 ;; + 2 * sum((-1)^(k+n/2)*bessel_k(2*k+1,z),k,0,n/2-1)
1380 (let*
1381 ((k (gensym))
1382 (answer `((mplus)
1383 ((mtimes) 2
1384 ((%sum)
1385 ((mtimes)
1386 ((%bessel_k) ((mplus) 1 ((mtimes) 2 ,k)) ,z)
1387 ((mexpt) -1
1388 ((mplus) ,k ((mtimes) ((rat) 1 2) ,n))))
1389 ,k 0 ((mplus) -1 ((mtimes) ((rat) 1 2) ,n))))
1390 ((mtimes)
1391 ((rat) 1 2)
1392 $%pi
1393 ((mexpt) -1 ((mtimes) ((rat) 1 2) ,n))
1394 z
1395 ((mplus)
1396 ((mtimes)
1397 ((%bessel_k) 0 ,z)
1398 ((%struve_l) -1 ,z))
1399 ((mtimes)
1400 ((%bessel_k) 1 ,z)
1401 ((%struve_l) 0 ,z)))))))
1402 ;; expand out the sum if n < 10. Otherwise fix up the indices
1403 (if (< n 10)
1404 (meval `(($ev) ,answer $sum)) ; Is there a better way?
1405 (simplify ($niceindices answer)))))))
1406 (t nil)))
1407
1408 (putprop '%bessel_k `((n z) nil ,#'bessel-k-integral-2) 'integral)
1409
1410 ;; Support a simplim%bessel_k function to handle specific limits
1411
1412 (defprop %bessel_k simplim%bessel_k simplim%function)
1413
1414 (defun simplim%bessel_k (expr var val)
1415 ;; Look for the limit of the arguments.
1416 (let ((v (limit (cadr expr) var val 'think))
1417 (z (limit (caddr expr) var val 'think)))
1418 (cond
1419 ;; Handle an argument 0 at this place.
1420 ((or (zerop1 z)
1421 (eq z '$zeroa)
1422 (eq z '$zerob))
1423 (cond ((zerop1 v)
1424 ;; bessel_k(0,0)
1425 '$inf)
1426 ((integerp v)
1427 ;; bessel_k(n,0), n an integer
1428 (cond ((evenp v) '$inf)
1429 (t (cond ((eq z '$zeroa) '$inf)
1430 ((eq z '$zerob) '$minf)
1431 (t '$infinity)))))
1432 ((not (zerop1 ($realpart v)))
1433 ;; bessel_k(v,0), Re(v)#0
1434 '$infinity)
1435 ((and (zerop1 ($realpart v))
1436 (not (zerop1 v)))
1437 ;; bessel_k(v,0), Re(v)=0 and v#0
1438 '$und)
1439 ;; Call the simplifier of the function.
1440 (t (simplify (list '(%bessel_k) v z)))))
1441 ((eq z '$inf)
1442 ;; bessel_k(v,inf)
1443 0)
1444 ((eq z '$minf)
1445 ;; bessel_k(v,minf)
1446 '$infinity)
1447 (t
1448 ;; All other cases are handled by the simplifier of the function.
1449 (simplify (list '(%bessel_k) v z))))))
1450
1451 #+nil
1452 (defun simp-bessel-k (expr ignored z)
1453 (declare (ignore ignored))
1454 (twoargcheck expr)
1455 (let ((order (simpcheck (cadr expr) z))
1456 (arg (simpcheck (caddr expr) z))
1457 (rat-order nil))
1458 (cond
1459 ((zerop1 arg)
1460 ;; Domain error for a zero argument.
1461 (simp-domain-error
1462 (intl:gettext "bessel_k: bessel_k(~:M,~:M) is undefined.") order arg))
1463
1464 ((complex-float-numerical-eval-p order arg)
1465 (cond ((= 0 ($imagpart order))
1466 ;; order is real, arg is real or complex
1467 (complexify (bessel-k ($float order)
1468 (complex ($float ($realpart arg))
1469 ($float ($imagpart arg))))))
1470 (t
1471 ;; order is complex, arg is real or complex
1472 (let (($numer t)
1473 ($float t)
1474 (order ($float order))
1475 (arg ($float arg)))
1476 ($float
1477 ($rectform
1478 (bessel-k-hypergeometric order arg)))))))
1479
1480 ((mminusp order)
1481 ;; A&S 9.6.6: K[-v](x) = K[v](x)
1482 (take '(%bessel_k) (mul -1 order) arg))
1483
1484 ((and $besselexpand
1485 (setq rat-order (max-numeric-ratio-p order 2)))
1486 ;; When order is a fraction with a denominator of 2, we
1487 ;; can express the result in terms of elementary
1488 ;; functions. From A&S 10.2.16 and 10.2.17:
1489 ;;
1490 ;; K[1/2](z) = sqrt(%pi/2/z)*exp(-z)
1491 ;; = K[-1/2](z)
1492 (bessel-k-half-order rat-order arg))
1493
1494 ((and $bessel_reduce
1495 (and (integerp order)
1496 (plusp order)
1497 (> order 1)))
1498 ;; Reduce a bessel function of order > 2 to order 1 and 0.
1499 ;; A&S 9.6.26: bessel_k(v,z) = 2*(v-1)/z*bessel_k(v-1,z)+bessel_k(v-2,z)
1500 (add (mul 2
1501 (- order 1)
1502 (inv arg)
1503 (take '(%bessel_k) (- order 1) arg))
1504 (take '(%bessel_k) (- order 2) arg)))
1505
1506 ($hypergeometric_representation
1507 (bessel-k-hypergeometric order arg))
1508
1509 (t
1510 (eqtest (list '(%bessel_k) order arg) expr)))))
1511
1512 (def-simplifier bessel_k (order arg)
1513 (let ((rat-order nil))
1514 (cond
1515 ((zerop1 arg)
1516 ;; Domain error for a zero argument.
1517 (simp-domain-error
1518 (intl:gettext "bessel_k: bessel_k(~:M,~:M) is undefined.") order arg))
1519
1520 ((complex-float-numerical-eval-p order arg)
1521 (cond ((= 0 ($imagpart order))
1522 ;; order is real, arg is real or complex
1523 (complexify (bessel-k ($float order)
1524 (complex ($float ($realpart arg))
1525 ($float ($imagpart arg))))))
1526 (t
1527 ;; order is complex, arg is real or complex
1528 (let (($numer t)
1529 ($float t)
1530 (order ($float order))
1531 (arg ($float arg)))
1532 ($float
1533 ($rectform
1534 (bessel-k-hypergeometric order arg)))))))
1535
1536 ((mminusp order)
1537 ;; A&S 9.6.6: K[-v](x) = K[v](x)
1538 (take '(%bessel_k) (mul -1 order) arg))
1539
1540 ((and $besselexpand
1541 (setq rat-order (max-numeric-ratio-p order 2)))
1542 ;; When order is a fraction with a denominator of 2, we
1543 ;; can express the result in terms of elementary
1544 ;; functions. From A&S 10.2.16 and 10.2.17:
1545 ;;
1546 ;; K[1/2](z) = sqrt(%pi/2/z)*exp(-z)
1547 ;; = K[-1/2](z)
1548 (bessel-k-half-order rat-order arg))
1549
1550 ((and $bessel_reduce
1551 (and (integerp order)
1552 (plusp order)
1553 (> order 1)))
1554 ;; Reduce a bessel function of order > 2 to order 1 and 0.
1555 ;; A&S 9.6.26: bessel_k(v,z) = 2*(v-1)/z*bessel_k(v-1,z)+bessel_k(v-2,z)
1556 (add (mul 2
1557 (- order 1)
1558 (inv arg)
1559 (take '(%bessel_k) (- order 1) arg))
1560 (take '(%bessel_k) (- order 2) arg)))
1561
1562 ($hypergeometric_representation
1563 (bessel-k-hypergeometric order arg))
1564
1565 (t
1566 (give-up)))))
1567
1568 ;; Returns the hypergeometric representation of bessel_k
1569 (defun bessel-k-hypergeometric (order arg)
1570 ;; http://functions.wolfram.com/03.04.26.0002.01
1571 ;; hypergeometric([],[1-v],z^2/4)*2^(v-1)*gamma(v)/z^v
1572 ;; + hypergeometric([],[v+1],z^2/4)*2^(-v-1)*gamma(-v)*z^v
1573 (add (mul (power 2 (sub order 1))
1574 (take '(%gamma) order)
1575 (inv (power arg order))
1576 (take '($hypergeometric)
1577 (list '(mlist))
1578 (list '(mlist) (sub 1 order))
1579 (div (mul arg arg) 4)))
1580 (mul (inv (power 2 (add order 1)))
1581 (power arg order)
1582 (take '(%gamma) (neg order))
1583 (take '($hypergeometric)
1584 (list '(mlist))
1585 (list '(mlist) (add order 1))
1586 (div (mul arg arg) 4)))))
1587
1588 ;; Compute value of Modified Bessel function of the second kind of order ORDER
1589 (defun bessel-k (order arg)
1590 (cond
1591 ((zerop (imagpart arg))
1592 ;; We have numeric args and the first arg is purely real. Call the
1593 ;; real-valued Bessel function. Handle orders 0 and 1 specially, when
1594 ;; possible.
1595 (let ((arg (realpart arg)))
1596 (cond
1597 ((< arg 0)
1598 ;; This is the extension for negative arg.
1599 ;; We use A&S 9.6.6 and 9.6.31 for evaluation:
1600 ;; K[v](-z) = K[|v|](-z)
1601 ;; = exp(-%i*%pi*|v|)*K[|v|](z) - %i*%pi*I[|v|](z)
1602 (let* ((dpi (coerce pi 'flonum))
1603 (s1 (cis (* dpi (- (abs order)))))
1604 (s2 (* (complex 0 -1) dpi))
1605 (result (+ (* s1 (bessel-k (abs order) (- arg)))
1606 (* s2 (bessel-i (abs order) (- arg)))))
1607 (rem (nth-value 1 (floor order))))
1608 (cond ((zerop rem)
1609 ;; order is an integer or a float representation of an
1610 ;; integer, the result is a general complex
1611 result)
1612 ((= rem 1/2)
1613 ;; order is half-integral-value or an float representation
1614 ;; and arg < 0, the result is pure imaginary
1615 (complex 0 (imagpart result)))
1616 ;; in all other cases general complex result
1617 (t result))))
1618 ((= order 0)
1619 (slatec:dbesk0 (float arg)))
1620 ((= order 1)
1621 (slatec:dbesk1 (float arg)))
1622 (t
1623 ;; From A&S 9.6.6, K[-v](z) = K[v](z), so take the
1624 ;; absolute value of the order.
1625 (multiple-value-bind (n alpha) (floor (abs (float order)))
1626 (let ((jvals (make-array (1+ n) :element-type 'flonum)))
1627 (slatec:dbesk (float arg) alpha 1 (1+ n) jvals 0)
1628 (aref jvals n)))))))
1629 (t
1630 ;; The arg is complex. Use the complex-valued Bessel function. From
1631 ;; A&S 9.6.6, K[-v](z) = K[v](z), so take the absolute value of the order.
1632 (multiple-value-bind (n alpha) (floor (abs (float order)))
1633 (let ((cyr (make-array (1+ n) :element-type 'flonum))
1634 (cyi (make-array (1+ n) :element-type 'flonum)))
1635 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n
1636 v-cyr v-cyi v-nz v-ierr)
1637 (slatec::zbesk (float (realpart arg))
1638 (float (imagpart arg))
1639 alpha 1 (1+ n) cyr cyi 0 0)
1640 (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz))
1641
1642 ;; We should check for errors here based on the value of v-ierr.
1643 (when (plusp v-ierr)
1644 (format t "zbesk ierr = ~A~%" v-ierr))
1645 (complex (aref cyr n) (aref cyi n))))))))
1646
1647 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1648 ;;;
1649 ;;; Expansion of Bessel J and Y functions of half-integral order.
1650 ;;;
1651 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1652 ;;
1653 ;; From A&S 10.1.1, we have
1654 ;;
1655 ;; J[n+1/2](z) = sqrt(2*z/%pi)*j[n](z)
1656 ;; Y[n+1/2](z) = sqrt(2*z/%pi)*y[n](z)
1657 ;;
1658 ;; where j[n](z) is the spherical bessel function of the first kind
1659 ;; and y[n](z) is the spherical bessel function of the second kind.
1660 ;;
1661 ;; A&S 10.1.8 and 10.1.9 give
1662 ;;
1663 ;; j[n](z) = 1/z*[P(n+1/2,z)*sin(z-n*%pi/2) + Q(n+1/2,z)*cos(z-n*%pi/2)]
1664 ;;
1665 ;; y[n](z) = (-1)^(n+1)*1/z*[P(n+1/2,z)*cos(z+n*%pi/2) - Q(n+1/2,z)*sin(z+n*%pi/2)]
1666 ;;
1667
1668 ;; A&S 10.1.10
1669 ;;
1670 ;; j[n](z) = f[n](z)*sin(z) + (-1)^(n+1)*f[-n-1](z)*cos(z)
1671 ;;
1672 ;; f[0](z) = 1/z, f[1](z) = 1/z^2
1673 ;;
1674 ;; f[n-1](z) + f[n+1](z) = (2*n+1)/z*f[n](z)
1675 ;;
1676 (defun f-fun (n z)
1677 (cond ((= n 0)
1678 (div 1 z))
1679 ((= n 1)
1680 (div 1 (mul z z)))
1681 ((= n -1)
1682 0)
1683 ((>= n 2)
1684 ;; f[n+1](z) = (2*n+1)/z*f[n](z) - f[n-1](z) or
1685 ;; f[n](z) = (2*n-1)/z*f[n-1](z) - f[n-2](z)
1686 (sub (mul (div (+ n n -1) z)
1687 (f-fun (1- n) z))
1688 (f-fun (- n 2) z)))
1689 (t
1690 ;; Negative n
1691 ;;
1692 ;; f[n-1](z) = (2*n+1)/z*f[n](z) - f[n+1](z) or
1693 ;; f[n](z) = (2*n+3)/z*f[n+1](z) - f[n+2](z)
1694 (sub (mul (div (+ n n 3) z)
1695 (f-fun (1+ n) z))
1696 (f-fun (+ n 2) z)))))
1697
1698 ;; Compute sqrt(2*z/%pi)
1699 (defun root-2z/pi (z)
1700 (power (mul 2 z (inv '$%pi)) '((rat simp) 1 2)))
1701
1702 (defun bessel-j-half-order (order arg)
1703 "Compute J[n+1/2](z)"
1704 (let* ((n (floor order))
1705 (sign (if (oddp n) -1 1))
1706 (jn (sub (mul ($expand (f-fun n arg))
1707 (take '(%sin) arg))
1708 (mul sign
1709 ($expand (f-fun (- (- n) 1) arg))
1710 (take '(%cos) arg)))))
1711 (mul (root-2z/pi arg)
1712 jn)))
1713
1714 (defun bessel-y-half-order (order arg)
1715 "Compute Y[n+1/2](z)"
1716 ;; A&S 10.1.1:
1717 ;; Y[n+1/2](z) = sqrt(2*z/%pi)*y[n](z)
1718 ;;
1719 ;; A&S 10.1.15:
1720 ;; y[n](z) = (-1)^(n+1)*j[-n-1](z)
1721 ;;
1722 ;; So
1723 ;; Y[n+1/2](z) = sqrt(2*z/%pi)*(-1)^(n+1)*j[-n-1](z)
1724 ;; = (-1)^(n+1)*sqrt(2*z/%pi)*j[-n-1](z)
1725 ;; = (-1)^(n+1)*J[-n-1/2](z)
1726 (let* ((n (floor order))
1727 (jn (bessel-j-half-order (- (- order) 1/2) arg)))
1728 (if (evenp n)
1729 (mul -1 jn)
1730 jn)))
1731
1732 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1733 ;;;
1734 ;;; Expansion of Bessel I and K functions of half-integral order.
1735 ;;;
1736 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1737
1738 ;; See A&S 10.2.12
1739 ;;
1740 ;; I[n+1/2](z) = sqrt(2*z/%pi)*[g[n](z)*sinh(z) + g[-n-1](z)*cosh(z)]
1741 ;;
1742 ;; g[0](z) = 1/z, g[1](z) = -1/z^2
1743 ;;
1744 ;; g[n-1](z) - g[n+1](z) = (2*n+1)/z*g[n](z)
1745 ;;
1746 ;;
1747 (defun g-fun (n z)
1748 (declare (type integer n))
1749 (cond ((= n 0)
1750 (div 1 z))
1751 ((= n 1)
1752 (div -1 (mul z z)))
1753 ((>= n 2)
1754 ;; g[n](z) = g[n-2](z) - (2*n-1)/z*g[n-1](z)
1755 (sub (g-fun (- n 2) z)
1756 (mul (div (+ n n -1) z)
1757 (g-fun (- n 1) z))))
1758 (t
1759 ;; n is negative
1760 ;;
1761 ;; g[n](z) = (2*n+3)/z*g[n+1](z) + g[n+2](z)
1762 (add (mul (div (+ n n 3) z)
1763 (g-fun (+ n 1) z))
1764 (g-fun (+ n 2) z)))))
1765
1766 ;; See A&S 10.2.12
1767 ;;
1768 ;; I[n+1/2](z) = sqrt(2*z/%pi)*[g[n](z)*sinh(z) + g[-n-1](z)*cosh(z)]
1769
1770 (defun bessel-i-half-order (order arg)
1771 (let ((order (floor order)))
1772 (mul (root-2z/pi arg)
1773 (add (mul ($expand (g-fun order arg))
1774 (take '(%sinh) arg))
1775 (mul ($expand (g-fun (- (+ order 1)) arg))
1776 (take '(%cosh) arg))))))
1777
1778 ;; See A&S 10.2.15
1779 ;;
1780 ;; sqrt(%pi/2/z)*K[n+1/2](z) = (%pi/2/z)*exp(-z)*sum((n+1/2,k)/(2*z)^k,k,0,n)
1781 ;;
1782 ;; or
1783 ;;
1784 ;; K[n+1/2](z) = sqrt(%pi/2)/sqrt(z)*exp(-z)*sum((n+1/2,k)/(2*z)^k,k,0,n)
1785 ;;
1786 ;; where (A&S 10.1.9)
1787 ;;
1788 ;; (n+1/2,k) = (n+k)!/k!/(n-k)!
1789
1790 (defun k-fun (n z)
1791 (declare (type unsigned-byte n))
1792 ;; Computes the sum above
1793 (let ((sum 1)
1794 (term 1))
1795 (loop for k from 0 upto n do
1796 (setf term (mul term
1797 (div (div (* (- n k) (+ n k 1))
1798 (+ k 1))
1799 (mul 2 z))))
1800 (setf sum (add sum term)))
1801 sum))
1802
1803 (defun bessel-k-half-order (order arg)
1804 (let ((order (truncate (abs order))))
1805 (mul (mul (power (div '$%pi 2) '((rat simp) 1 2))
1806 (power arg '((rat simp) -1 2))
1807 (power '$%e (neg arg)))
1808 (k-fun (abs order) arg))))
1809
1810 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1811 ;;;
1812 ;;; Realpart und imagpart of Bessel functions
1813 ;;;
1814 ;;; We handle the special cases when we know that the Bessel function
1815 ;;; is pure real or pure imaginary. In all other cases Maxima generates
1816 ;;; a general noun form as result.
1817 ;;;
1818 ;;; To get the complex sign of the order an argument of the Bessel function
1819 ;;; the function $csign is used which calls $sign in a complex mode.
1820 ;;;
1821 ;;; This is an extension of of the algorithm of the function risplit in
1822 ;;; the file rpart.lisp. risplit looks for a risplit-function on the
1823 ;;; property list and call it if available.
1824 ;;;
1825 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1826
1827 (defvar *debug-bessel* nil)
1828
1829 ;;; Put the risplit-function for Bessel J and Bessel I on the property list
1830
1831 (defprop %bessel_j risplit-bessel-j-or-i risplit-function)
1832 (defprop %bessel_i risplit-bessel-j-or-i risplit-function)
1833
1834 ;;; realpart und imagpart for Bessel J and Bessel I function
1835
1836 (defun risplit-bessel-j-or-i (expr)
1837 (when *debug-bessel*
1838 (format t "~&RISPLIT-BESSEL-J with ~A~%" expr))
1839 (let ((order (cadr expr))
1840 (arg (caddr expr))
1841 (sign-order ($csign (cadr expr)))
1842 (sign-arg ($csign (caddr expr))))
1843
1844 (when *debug-bessel*
1845 (format t "~& : order = ~A~%" order)
1846 (format t "~& : arg = ~A~%" arg)
1847 (format t "~& : sign-order = ~A~%" sign-order)
1848 (format t "~& : sign-arg = ~A~%" sign-arg))
1849
1850 (cond
1851 ((or (member sign-order '($complex $imaginary))
1852 (eq sign-arg '$complex))
1853 ;; order or arg are complex, return general noun-form
1854 (risplit-noun expr))
1855 ((eq sign-arg '$imaginary)
1856 ;; arg is pure imaginary
1857 (cond
1858 ((or ($oddp order)
1859 ($featurep order '$odd))
1860 ;; order is an odd integer, pure imaginary noun-form
1861 (cons 0 (mul -1 '$%i expr)))
1862 ((or ($evenp order)
1863 ($featurep order '$even))
1864 ;; order is an even integer, real noun-form
1865 (cons expr 0))
1866 (t
1867 ;; order is not an odd or even integer, or Maxima can not
1868 ;; determine it, return general noun-form
1869 (risplit-noun expr))))
1870 ;; At this point order and arg are real.
1871 ;; We have to look for some special cases involing a negative arg
1872 ((or (maxima-integerp order)
1873 (member sign-arg '($pos $pz)))
1874 ;; arg is positive or order an integer, real noun-form
1875 (cons expr 0))
1876 ;; At this point we know that arg is negative or the sign is not known
1877 ((zerop1 (sub ($truncate ($multthru 2 order)) ($multthru 2 order)))
1878 ;; order is half integral
1879 (cond
1880 ((eq sign-arg '$neg)
1881 ;; order is half integral and arg negative, imaginary noun-form
1882 (cons 0 (mul -1 '$%i expr)))
1883 (t
1884 ;; the sign of arg or order is not known
1885 (risplit-noun expr))))
1886 (t
1887 ;; the sign of arg or order is not known
1888 (risplit-noun expr)))))
1889
1890 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1891
1892 ;;; Put the risplit-function for Bessel K and Bessel Y on the property list
1893
1894 (defprop %bessel_k risplit-bessel-k-or-y risplit-function)
1895 (defprop %bessel_y risplit-bessel-k-or-y risplit-function)
1896
1897 ;;; realpart und imagpart for Bessel K and Bessel Y function
1898
1899 (defun risplit-bessel-k-or-y (expr)
1900 (when *debug-bessel*
1901 (format t "~&RISPLIT-BESSEL-K with ~A~%" expr))
1902 (let ((order (cadr expr))
1903 (arg (caddr expr))
1904 (sign-order ($csign (cadr expr)))
1905 (sign-arg ($csign (caddr expr))))
1906
1907 (when *debug-bessel*
1908 (format t "~& : order = ~A~%" order)
1909 (format t "~& : arg = ~A~%" arg)
1910 (format t "~& : sign-order = ~A~%" sign-order)
1911 (format t "~& : sign-arg = ~A~%" sign-arg))
1912
1913 (cond
1914 ((or (member sign-order '($complex $imaginary))
1915 (member sign-arg '($complex '$imaginary)))
1916 (risplit-noun expr))
1917 ;; At this point order and arg are real valued.
1918 ;; We have to look for some special cases involing a negative arg
1919 ((member sign-arg '($pos $pz))
1920 ;; arg is positive
1921 (cons expr 0))
1922 ;; At this point we know that arg is negative or the sign is not known
1923 ((and (not (maxima-integerp order))
1924 (zerop1 (sub ($truncate ($multthru 2 order)) ($multthru 2 order))))
1925 ;; order is half integral
1926 (cond
1927 ((eq sign-arg '$neg)
1928 (cons 0 (mul -1 '$%i expr)))
1929 (t
1930 ;; the sign of arg is not known
1931 (risplit-noun expr))))
1932 (t
1933 ;; the sign of arg is not known
1934 (risplit-noun expr)))))
1935
1936 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1937 ;;;
1938 ;;; Implementation of scaled Bessel I functions
1939 ;;;
1940 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1941
1942 ;; FIXME: The following scaled functions need work. They should be
1943 ;; extended to match bessel_i, but still carefully compute the scaled
1944 ;; value.
1945
1946 ;; I think g0(x) = exp(-x)*I[0](x), g1(x) = exp(-x)*I[1](x), and
1947 ;; gn(x,n) = exp(-x)*I[n](x), based on some simple numerical
1948 ;; evaluations.
1949
1950 (defmfun $scaled_bessel_i0 ($x)
1951 (cond ((mnump $x)
1952 ;; XXX Should we return noun forms if $x is rational?
1953 (slatec:dbsi0e ($float $x)))
1954 (t
1955 (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil)))
1956 `((%bessel_i) 0 ,$x)))))
1957
1958 (defmfun $scaled_bessel_i1 ($x)
1959 (cond ((mnump $x)
1960 ;; XXX Should we return noun forms if $x is rational?
1961 (slatec:dbsi1e ($float $x)))
1962 (t
1963 (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil)))
1964 `((%bessel_i) 1 ,$x)))))
1965
1966 (defmfun $scaled_bessel_i ($n $x)
1967 (cond ((and (mnump $x) (mnump $n))
1968 ;; XXX Should we return noun forms if $n and $x are rational?
1969 (multiple-value-bind (n alpha) (floor ($float $n))
1970 (let ((iarray (make-array (1+ n) :element-type 'flonum)))
1971 (slatec:dbesi ($float $x) alpha 2 (1+ n) iarray 0)
1972 (aref iarray n))))
1973 (t
1974 (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil)))
1975 ($bessel_i $n $x)))))
1976
1977 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1978 ;;;
1979 ;;; Implementation of the Hankel 1 function
1980 ;;;
1981 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1982
1983 (defmfun $hankel_1 (v z)
1984 (simplify (list '(%hankel_1) v z)))
1985
1986 #||
1987 (defprop $hankel_1 %hankel_1 alias)
1988 (defprop $hankel_1 %hankel_1 verb)
1989 (defprop %hankel_1 $hankel_1 reversealias)
1990 (defprop %hankel_1 $hankel_1 noun)
1991 ||#
1992
1993 ;; hankel_1 distributes over lists, matrices, and equations
1994
1995 (defprop %hankel_1 (mlist $matrix mequal) distribute_over)
1996
1997 #+nil
1998 (defprop %hankel_1 simp-hankel-1 operators)
1999
2000 ; Derivatives of the Hankel 1 function
2001 (defprop %hankel_1
2002 ((n x)
2003 nil
2004
2005 ;; Derivative wrt arg x. A&S 9.1.27.
2006 ((mtimes)
2007 ((mplus)
2008 ((%hankel_1) ((mplus) -1 n) x)
2009 ((mtimes) -1 ((%hankel_1) ((mplus) 1 n) x)))
2010 ((rat) 1 2)))
2011 grad)
2012
2013 #+nil
2014 (defun simp-hankel-1 (expr ignored z)
2015 (declare (ignore ignored))
2016 (twoargcheck expr)
2017 (let ((order (simpcheck (cadr expr) z))
2018 (arg (simpcheck (caddr expr) z))
2019 rat-order)
2020 (cond
2021 ((zerop1 arg)
2022 (simp-domain-error
2023 (intl:gettext "hankel_1: hankel_1(~:M,~:M) is undefined.")
2024 order arg))
2025 ((complex-float-numerical-eval-p order arg)
2026 (cond ((= 0 ($imagpart order))
2027 ;; order is real, arg is real or complex
2028 (complexify (hankel-1 ($float order)
2029 (complex ($float ($realpart arg))
2030 ($float ($imagpart arg))))))
2031 (t
2032 ;; The order is complex. Use
2033 ;; hankel_1(v,z) = bessel_j(v,z) + %i*bessel_y(v,z)
2034 ;; and evaluate using the hypergeometric function
2035 (let (($numer t)
2036 ($float t)
2037 (order ($float order))
2038 (arg ($float arg)))
2039 ($float
2040 ($rectform
2041 (add (bessel-j-hypergeometric order arg)
2042 (mul '$%i
2043 (bessel-y-hypergeometric order arg)))))))))
2044 ((and $besselexpand
2045 (setq rat-order (max-numeric-ratio-p order 2)))
2046 ;; When order is a fraction with a denominator of 2, we can express
2047 ;; the result in terms of elementary functions.
2048 ;; Use the definition hankel_1(v,z) = bessel_j(v,z)+%i*bessel_y(v,z)
2049 (sratsimp
2050 (add (bessel-j-half-order rat-order arg)
2051 (mul '$%i
2052 (bessel-y-half-order rat-order arg)))))
2053 (t (eqtest (list '(%hankel_1) order arg) expr)))))
2054
2055 (def-simplifier hankel_1 (order arg)
2056 (let (rat-order)
2057 (cond
2058 ((zerop1 arg)
2059 (simp-domain-error
2060 (intl:gettext "hankel_1: hankel_1(~:M,~:M) is undefined.")
2061 order arg))
2062 ((complex-float-numerical-eval-p order arg)
2063 (cond ((= 0 ($imagpart order))
2064 ;; order is real, arg is real or complex
2065 (complexify (hankel-1 ($float order)
2066 (complex ($float ($realpart arg))
2067 ($float ($imagpart arg))))))
2068 (t
2069 ;; The order is complex. Use
2070 ;; hankel_1(v,z) = bessel_j(v,z) + %i*bessel_y(v,z)
2071 ;; and evaluate using the hypergeometric function
2072 (let (($numer t)
2073 ($float t)
2074 (order ($float order))
2075 (arg ($float arg)))
2076 ($float
2077 ($rectform
2078 (add (bessel-j-hypergeometric order arg)
2079 (mul '$%i
2080 (bessel-y-hypergeometric order arg)))))))))
2081 ((and $besselexpand
2082 (setq rat-order (max-numeric-ratio-p order 2)))
2083 ;; When order is a fraction with a denominator of 2, we can express
2084 ;; the result in terms of elementary functions.
2085 ;; Use the definition hankel_1(v,z) = bessel_j(v,z)+%i*bessel_y(v,z)
2086 (sratsimp
2087 (add (bessel-j-half-order rat-order arg)
2088 (mul '$%i
2089 (bessel-y-half-order rat-order arg)))))
2090 (t (give-up)))))
2091
2092 ;; Numerically compute H1[v](z).
2093 ;;
2094 ;; A&S 9.1.3 says H1[v](z) = J[v](z) + %i * Y[v](z)
2095 ;;
2096 (defun hankel-1 (v z)
2097 (let ((v (float v))
2098 (z (coerce z '(complex flonum))))
2099 (cond ((minusp v)
2100 ;; A&S 9.1.6:
2101 ;;
2102 ;; H1[-v](z) = exp(v*%pi*%i)*H1[v](z)
2103 ;;
2104 ;; or
2105 ;;
2106 ;; H1[v](z) = exp(-v*%pi*%i)*H1[-v](z)
2107
2108 (* (cis (* (float pi) (- v))) (hankel-1 (- v) z)))
2109 (t
2110 (multiple-value-bind (n fnu) (floor v)
2111 (let ((zr (realpart z))
2112 (zi (imagpart z))
2113 (cyr (make-array (1+ n) :element-type 'flonum))
2114 (cyi (make-array (1+ n) :element-type 'flonum)))
2115 (multiple-value-bind (dzr dzi df dk dm dn dcyr dcyi nz ierr)
2116 (slatec::zbesh zr zi fnu 1 1 (1+ n) cyr cyi 0 0)
2117 (declare (ignore dzr dzi df dk dm dn dcyr dcyi nz ierr))
2118 (complex (aref cyr n) (aref cyi n)))))))))
2119
2120 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2121 ;;;
2122 ;;; Implementation of the Hankel 2 function
2123 ;;;
2124 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2125
2126 (defmfun $hankel_2 (v z)
2127 (simplify (list '(%hankel_2) v z)))
2128
2129 #||
2130 (defprop $hankel_2 %hankel_2 alias)
2131 (defprop $hankel_2 %hankel_2 verb)
2132 (defprop %hankel_2 $hankel_2 reversealias)
2133 (defprop %hankel_2 $hankel_2 noun)
2134 ||#
2135
2136 ;; hankel_2 distributes over lists, matrices, and equations
2137
2138 (defprop %hankel_2 (mlist $matrix mequal) distribute_over)
2139
2140 #+nil
2141 (defprop %hankel_2 simp-hankel-2 operators)
2142
2143 ; Derivatives of the Hankel 2 function
2144 (defprop %hankel_2
2145 ((n x)
2146 nil
2147
2148 ;; Derivative wrt arg x. A&S 9.1.27.
2149 ((mtimes)
2150 ((mplus)
2151 ((%hankel_2) ((mplus) -1 n) x)
2152 ((mtimes) -1 ((%hankel_2) ((mplus) 1 n) x)))
2153 ((rat) 1 2)))
2154 grad)
2155
2156 #+nil
2157 (defun simp-hankel-2 (expr ignored z)
2158 (declare (ignore ignored))
2159 (twoargcheck expr)
2160 (let ((order (simpcheck (cadr expr) z))
2161 (arg (simpcheck (caddr expr) z))
2162 rat-order)
2163 (cond
2164 ((zerop1 arg)
2165 (simp-domain-error
2166 (intl:gettext "hankel_2: hankel_2(~:M,~:M) is undefined.")
2167 order arg))
2168 ((complex-float-numerical-eval-p order arg)
2169 (cond ((= 0 ($imagpart order))
2170 ;; order is real, arg is real or complex
2171 (complexify (hankel-2 ($float order)
2172 (complex ($float ($realpart arg))
2173 ($float ($imagpart arg))))))
2174 (t
2175 ;; The order is complex. Use
2176 ;; hankel_2(v,z) = bessel_j(v,z) - %i*bessel_y(v,z)
2177 ;; and evaluate using the hypergeometric function
2178 (let (($numer t)
2179 ($float t)
2180 (order ($float order))
2181 (arg ($float arg)))
2182 ($float
2183 ($rectform
2184 (sub (bessel-j-hypergeometric order arg)
2185 (mul '$%i
2186 (bessel-y-hypergeometric order arg)))))))))
2187 ((and $besselexpand
2188 (setq rat-order (max-numeric-ratio-p order 2)))
2189 ;; When order is a fraction with a denominator of 2, we can express
2190 ;; the result in terms of elementary functions.
2191 ;; Use the definition hankel_2(v,z) = bessel_j(v,z)-%i*bessel_y(v,z)
2192 (sratsimp
2193 (sub (bessel-j-half-order rat-order arg)
2194 (mul '$%i
2195 (bessel-y-half-order rat-order arg)))))
2196 (t (eqtest (list '(%hankel_2) order arg) expr)))))
2197
2198 (def-simplifier hankel_2 (order arg)
2199 (let (rat-order)
2200 (cond
2201 ((zerop1 arg)
2202 (simp-domain-error
2203 (intl:gettext "hankel_2: hankel_2(~:M,~:M) is undefined.")
2204 order arg))
2205 ((complex-float-numerical-eval-p order arg)
2206 (cond ((= 0 ($imagpart order))
2207 ;; order is real, arg is real or complex
2208 (complexify (hankel-2 ($float order)
2209 (complex ($float ($realpart arg))
2210 ($float ($imagpart arg))))))
2211 (t
2212 ;; The order is complex. Use
2213 ;; hankel_2(v,z) = bessel_j(v,z) - %i*bessel_y(v,z)
2214 ;; and evaluate using the hypergeometric function
2215 (let (($numer t)
2216 ($float t)
2217 (order ($float order))
2218 (arg ($float arg)))
2219 ($float
2220 ($rectform
2221 (sub (bessel-j-hypergeometric order arg)
2222 (mul '$%i
2223 (bessel-y-hypergeometric order arg)))))))))
2224 ((and $besselexpand
2225 (setq rat-order (max-numeric-ratio-p order 2)))
2226 ;; When order is a fraction with a denominator of 2, we can express
2227 ;; the result in terms of elementary functions.
2228 ;; Use the definition hankel_2(v,z) = bessel_j(v,z)-%i*bessel_y(v,z)
2229 (sratsimp
2230 (sub (bessel-j-half-order rat-order arg)
2231 (mul '$%i
2232 (bessel-y-half-order rat-order arg)))))
2233 (t (give-up)))))
2234
2235 ;; Numerically compute H2[v](z).
2236 ;;
2237 ;; A&S 9.1.4 says H2[v](z) = J[v](z) - %i * Y[v](z)
2238 ;;
2239 (defun hankel-2 (v z)
2240 (let ((v (float v))
2241 (z (coerce z '(complex flonum))))
2242 (cond ((minusp v)
2243 ;; A&S 9.1.6:
2244 ;;
2245 ;; H2[-v](z) = exp(-v*%pi*%i)*H2[v](z)
2246 ;;
2247 ;; or
2248 ;;
2249 ;; H2[v](z) = exp(v*%pi*%i)*H2[-v](z)
2250
2251 (* (cis (* (float pi) v)) (hankel-2 (- v) z)))
2252 (t
2253 (multiple-value-bind (n fnu) (floor v)
2254 (let ((zr (realpart z))
2255 (zi (imagpart z))
2256 (cyr (make-array (1+ n) :element-type 'flonum))
2257 (cyi (make-array (1+ n) :element-type 'flonum)))
2258 (multiple-value-bind (dzr dzi df dk dm dn dcyr dcyi nz ierr)
2259 (slatec::zbesh zr zi fnu 1 2 (1+ n) cyr cyi 0 0)
2260 (declare (ignore dzr dzi df dk dm dn dcyr dcyi nz ierr))
2261 (complex (aref cyr n) (aref cyi n)))))))))
2262
2263 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2264 ;;;
2265 ;;; Implementation of Struve H function
2266 ;;;
2267 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2268
2269 (defmfun $struve_h (v z)
2270 (simplify (list '(%struve_h) v z)))
2271
2272 #||
2273 (defprop $struve_h %struve_h alias)
2274 (defprop $struve_h %struve_h verb)
2275 (defprop %struve_h $struve_h reversealias)
2276 (defprop %struve_h $struve_h noun)
2277 ||#
2278
2279 ;; struve_h distributes over lists, matrices, and equations
2280
2281 (defprop %struve_h (mlist $matrix mequal) distribute_over)
2282
2283 #+nil
2284 (defprop %struve_h simp-struve-h operators)
2285
2286 ; Derivatives of the Struve H function
2287 (defprop %struve_h
2288 ((v z)
2289 nil
2290
2291 ;; Derivative wrt arg z. A&S 12.1.10.
2292 ((mtimes)
2293 ((rat) 1 2)
2294 ((mplus)
2295 ((%struve_h) ((mplus) -1 v) z)
2296 ((mtimes) -1 ((%struve_h) ((mplus) 1 v) z))
2297 ((mtimes)
2298 ((mexpt) $%pi ((rat) -1 2))
2299 ((mexpt) 2 ((mtimes) -1 v))
2300 ((mexpt) ((%gamma) ((mplus) ((rat) 3 2) v)) -1)
2301 ((mexpt) z v)))))
2302 grad)
2303
2304 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2305
2306 #+nil
2307 (defun simp-struve-h (expr ignored z)
2308 (declare (ignore ignored))
2309 (twoargcheck expr)
2310 (let ((order (simpcheck (cadr expr) z))
2311 (arg (simpcheck (caddr expr) z)))
2312 (cond
2313
2314 ;; Check for special values
2315
2316 ((zerop1 arg)
2317 (cond ((eq ($csign (add order 1)) '$zero)
2318 ; http://functions.wolfram.com/03.09.03.0018.01
2319 (cond ((or ($bfloatp order)
2320 ($bfloatp arg))
2321 ($bfloat (div 2 '$%pi)))
2322 ((or (floatp order)
2323 (floatp arg))
2324 ($float (div 2 '$%pi)))
2325 (t
2326 (div 2 '$%pi))))
2327 ((eq ($sign (add ($realpart order) 1)) '$pos)
2328 ; http://functions.wolfram.com/03.09.03.0001.01
2329 arg)
2330 ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz))
2331 (simp-domain-error
2332 (intl:gettext "struve_h: struve_h(~:M,~:M) is undefined.")
2333 order arg))
2334 (t
2335 (eqtest (list '(%struve_h) order arg) expr))))
2336
2337 ;; Check for numerical evaluation
2338
2339 ((complex-float-numerical-eval-p order arg)
2340 ; A&S 12.1.21
2341 (let (($numer t) ($float t))
2342 ($rectform
2343 ($float
2344 (mul
2345 ($rectform (power arg (add order 1.0)))
2346 ($rectform (inv (power 2.0 order)))
2347 (inv (power ($float '$%pi) 0.5))
2348 (inv (simplify (list '(%gamma) (add order 1.5))))
2349 (simplify (list '($hypergeometric)
2350 (list '(mlist) 1)
2351 (list '(mlist) '((rat simp) 3 2)
2352 (add order '((rat simp) 3 2)))
2353 (div (mul arg arg) -4.0))))))))
2354
2355 ((complex-bigfloat-numerical-eval-p order arg)
2356 ; A&S 12.1.21
2357 (let (($ratprint nil)
2358 (arg ($bfloat arg))
2359 (order ($bfloat order)))
2360 ($rectform
2361 ($bfloat
2362 (mul
2363 ($rectform (power arg (add order 1)))
2364 ($rectform (inv (power 2 order)))
2365 (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2))))
2366 (inv (simplify (list '(%gamma)
2367 (add order ($bfloat '((rat simp) 3 2))))))
2368 (simplify (list '($hypergeometric)
2369 (list '(mlist) 1)
2370 (list '(mlist) '((rat simp) 3 2)
2371 (add order '((rat simp) 3 2)))
2372 (div (mul arg arg) ($bfloat -4)))))))))
2373
2374 ;; Transformations and argument simplifications
2375
2376 ((and $besselexpand
2377 (ratnump order)
2378 (integerp (mul 2 order)))
2379 (cond
2380 ((eq ($sign order) '$pos)
2381 ;; Expansion of Struve H for a positive half integral order.
2382 (sratsimp
2383 (add
2384 (mul
2385 (inv (simplify (list '(mfactorial) (sub order
2386 '((rat simp) 1 2)))))
2387 (inv (power '$%pi '((rat simp) 1 2 )))
2388 (power (div arg 2) (add order -1))
2389 (let ((index (gensumindex)))
2390 (dosum
2391 (mul
2392 (simplify (list '($pochhammer) '((rat simp) 1 2) index))
2393 (simplify (list '($pochhammer)
2394 (sub '((rat simp) 1 2) order)
2395 index))
2396 (power (mul -1 arg arg (inv 4)) (mul -1 index)))
2397 index 0 (sub order '((rat simp) 1 2)) t)))
2398 (mul
2399 (power (div 2 '$%pi) '((rat simp) 1 2))
2400 (power -1 (add order '((rat simp) 1 2)))
2401 (inv (power arg '((rat simp) 1 2)))
2402 (add
2403 (mul
2404 (simplify
2405 (list '(%sin)
2406 (add (mul '((rat simp) 1 2)
2407 '$%pi
2408 (add order '((rat simp) 1 2)))
2409 arg)))
2410 (let ((index (gensumindex)))
2411 (dosum
2412 (mul
2413 (power -1 index)
2414 (simplify (list '(mfactorial)
2415 (add (mul 2 index)
2416 order
2417 '((rat simp) -1 2))))
2418 (inv (simplify (list '(mfactorial) (mul 2 index))))
2419 (inv (simplify (list '(mfactorial)
2420 (add (mul -2 index)
2421 order
2422 '((rat simp) -1 2)))))
2423 (inv (power (mul 2 arg) (mul 2 index))))
2424 index 0
2425 (simplify (list '($floor)
2426 (div (sub (mul 2 order) 1) 4)))
2427 t)))
2428 (mul
2429 (simplify (list '(%cos)
2430 (add (mul '((rat simp) 1 2)
2431 '$%pi
2432 (add order '((rat simp) 1 2)))
2433 arg)))
2434 (let ((index (gensumindex)))
2435 (dosum
2436 (mul
2437 (power -1 index)
2438 (simplify (list '(mfactorial)
2439 (add (mul 2 index)
2440 order
2441 '((rat simp) 1 2))))
2442 (power (mul 2 arg) (mul -1 (add (mul 2 index) 1)))
2443 (inv (simplify (list '(mfactorial)
2444 (add (mul 2 index) 1))))
2445 (inv (simplify (list '(mfactorial)
2446 (add (mul -2 index)
2447 order
2448 '((rat simp) -3 2))))))
2449 index 0
2450 (simplify (list '($floor)
2451 (div (sub (mul 2 order) 3) 4)))
2452 t))))))))
2453
2454 ((eq ($sign order) '$neg)
2455 ;; Expansion of Struve H for a negative half integral order.
2456 (sratsimp
2457 (add
2458 (mul
2459 (power (div 2 '$%pi) '((rat simp) 1 2))
2460 (power -1 (add order '((rat simp) 1 2)))
2461 (inv (power arg '((rat simp) 1 2)))
2462 (add
2463 (mul
2464 (simplify (list '(%sin)
2465 (add
2466 (mul
2467 '((rat simp) 1 2)
2468 '$%pi
2469 (add order '((rat simp) 1 2)))
2470 arg)))
2471 (let ((index (gensumindex)))
2472 (dosum
2473 (mul
2474 (power -1 index)
2475 (simplify (list '(mfactorial)
2476 (add (mul 2 index)
2477 (neg order)
2478 '((rat simp) -1 2))))
2479 (inv (simplify (list '(mfactorial) (mul 2 index))))
2480 (inv (simplify (list '(mfactorial)
2481 (add (mul -2 index)
2482 (neg order)
2483 '((rat simp) -1 2)))))
2484 (inv (power (mul 2 arg) (mul 2 index))))
2485 index 0
2486 (simplify (list '($floor)
2487 (div (add (mul 2 order) 1) -4)))
2488 t)))
2489 (mul
2490 (simplify (list '(%cos)
2491 (add
2492 (mul
2493 '((rat simp) 1 2)
2494 '$%pi
2495 (add order '((rat simp) 1 2)))
2496 arg)))
2497 (let ((index (gensumindex)))
2498 (dosum
2499 (mul
2500 (power -1 index)
2501 (simplify (list '(mfactorial)
2502 (add (mul 2 index)
2503 (neg order)
2504 '((rat simp) 1 2))))
2505 (power (mul 2 arg) (mul -1 (add (mul 2 index) 1)))
2506 (inv (simplify (list '(mfactorial)
2507 (add (mul 2 index) 1))))
2508 (inv (simplify (list '(mfactorial)
2509 (add (mul -2 index)
2510 (neg order)
2511 '((rat simp) -3 2))))))
2512 index 0
2513 (simplify (list '($floor)
2514 (div (add (mul 2 order) 3) -4)))
2515 t))))))))))
2516 (t
2517 (eqtest (list '(%struve_h) order arg) expr)))))
2518
2519 (def-simplifier struve_h (order arg)
2520 (cond
2521
2522 ;; Check for special values
2523
2524 ((zerop1 arg)
2525 (cond ((eq ($csign (add order 1)) '$zero)
2526 ; http://functions.wolfram.com/03.09.03.0018.01
2527 (cond ((or ($bfloatp order)
2528 ($bfloatp arg))
2529 ($bfloat (div 2 '$%pi)))
2530 ((or (floatp order)
2531 (floatp arg))
2532 ($float (div 2 '$%pi)))
2533 (t
2534 (div 2 '$%pi))))
2535 ((eq ($sign (add ($realpart order) 1)) '$pos)
2536 ; http://functions.wolfram.com/03.09.03.0001.01
2537 arg)
2538 ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz))
2539 (simp-domain-error
2540 (intl:gettext "struve_h: struve_h(~:M,~:M) is undefined.")
2541 order arg))
2542 (t
2543 (give-up))))
2544
2545 ;; Check for numerical evaluation
2546
2547 ((complex-float-numerical-eval-p order arg)
2548 ; A&S 12.1.21
2549 (let (($numer t) ($float t))
2550 ($rectform
2551 ($float
2552 (mul
2553 ($rectform (power arg (add order 1.0)))
2554 ($rectform (inv (power 2.0 order)))
2555 (inv (power ($float '$%pi) 0.5))
2556 (inv (simplify (list '(%gamma) (add order 1.5))))
2557 (simplify (list '($hypergeometric)
2558 (list '(mlist) 1)
2559 (list '(mlist) '((rat simp) 3 2)
2560 (add order '((rat simp) 3 2)))
2561 (div (mul arg arg) -4.0))))))))
2562
2563 ((complex-bigfloat-numerical-eval-p order arg)
2564 ; A&S 12.1.21
2565 (let (($ratprint nil)
2566 (arg ($bfloat arg))
2567 (order ($bfloat order)))
2568 ($rectform
2569 ($bfloat
2570 (mul
2571 ($rectform (power arg (add order 1)))
2572 ($rectform (inv (power 2 order)))
2573 (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2))))
2574 (inv (simplify (list '(%gamma)
2575 (add order ($bfloat '((rat simp) 3 2))))))
2576 (simplify (list '($hypergeometric)
2577 (list '(mlist) 1)
2578 (list '(mlist) '((rat simp) 3 2)
2579 (add order '((rat simp) 3 2)))
2580 (div (mul arg arg) ($bfloat -4)))))))))
2581
2582 ;; Transformations and argument simplifications
2583
2584 ((and $besselexpand
2585 (ratnump order)
2586 (integerp (mul 2 order)))
2587 (cond
2588 ((eq ($sign order) '$pos)
2589 ;; Expansion of Struve H for a positive half integral order.
2590 (sratsimp
2591 (add
2592 (mul
2593 (inv (simplify (list '(mfactorial) (sub order
2594 '((rat simp) 1 2)))))
2595 (inv (power '$%pi '((rat simp) 1 2 )))
2596 (power (div arg 2) (add order -1))
2597 (let ((index (gensumindex)))
2598 (dosum
2599 (mul
2600 (simplify (list '($pochhammer) '((rat simp) 1 2) index))
2601 (simplify (list '($pochhammer)
2602 (sub '((rat simp) 1 2) order)
2603 index))
2604 (power (mul -1 arg arg (inv 4)) (mul -1 index)))
2605 index 0 (sub order '((rat simp) 1 2)) t)))
2606 (mul
2607 (power (div 2 '$%pi) '((rat simp) 1 2))
2608 (power -1 (add order '((rat simp) 1 2)))
2609 (inv (power arg '((rat simp) 1 2)))
2610 (add
2611 (mul
2612 (simplify
2613 (list '(%sin)
2614 (add (mul '((rat simp) 1 2)
2615 '$%pi
2616 (add order '((rat simp) 1 2)))
2617 arg)))
2618 (let ((index (gensumindex)))
2619 (dosum
2620 (mul
2621 (power -1 index)
2622 (simplify (list '(mfactorial)
2623 (add (mul 2 index)
2624 order
2625 '((rat simp) -1 2))))
2626 (inv (simplify (list '(mfactorial) (mul 2 index))))
2627 (inv (simplify (list '(mfactorial)
2628 (add (mul -2 index)
2629 order
2630 '((rat simp) -1 2)))))
2631 (inv (power (mul 2 arg) (mul 2 index))))
2632 index 0
2633 (simplify (list '($floor)
2634 (div (sub (mul 2 order) 1) 4)))
2635 t)))
2636 (mul
2637 (simplify (list '(%cos)
2638 (add (mul '((rat simp) 1 2)
2639 '$%pi
2640 (add order '((rat simp) 1 2)))
2641 arg)))
2642 (let ((index (gensumindex)))
2643 (dosum
2644 (mul
2645 (power -1 index)
2646 (simplify (list '(mfactorial)
2647 (add (mul 2 index)
2648 order
2649 '((rat simp) 1 2))))
2650 (power (mul 2 arg) (mul -1 (add (mul 2 index) 1)))
2651 (inv (simplify (list '(mfactorial)
2652 (add (mul 2 index) 1))))
2653 (inv (simplify (list '(mfactorial)
2654 (add (mul -2 index)
2655 order
2656 '((rat simp) -3 2))))))
2657 index 0
2658 (simplify (list '($floor)
2659 (div (sub (mul 2 order) 3) 4)))
2660 t))))))))
2661
2662 ((eq ($sign order) '$neg)
2663 ;; Expansion of Struve H for a negative half integral order.
2664 (sratsimp
2665 (add
2666 (mul
2667 (power (div 2 '$%pi) '((rat simp) 1 2))
2668 (power -1 (add order '((rat simp) 1 2)))
2669 (inv (power arg '((rat simp) 1 2)))
2670 (add
2671 (mul
2672 (simplify (list '(%sin)
2673 (add
2674 (mul
2675 '((rat simp) 1 2)
2676 '$%pi
2677 (add order '((rat simp) 1 2)))
2678 arg)))
2679 (let ((index (gensumindex)))
2680 (dosum
2681 (mul
2682 (power -1 index)
2683 (simplify (list '(mfactorial)
2684 (add (mul 2 index)
2685 (neg order)
2686 '((rat simp) -1 2))))
2687 (inv (simplify (list '(mfactorial) (mul 2 index))))
2688 (inv (simplify (list '(mfactorial)
2689 (add (mul -2 index)
2690 (neg order)
2691 '((rat simp) -1 2)))))
2692 (inv (power (mul 2 arg) (mul 2 index))))
2693 index 0
2694 (simplify (list '($floor)
2695 (div (add (mul 2 order) 1) -4)))
2696 t)))
2697 (mul
2698 (simplify (list '(%cos)
2699 (add
2700 (mul
2701 '((rat simp) 1 2)
2702 '$%pi
2703 (add order '((rat simp) 1 2)))
2704 arg)))
2705 (let ((index (gensumindex)))
2706 (dosum
2707 (mul
2708 (power -1 index)
2709 (simplify (list '(mfactorial)
2710 (add (mul 2 index)
2711 (neg order)
2712 '((rat simp) 1 2))))
2713 (power (mul 2 arg) (mul -1 (add (mul 2 index) 1)))
2714 (inv (simplify (list '(mfactorial)
2715 (add (mul 2 index) 1))))
2716 (inv (simplify (list '(mfactorial)
2717 (add (mul -2 index)
2718 (neg order)
2719 '((rat simp) -3 2))))))
2720 index 0
2721 (simplify (list '($floor)
2722 (div (add (mul 2 order) 3) -4)))
2723 t))))))))))
2724 (t
2725 (give-up))))
2726
2727 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2728 ;;;
2729 ;;; Implementation of Struve L function
2730 ;;;
2731 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2732
2733 (defmfun $struve_l (v z)
2734 (simplify (list '(%struve_l) v z)))
2735
2736 #||
2737 (defprop $struve_l %struve_l alias)
2738 (defprop $struve_l %struve_l verb)
2739 (defprop %struve_l $struve_l reversealias)
2740 (defprop %struve_l $struve_l noun)
2741
2742 (defprop %struve_l simp-struve-l operators)
2743 ||#
2744
2745 ;; struve_l distributes over lists, matrices, and equations
2746
2747 (defprop %struve_l (mlist $matrix mequal) distribute_over)
2748
2749 ; Derivatives of the Struve L function
2750 (defprop %struve_l
2751 ((v z)
2752 nil
2753
2754 ;; Derivative wrt arg z. A&S 12.2.5.
2755 ((mtimes)
2756 ((rat) 1 2)
2757 ((mplus)
2758 ((%struve_l) ((mplus) -1 v) z)
2759 ((%struve_l) ((mplus) 1 v) z)
2760 ((mtimes)
2761 ((mexpt) $%pi ((rat) -1 2))
2762 ((mexpt) 2 ((mtimes) -1 v))
2763 ((mexpt) ((%gamma) ((mplus) ((rat) 3 2) v)) -1)
2764 ((mexpt) z v)))))
2765 grad)
2766
2767 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2768
2769 #+nil
2770 (defun simp-struve-l (expr ignored z)
2771 (declare (ignore ignored))
2772 (twoargcheck expr)
2773 (let ((order (simpcheck (cadr expr) z))
2774 (arg (simpcheck (caddr expr) z)))
2775 (cond
2776
2777 ;; Check for special values
2778
2779 ((zerop1 arg)
2780 (cond ((eq ($csign (add order 1)) '$zero)
2781 ; http://functions.wolfram.com/03.10.03.0018.01
2782 (cond ((or ($bfloatp order)
2783 ($bfloatp arg))
2784 ($bfloat (div 2 '$%pi)))
2785 ((or (floatp order)
2786 (floatp arg))
2787 ($float (div 2 '$%pi)))
2788 (t
2789 (div 2 '$%pi))))
2790 ((eq ($sign (add ($realpart order) 1)) '$pos)
2791 ; http://functions.wolfram.com/03.10.03.0001.01
2792 arg)
2793 ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz))
2794 (simp-domain-error
2795 (intl:gettext "struve_l: struve_l(~:M,~:M) is undefined.")
2796 order arg))
2797 (t
2798 (eqtest (list '(%struve_l) order arg) expr))))
2799
2800 ;; Check for numerical evaluation
2801
2802 ((complex-float-numerical-eval-p order arg)
2803 ; http://functions.wolfram.com/03.10.26.0002.01
2804 (let (($numer t) ($float t))
2805 ($rectform
2806 ($float
2807 (mul
2808 ($rectform (power arg (add order 1.0)))
2809 ($rectform (inv (power 2.0 order)))
2810 (inv (power ($float '$%pi) 0.5))
2811 (inv (simplify (list '(%gamma) (add order 1.5))))
2812 (simplify (list '($hypergeometric)
2813 (list '(mlist) 1)
2814 (list '(mlist) '((rat simp) 3 2)
2815 (add order '((rat simp) 3 2)))
2816 (div (mul arg arg) 4.0))))))))
2817
2818 ((complex-bigfloat-numerical-eval-p order arg)
2819 ; http://functions.wolfram.com/03.10.26.0002.01
2820 (let (($ratprint nil))
2821 ($rectform
2822 ($bfloat
2823 (mul
2824 ($rectform (power arg (add order 1)))
2825 ($rectform (inv (power 2 order)))
2826 (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2))))
2827 (inv (simplify (list '(%gamma) (add order '((rat simp) 3 2)))))
2828 (simplify (list '($hypergeometric)
2829 (list '(mlist) 1)
2830 (list '(mlist) '((rat simp) 3 2)
2831 (add order '((rat simp) 3 2)))
2832 (div (mul arg arg) 4))))))))
2833
2834 ;; Transformations and argument simplifications
2835
2836 ((and $besselexpand
2837 (ratnump order)
2838 (integerp (mul 2 order)))
2839 (cond
2840 ((eq ($sign order) '$pos)
2841 ;; Expansion of Struve L for a positive half integral order.
2842 (sratsimp
2843 (add
2844 (mul -1
2845 (power 2 (sub 1 order))
2846 (power arg (sub order 1))
2847 (inv (power '$%pi '((rat simp) 1 2)))
2848 (inv (simplify (list '(mfactorial) (sub order
2849 '((rat simp) 1 2)))))
2850 (let ((index (gensumindex)))
2851 (dosum
2852 (mul
2853 (simplify (list '($pochhammer) '((rat simp) 1 2) index))
2854 (simplify (list '($pochhammer)
2855 (sub '((rat simp) 1 2) order)
2856 index))
2857 (power (mul arg arg (inv 4)) (mul -1 index)))
2858 index 0 (sub order '((rat simp) 1 2)) t)))
2859 (mul -1
2860 (power (div 2 '$%pi) '((rat simp) 1 2))
2861 (inv (power arg '((rat simp) 1 2)))
2862 ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2))
2863 (add
2864 (mul
2865 (let (($trigexpand t))
2866 (simplify (list '(%sinh)
2867 (sub (mul '((rat simp) 1 2)
2868 '$%pi
2869 '$%i
2870 (add order '((rat simp) 1 2)))
2871 arg))))
2872 (let ((index (gensumindex)))
2873 (dosum
2874 (mul
2875 (simplify (list '(mfactorial)
2876 (add (mul 2 index)
2877 (simplify (list '(mabs) order))
2878 '((rat simp) -1 2))))
2879 (inv (simplify (list '(mfactorial) (mul 2 index))))
2880 (inv (simplify (list '(mfactorial)
2881 (add (simplify (list '(mabs)
2882 order))
2883 (mul -2 index)
2884 '((rat simp) -1 2)))))
2885 (inv (power (mul 2 arg) (mul 2 index))))
2886 index 0
2887 (simplify (list '($floor)
2888 (div (sub (mul 2
2889 (simplify (list '(mabs)
2890 order)))
2891 1)
2892 4)))
2893 t)))
2894 (mul
2895 (let (($trigexpand t))
2896 (simplify (list '(%cosh)
2897 (sub (mul '((rat simp) 1 2)
2898 '$%pi
2899 '$%i
2900 (add order '((rat simp) 1 2)))
2901 arg))))
2902 (let ((index (gensumindex)))
2903 (dosum
2904 (mul
2905 (simplify (list '(mfactorial)
2906 (add (mul 2 index)
2907 (simplify (list '(mabs) order))
2908 '((rat simp) 1 2))))
2909 (power (mul 2 arg) (neg (add (mul 2 index) 1)))
2910 (inv (simplify (list '(mfactorial)
2911 (add (mul 2 index) 1))))
2912 (inv (simplify (list '(mfactorial)
2913 (add (simplify (list '(mabs)
2914 order))
2915 (mul -2 index)
2916 '((rat simp) -3 2))))))
2917 index 0
2918 (simplify (list '($floor)
2919 (div (sub (mul 2
2920 (simplify (list '(mabs)
2921 order)))
2922 3)
2923 4)))
2924 t))))))))
2925 ((eq ($sign order) '$neg)
2926 ;; Expansion of Struve L for a negative half integral order.
2927 (sratsimp
2928 (add
2929 (mul -1
2930 (power (div 2 '$%pi) '((rat simp) 1 2))
2931 (inv (power arg '((rat simp) 1 2)))
2932 ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2))
2933 (add
2934 (mul
2935 (let (($trigexpand t))
2936 (simplify (list '(%sinh)
2937 (sub (mul '((rat simp) 1 2)
2938 '$%pi
2939 '$%i
2940 (add order '((rat simp) 1 2)))
2941 arg))))
2942 (let ((index (gensumindex)))
2943 (dosum
2944 (mul
2945 (simplify (list '(mfactorial)
2946 (add (mul 2 index)
2947 (neg order)
2948 '((rat simp) -1 2))))
2949 (inv (simplify (list '(mfactorial) (mul 2 index))))
2950 (inv (simplify (list '(mfactorial)
2951 (add (neg order)
2952 (mul -2 index)
2953 '((rat simp) -1 2)))))
2954 (inv (power (mul 2 arg) (mul 2 index))))
2955 index 0
2956 (simplify (list '($floor)
2957 (div (add (mul 2 order) 1) -4)))
2958 t)))
2959 (mul
2960 (let (($trigexpand t))
2961 (simplify (list '(%cosh)
2962 (sub (mul '((rat simp) 1 2)
2963 '$%pi
2964 '$%i
2965 (add order '((rat simp) 1 2)))
2966 arg))))
2967 (let ((index (gensumindex)))
2968 (dosum
2969 (mul
2970 (simplify (list '(mfactorial)
2971 (add (mul 2 index)
2972 (neg order)
2973 '((rat simp) 1 2))))
2974 (power (mul 2 arg) (neg (add (mul 2 index) 1)))
2975 (inv (simplify (list '(mfactorial)
2976 (add (mul 2 index) 1))))
2977 (inv (simplify (list '(mfactorial)
2978 (add (neg order)
2979 (mul -2 index)
2980 '((rat simp) -3 2))))))
2981 index 0
2982 (simplify (list '($floor)
2983 (div (add (mul 2 order) 3) -4)))
2984 t))))))))))
2985 (t
2986 (eqtest (list '(%struve_l) order arg) expr)))))
2987
2988 (def-simplifier struve_l (order arg)
2989 (cond
2990
2991 ;; Check for special values
2992
2993 ((zerop1 arg)
2994 (cond ((eq ($csign (add order 1)) '$zero)
2995 ; http://functions.wolfram.com/03.10.03.0018.01
2996 (cond ((or ($bfloatp order)
2997 ($bfloatp arg))
2998 ($bfloat (div 2 '$%pi)))
2999 ((or (floatp order)
3000 (floatp arg))
3001 ($float (div 2 '$%pi)))
3002 (t
3003 (div 2 '$%pi))))
3004 ((eq ($sign (add ($realpart order) 1)) '$pos)
3005 ; http://functions.wolfram.com/03.10.03.0001.01
3006 arg)
3007 ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz))
3008 (simp-domain-error
3009 (intl:gettext "struve_l: struve_l(~:M,~:M) is undefined.")
3010 order arg))
3011 (t
3012 (give-up))))
3013
3014 ;; Check for numerical evaluation
3015
3016 ((complex-float-numerical-eval-p order arg)
3017 ; http://functions.wolfram.com/03.10.26.0002.01
3018 (let (($numer t) ($float t))
3019 ($rectform
3020 ($float
3021 (mul
3022 ($rectform (power arg (add order 1.0)))
3023 ($rectform (inv (power 2.0 order)))
3024 (inv (power ($float '$%pi) 0.5))
3025 (inv (simplify (list '(%gamma) (add order 1.5))))
3026 (simplify (list '($hypergeometric)
3027 (list '(mlist) 1)
3028 (list '(mlist) '((rat simp) 3 2)
3029 (add order '((rat simp) 3 2)))
3030 (div (mul arg arg) 4.0))))))))
3031
3032 ((complex-bigfloat-numerical-eval-p order arg)
3033 ; http://functions.wolfram.com/03.10.26.0002.01
3034 (let (($ratprint nil))
3035 ($rectform
3036 ($bfloat
3037 (mul
3038 ($rectform (power arg (add order 1)))
3039 ($rectform (inv (power 2 order)))
3040 (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2))))
3041 (inv (simplify (list '(%gamma) (add order '((rat simp) 3 2)))))
3042 (simplify (list '($hypergeometric)
3043 (list '(mlist) 1)
3044 (list '(mlist) '((rat simp) 3 2)
3045 (add order '((rat simp) 3 2)))
3046 (div (mul arg arg) 4))))))))
3047
3048 ;; Transformations and argument simplifications
3049
3050 ((and $besselexpand
3051 (ratnump order)
3052 (integerp (mul 2 order)))
3053 (cond
3054 ((eq ($sign order) '$pos)
3055 ;; Expansion of Struve L for a positive half integral order.
3056 (sratsimp
3057 (add
3058 (mul -1
3059 (power 2 (sub 1 order))
3060 (power arg (sub order 1))
3061 (inv (power '$%pi '((rat simp) 1 2)))
3062 (inv (simplify (list '(mfactorial) (sub order
3063 '((rat simp) 1 2)))))
3064 (let ((index (gensumindex)))
3065 (dosum
3066 (mul
3067 (simplify (list '($pochhammer) '((rat simp) 1 2) index))
3068 (simplify (list '($pochhammer)
3069 (sub '((rat simp) 1 2) order)
3070 index))
3071 (power (mul arg arg (inv 4)) (mul -1 index)))
3072 index 0 (sub order '((rat simp) 1 2)) t)))
3073 (mul -1
3074 (power (div 2 '$%pi) '((rat simp) 1 2))
3075 (inv (power arg '((rat simp) 1 2)))
3076 ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2))
3077 (add
3078 (mul
3079 (let (($trigexpand t))
3080 (simplify (list '(%sinh)
3081 (sub (mul '((rat simp) 1 2)
3082 '$%pi
3083 '$%i
3084 (add order '((rat simp) 1 2)))
3085 arg))))
3086 (let ((index (gensumindex)))
3087 (dosum
3088 (mul
3089 (simplify (list '(mfactorial)
3090 (add (mul 2 index)
3091 (simplify (list '(mabs) order))
3092 '((rat simp) -1 2))))
3093 (inv (simplify (list '(mfactorial) (mul 2 index))))
3094 (inv (simplify (list '(mfactorial)
3095 (add (simplify (list '(mabs)
3096 order))
3097 (mul -2 index)
3098 '((rat simp) -1 2)))))
3099 (inv (power (mul 2 arg) (mul 2 index))))
3100 index 0
3101 (simplify (list '($floor)
3102 (div (sub (mul 2
3103 (simplify (list '(mabs)
3104 order)))
3105 1)
3106 4)))
3107 t)))
3108 (mul
3109 (let (($trigexpand t))
3110 (simplify (list '(%cosh)
3111 (sub (mul '((rat simp) 1 2)
3112 '$%pi
3113 '$%i
3114 (add order '((rat simp) 1 2)))
3115 arg))))
3116 (let ((index (gensumindex)))
3117 (dosum
3118 (mul
3119 (simplify (list '(mfactorial)
3120 (add (mul 2 index)
3121 (simplify (list '(mabs) order))
3122 '((rat simp) 1 2))))
3123 (power (mul 2 arg) (neg (add (mul 2 index) 1)))
3124 (inv (simplify (list '(mfactorial)
3125 (add (mul 2 index) 1))))
3126 (inv (simplify (list '(mfactorial)
3127 (add (simplify (list '(mabs)
3128 order))
3129 (mul -2 index)
3130 '((rat simp) -3 2))))))
3131 index 0
3132 (simplify (list '($floor)
3133 (div (sub (mul 2
3134 (simplify (list '(mabs)
3135 order)))
3136 3)
3137 4)))
3138 t))))))))
3139 ((eq ($sign order) '$neg)
3140 ;; Expansion of Struve L for a negative half integral order.
3141 (sratsimp
3142 (add
3143 (mul -1
3144 (power (div 2 '$%pi) '((rat simp) 1 2))
3145 (inv (power arg '((rat simp) 1 2)))
3146 ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2))
3147 (add
3148 (mul
3149 (let (($trigexpand t))
3150 (simplify (list '(%sinh)
3151 (sub (mul '((rat simp) 1 2)
3152 '$%pi
3153 '$%i
3154 (add order '((rat simp) 1 2)))
3155 arg))))
3156 (let ((index (gensumindex)))
3157 (dosum
3158 (mul
3159 (simplify (list '(mfactorial)
3160 (add (mul 2 index)
3161 (neg order)
3162 '((rat simp) -1 2))))
3163 (inv (simplify (list '(mfactorial) (mul 2 index))))
3164 (inv (simplify (list '(mfactorial)
3165 (add (neg order)
3166 (mul -2 index)
3167 '((rat simp) -1 2)))))
3168 (inv (power (mul 2 arg) (mul 2 index))))
3169 index 0
3170 (simplify (list '($floor)
3171 (div (add (mul 2 order) 1) -4)))
3172 t)))
3173 (mul
3174 (let (($trigexpand t))
3175 (simplify (list '(%cosh)
3176 (sub (mul '((rat simp) 1 2)
3177 '$%pi
3178 '$%i
3179 (add order '((rat simp) 1 2)))
3180 arg))))
3181 (let ((index (gensumindex)))
3182 (dosum
3183 (mul
3184 (simplify (list '(mfactorial)
3185 (add (mul 2 index)
3186 (neg order)
3187 '((rat simp) 1 2))))
3188 (power (mul 2 arg) (neg (add (mul 2 index) 1)))
3189 (inv (simplify (list '(mfactorial)
3190 (add (mul 2 index) 1))))
3191 (inv (simplify (list '(mfactorial)
3192 (add (neg order)
3193 (mul -2 index)
3194 '((rat simp) -3 2))))))
3195 index 0
3196 (simplify (list '($floor)
3197 (div (add (mul 2 order) 3) -4)))
3198 t))))))))))
3199 (t
3200 (give-up))))
3201
3202 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3203
3204 (defmspec $gauss (form)
3205 (format t
3206 "NOTE: The gauss function is superseded by random_normal in the `distrib' package.
3207 Perhaps you meant to enter `~a'.~%"
3208 (print-invert-case (implode (mstring `(($random_normal) ,@ (cdr form))))))
3209 '$done)
Maxima fork w/ changes for Cygwin