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(defprop $bessel_j %bessel_j alias) 56(defprop $bessel_j %bessel_j verb) 57(defprop %bessel_j $bessel_j reversealias) 58(defprop %bessel_j $bessel_j noun) 59 60;; Bessel J is a simplifying function. 61 62(defprop %bessel_j simp-bessel-j operators) 63 64;; Bessel J distributes over lists, matrices, and equations 65 66(defprop %bessel_j (mlist $matrix mequal) distribute_over) 67 68;; Derivatives of the Bessel function. 69(defprop %bessel_j 70 ((n x) 71 ;; Derivative wrt to order n. A&S 9.1.64. Do we really want to 72 ;; do this? It's quite messy. 73 ;; 74 ;; J[n](x)*log(x/2) 75 ;; - (x/2)^n*sum((-1)^k*psi(n+k+1)/gamma(n+k+1)*(x^2/4)^k/k!,k,0,inf) 76 ((mplus) 77 ((mtimes) 78 ((%bessel_j) n x) 79 ((%log) ((mtimes) ((rat) 1 2) x))) 80 ((mtimes) -1 81 ((mexpt) ((mtimes) x ((rat) 1 2)) n) 82 ((%sum) 83 ((mtimes) ((mexpt) -1 $%k) 84 ((mexpt) ((mfactorial) $%k) -1) 85 ((mqapply) (($psi array) 0) ((mplus) 1 $%k n)) 86 ((mexpt) ((%gamma) ((mplus) 1 $%k n)) -1) 87 ((mexpt) ((mtimes) x x ((rat) 1 4)) $%k)) 88 $%k 0 $inf))) 89 90 ;; Derivative wrt to arg x. 91 ;; A&S 9.1.27; changed from 9.1.30 so that taylor works on Bessel functions 92 ((mtimes) 93 ((mplus) 94 ((%bessel_j) ((mplus) -1 n) x) 95 ((mtimes) -1 ((%bessel_j) ((mplus) 1 n) x))) 96 ((rat) 1 2))) 97 grad) 98 99;; Integral of the Bessel function wrt z 100(defun bessel-j-integral-2 (v z) 101 (case v 102 (0 103 ;; integrate(bessel_j(0,z),z) 104 ;; = (1/2)*z*(%pi*bessel_j(1,z)*struve_h(0,z) 105 ;; +bessel_j(0,z)*(2-%pi*struve_h(1,z))) 106 `((mtimes) ((rat) 1 2) ,z 107 ((mplus) 108 ((mtimes) $%pi 109 ((%bessel_j) 1 ,z) 110 ((%struve_h) 0 ,z)) 111 ((mtimes) 112 ((%bessel_j) 0 ,z) 113 ((mplus) 2 ((mtimes) -1 $%pi ((%struve_h) 1 ,z))))))) 114 (1 115 ;; integrate(bessel_j(1,z),z) = -bessel_j(0,z) 116 `((mtimes) -1 ((%bessel_j) 0 ,z))) 117 (otherwise 118 ;; http://functions.wolfram.com/03.01.21.0002.01 119 ;; integrate(bessel_j(v,z),z) 120 ;; = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2) 121 ;; * hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],-z^2/4) 122 ;; = 2^(-v)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],-z^2/4) 123 ;; / gamma(v+2) 124 `((mtimes) 125 (($hypergeometric) 126 ((mlist) 127 ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) ,v))) 128 ((mlist) 129 ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) ,v)) 130 ((mplus) 1 ,v)) 131 ((mtimes) ((rat) -1 4) ((mexpt) ,z 2))) 132 ((mexpt) 2 ((mtimes) -1 ,v)) 133 ((mexpt) ((%gamma) ((mplus) 2 ,v)) -1) 134 ((mexpt) z ((mplus) 1 ,v)))))) 135 136(putprop '%bessel_j `((v z) nil ,#'bessel-j-integral-2) 'integral) 137 138;; Support a simplim%bessel_j function to handle specific limits 139 140(defprop %bessel_j simplim%bessel_j simplim%function) 141 142(defun simplim%bessel_j (expr var val) 143 ;; Look for the limit of the arguments. 144 (let ((v (limit (cadr expr) var val 'think)) 145 (z (limit (caddr expr) var val 'think))) 146 (cond 147 ;; Handle an argument 0 at this place. 148 ((or (zerop1 z) 149 (eq z '$zeroa) 150 (eq z '$zerob)) 151 (let ((sgn ($sign ($realpart v)))) 152 (cond ((and (eq sgn '$neg) 153 (not (maxima-integerp v))) 154 ;; bessel_j(v,0), Re(v)<0 and v not an integer 155 '$infinity) 156 ((and (eq sgn '$zero) 157 (not (zerop1 v))) 158 ;; bessel_j(v,0), Re(v)=0 and v #0 159 '$und) 160 ;; Call the simplifier of the function. 161 (t (simplify (list '(%bessel_j) v z)))))) 162 ((or (eq z '$inf) 163 (eq z '$minf)) 164 ;; bessel_j(v,inf) or bessel_j(v,minf) 165 0) 166 (t 167 ;; All other cases are handled by the simplifier of the function. 168 (simplify (list '(%bessel_j) v z)))))) 169 170(defun simp-bessel-j (expr ignored z) 171 (declare (ignore ignored)) 172 (twoargcheck expr) 173 (let ((order (simpcheck (cadr expr) z)) 174 (arg (simpcheck (caddr expr) z)) 175 (rat-order nil)) 176 (cond 177 ((zerop1 arg) 178 ;; We handle the different case for zero arg carefully. 179 (let ((sgn ($sign ($realpart order)))) 180 (cond ((and (eq sgn '$zero) 181 (zerop1 ($imagpart order))) 182 ;; bessel_j(0,0) = 1 183 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1)) 184 ((or (floatp order) (floatp arg)) 1.0) 185 (t 1))) 186 ((or (eq sgn '$pos) 187 (maxima-integerp order)) 188 ;; bessel_j(v,0) and Re(v)>0 or v an integer 189 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0)) 190 ((or (floatp order) (floatp arg)) 0.0) 191 (t 0))) 192 ((and (eq sgn '$neg) 193 (not (maxima-integerp order))) 194 ;; bessel_j(v,0) and Re(v)<0 and v not an integer 195 (simp-domain-error 196 (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.") 197 order arg)) 198 ((and (eq sgn '$zero) 199 (not (zerop1 ($imagpart order)))) 200 ;; bessel_j(v,0) and Re(v)=0 and v # 0 201 (simp-domain-error 202 (intl:gettext "bessel_j: bessel_j(~:M,~:M) is undefined.") 203 order arg)) 204 (t 205 ;; No information about the sign of the order 206 (eqtest (list '(%bessel_j) order arg) expr))))) 207 208 ((complex-float-numerical-eval-p order arg) 209 ;; We have numeric order and arg and $numer is true, or we have either 210 ;; the order or arg being floating-point, so let's evaluate it 211 ;; numerically. 212 ;; The numerical routine bessel-j returns a CL number, so we have 213 ;; to add the conversion to a Maxima-complex-number. 214 (cond ((= 0 ($imagpart order)) 215 ;; order is real, arg is real or complex 216 (complexify (bessel-j ($float order) 217 (complex ($float ($realpart arg)) 218 ($float ($imagpart arg)))))) 219 (t 220 ;; order is complex, arg is real or complex 221 (let (($numer t) 222 ($float t) 223 (order ($float order)) 224 (arg ($float arg))) 225 ($float 226 ($rectform 227 (bessel-j-hypergeometric order arg))))))) 228 229 ((and (integerp order) (minusp order)) 230 ;; Some special cases when the order is an integer. 231 ;; A&S 9.1.5: J[-n](x) = (-1)^n*J[n](x) 232 (if (evenp order) 233 (take '(%bessel_j) (- order) arg) 234 (mul -1 (take '(%bessel_j) (- order) arg)))) 235 236 ((and $besselexpand 237 (setq rat-order (max-numeric-ratio-p order 2))) 238 ;; When order is a fraction with a denominator of 2, we 239 ;; can express the result in terms of elementary functions. 240 (bessel-j-half-order rat-order arg)) 241 242 ((and $bessel_reduce 243 (and (integerp order) 244 (plusp order) 245 (> order 1))) 246 ;; Reduce a bessel function of order > 2 to order 1 and 0. 247 ;; A&S 9.1.27: bessel_j(v,z) = 2*(v-1)/z*bessel_j(v-1,z)-bessel_j(v-2,z) 248 (sub (mul 2 249 (- order 1) 250 (inv arg) 251 (take '(%bessel_j) (- order 1) arg)) 252 (take '(%bessel_j) (- order 2) arg))) 253 254 ((and $%iargs (multiplep arg '$%i)) 255 ;; bessel_j(v, %i*x) = (%i*x)^v/(x^v) * bessel_i(v, x) 256 ;; (From http://functions.wolfram.com/03.01.27.0002.01) 257 (let ((x (coeff arg '$%i 1))) 258 (mul (power (mul '$%i x) order) 259 (inv (power x order)) 260 (take '(%bessel_i) order x)))) 261 262 ($hypergeometric_representation 263 (bessel-j-hypergeometric order arg)) 264 265 (t 266 (eqtest (list '(%bessel_j) order arg) expr))))) 267 268;; Returns the hypergeometric representation of bessel_j 269(defun bessel-j-hypergeometric (order arg) 270 ;; A&S 9.1.69 / http://functions.wolfram.com/03.01.26.0002.01 271 ;; hypergeometric([],[v+1],-z^2/4)*(z/2)^v/gamma(v+1) 272 (mul (inv (take '(%gamma) (add order 1))) 273 (power (div arg 2) order) 274 (take '($hypergeometric) 275 (list '(mlist)) 276 (list '(mlist) (add order 1)) 277 (neg (div (mul arg arg) 4))))) 278 279;; Compute value of Bessel function of the first kind of order ORDER. 280(defun bessel-j (order arg) 281 (cond 282 ((zerop (imagpart arg)) 283 ;; We have numeric args and the arg is purely real. 284 ;; Call the real-valued Bessel function when possible. 285 (let ((arg (realpart arg))) 286 (cond 287 ((= order 0) 288 (slatec:dbesj0 (float arg))) 289 ((= order 1) 290 (slatec:dbesj1 (float arg))) 291 ((minusp order) 292 (cond ((zerop (nth-value 1 (truncate order))) 293 ;; The order is a negative integer. We use A&S 9.1.5: 294 ;; J[-n](z) = (-1)^n*J[n](z) 295 ;; and not the Hankel functions. 296 (if (evenp (floor order)) 297 (bessel-j (- order) arg) 298 (- (bessel-j (- order) arg)))) 299 (t 300 ;; Bessel function of negative order. We use the Hankel 301 ;; functions to compute this: 302 ;; J[v](z) = (H1[v](z) + H2[v](z)) / 2. 303 ;; This works for negative and positive arg and handles 304 ;; special cases correctly. 305 (let ((result (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg))))) 306 (cond ((= (nth-value 1 (floor order)) 1/2) 307 ;; ORDER is a half-integral-value or a float 308 ;; representation, thereof. 309 (if (minusp arg) 310 ;; arg is negative, the result is purely imaginary 311 (complex 0 (imagpart result)) 312 ;; arg is positive, the result is purely real 313 (realpart result))) 314 ;; in all other cases general complex result 315 (t result)))))) 316 (t 317 ;; We have a real arg and order > 0 and order not 0 or 1 318 ;; for this case we can call the function dbesj 319 (multiple-value-bind (n alpha) (floor (float order)) 320 (let ((jvals (make-array (1+ n) :element-type 'flonum))) 321 (slatec:dbesj (abs (float arg)) alpha (1+ n) jvals 0) 322 (cond ((>= arg 0) 323 (aref jvals n)) 324 (t 325 ;; Use analytic continuation formula A&S 9.1.35: 326 ;; J[v](z*exp(m*%pi*%i)) = exp(m*%pi*%i*v)*J[v](z) 327 ;; for an integer m. In particular, for m = 1: 328 ;; J[v](-x) = exp(v*%pi*%i)*J[v](x) 329 ;; and handle special cases 330 (cond ((zerop (nth-value 1 (truncate order))) 331 ;; order is an integer 332 (if (evenp (floor order)) 333 (aref jvals n) 334 (- (aref jvals n)))) 335 ((= (nth-value 1 (floor order)) 1/2) 336 ;; Order is a half-integral-value and we know that 337 ;; arg < 0, so the result is purely imginary. 338 (if (evenp (floor order)) 339 (complex 0 (aref jvals n)) 340 (complex 0 (- (aref jvals n))))) 341 ;; In all other cases a general complex result 342 (t 343 (* (cis (* order pi)) 344 (aref jvals n)))))))))))) 345 (t 346 ;; The arg is complex. Use the complex-valued Bessel function. 347 (cond ((mminusp order) 348 ;; Bessel function of negative order. We use the Hankel function to 349 ;; compute this, because A&S 9.1.3 and 9.1.4 say 350 ;; H1[v](z) = J[v](z) + %i * Y[v](z) 351 ;; and 352 ;; H2[v](z) = J[v](z) - %i * Y[v](z). 353 ;; Thus 354 ;; J[v](z) = (H1[v](z) + H2[v](z)) / 2. 355 ;; Not the most efficient way, but perhaps good enough for maxima. 356 (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg)))) 357 (t 358 (multiple-value-bind (n alpha) (floor (float order)) 359 (let ((cyr (make-array (1+ n) :element-type 'flonum)) 360 (cyi (make-array (1+ n) :element-type 'flonum))) 361 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n 362 v-cyr v-cyi v-nz v-ierr) 363 (slatec:zbesj (float (realpart arg)) 364 (float (imagpart arg)) 365 alpha 1 (1+ n) cyr cyi 0 0) 366 (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz)) 367 368 ;; Should check the return status in v-ierr of this routine. 369 (when (plusp v-ierr) 370 (format t "zbesj ierr = ~A~%" v-ierr)) 371 (complex (aref cyr n) (aref cyi n)))))))))) 372 373;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 374;;; 375;;; Implementation of the Bessel Y function 376;;; 377;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 378 379(defmfun $bessel_y (v z) 380 (simplify (list '(%bessel_y) v z))) 381 382(defprop $bessel_y %bessel_y alias) 383(defprop $bessel_y %bessel_y verb) 384(defprop %bessel_y $bessel_y reversealias) 385(defprop %bessel_y $bessel_y noun) 386 387(defprop %bessel_y simp-bessel-y operators) 388 389;; Bessel Y distributes over lists, matrices, and equations 390 391(defprop %bessel_y (mlist $matrix mequal) distribute_over) 392 393(defprop %bessel_y 394 ((n x) 395 ;; Derivative wrt order n. A&S 9.1.65. 396 ;; 397 ;; cot(n*%pi)*[diff(bessel_j(n,x),n)-%pi*bessel_y(n,x)] 398 ;; - csc(n*%pi)*diff(bessel_j(-n,x),n)-%pi*bessel_j(n,x) 399 ((mplus) 400 ((mtimes) -1 $%pi ((%bessel_j) n x)) 401 ((mtimes) 402 -1 403 ((%csc) ((mtimes) $%pi n)) 404 ((%derivative) ((%bessel_j) ((mtimes) -1 n) x) n 1)) 405 ((mtimes) 406 ((%cot) ((mtimes) $%pi n)) 407 ((mplus) 408 ((mtimes) -1 $%pi ((%bessel_y) n x)) 409 ((%derivative) ((%bessel_j) n x) n 1)))) 410 411 ;; Derivative wrt to arg x. A&S 9.1.27; changed from A&S 9.1.30 412 ;; to be consistent with bessel_j. 413 ((mtimes) 414 ((mplus) 415 ((%bessel_y)((mplus) -1 n) x) 416 ((mtimes) -1 ((%bessel_y) ((mplus) 1 n) x))) 417 ((rat) 1 2))) 418 grad) 419 420;; Integral of the Bessel Y function wrt z 421;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/21/01/01/ 422(defun bessel-y-integral-2 (n z) 423 (cond 424 ((and ($integerp n) (<= 0 n)) 425 (cond 426 (($oddp n) 427 ;; integrate(bessel_y(2*N+1,z),z) , N > 0 428 ;; = -bessel_y(0,z) - 2 * sum(bessel_y(2*k,z),k,1,N) 429 (let* ((k (gensym)) 430 (answer `((mplus) ((mtimes) -1 ((%bessel_y) 0 ,z)) 431 ((mtimes) -2 432 ((%sum) ((%bessel_y) ((mtimes) 2 ,k) ,z) ,k 1 433 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n))))))) 434 ;; Expand out the sum if n < 10. Otherwise fix up the indices 435 (if (< n 10) 436 (meval `(($ev) ,answer $sum)) ; Is there a better way? 437 (simplify ($niceindices answer))))) 438 (($evenp n) 439 ;; integrate(bessel_y(2*N,z),z) , N > 0 440 ;; = (1/2)*%pi*z*(bessel_y(0,z)*struve_h(-1,z) 441 ;; +bessel_y(1,z)*struve_h(0,z)) 442 ;; - 2 * sum(bessel_y(2*k+1,z),k,0,N-1) 443 (let* 444 ((k (gensym)) 445 (answer `((mplus) 446 ((mtimes) -2 447 ((%sum) 448 ((%bessel_y) ((mplus) 1 ((mtimes) 2 ,k)) ,z) 449 ,k 0 450 ((mplus) 451 -1 452 ((mtimes) ((rat) 1 2) ,n)))) 453 ((mtimes) ((rat) 1 2) $%pi ,z 454 ((mplus) 455 ((mtimes) 456 ((%bessel_y) 0 ,z) 457 ((%struve_h) -1 ,z)) 458 ((mtimes) 459 ((%bessel_y) 1 ,z) 460 ((%struve_h) 0 ,z))))))) 461 ;; Expand out the sum if n < 10. Otherwise fix up the indices 462 (if (< n 10) 463 (meval `(($ev) ,answer $sum)) ; Is there a better way? 464 (simplify ($niceindices answer))))))) 465 (t nil))) 466 467(putprop '%bessel_y `((n z) nil ,#'bessel-y-integral-2) 'integral) 468 469;; Support a simplim%bessel_y function to handle specific limits 470 471(defprop %bessel_y simplim%bessel_y simplim%function) 472 473(defun simplim%bessel_y (expr var val) 474 ;; Look for the limit of the arguments. 475 (let ((v (limit (cadr expr) var val 'think)) 476 (z (limit (caddr expr) var val 'think))) 477 (cond 478 ;; Handle an argument 0 at this place. 479 ((or (zerop1 z) 480 (eq z '$zeroa) 481 (eq z '$zerob)) 482 (cond ((zerop1 v) 483 ;; bessel_y(0,0) 484 '$minf) 485 ((integerp v) 486 ;; bessel_y(n,0), n an integer 487 (cond ((evenp v) '$minf) 488 (t (cond ((eq z '$zeroa) '$minf) 489 ((eq z '$zerob) '$inf) 490 (t '$infinity))))) 491 ((not (zerop1 ($realpart v))) 492 ;; bessel_y(v,0), Re(v)#0 493 '$infinity) 494 ((and (zerop1 ($realpart v)) 495 (not (zerop1 v))) 496 ;; bessel_y(v,0), Re(v)=0 and v#0 497 '$und) 498 ;; Call the simplifier of the function. 499 (t (simplify (list '(%bessel_y) v z))))) 500 ((or (eq z '$inf) 501 (eq z '$minf)) 502 ;; bessel_y(v,inf) or bessel_y(v,minf) 503 0) 504 (t 505 ;; All other cases are handled by the simplifier of the function. 506 (simplify (list '(%bessel_y) v z)))))) 507 508(defun simp-bessel-y (expr ignored z) 509 (declare (ignore ignored)) 510 (twoargcheck expr) 511 (let ((order (simpcheck (cadr expr) z)) 512 (arg (simpcheck (caddr expr) z)) 513 (rat-order nil)) 514 (cond 515 ((zerop1 arg) 516 ;; Domain error for a zero argument. 517 (simp-domain-error 518 (intl:gettext "bessel_y: bessel_y(~:M,~:M) is undefined.") order arg)) 519 520 ((complex-float-numerical-eval-p order arg) 521 ;; We have numeric order and arg and $numer is true, or 522 ;; we have either the order or arg being floating-point, 523 ;; so let's evaluate it numerically. 524 (cond ((= 0 ($imagpart order)) 525 ;; order is real, arg is real or complex 526 (complexify (bessel-y ($float order) 527 (complex ($float ($realpart arg)) 528 ($float ($imagpart arg)))))) 529 (t 530 ;; order is complex, arg is real or complex 531 (let (($numer t) 532 ($float t) 533 (order ($float order)) 534 (arg ($float arg))) 535 ($float 536 ($rectform 537 (bessel-y-hypergeometric order arg))))))) 538 539 ((and (integerp order) (minusp order)) 540 ;; Special case when the order is an integer. 541 ;; A&S 9.1.5: Y[-n](x) = (-1)^n*Y[n](x) 542 (if (evenp order) 543 (take '(%bessel_y) (- order) arg) 544 (mul -1 (take '(%bessel_y) (- order) arg)))) 545 546 ((and $besselexpand 547 (setq rat-order (max-numeric-ratio-p order 2))) 548 ;; When order is a fraction with a denominator of 2, we 549 ;; can express the result in terms of elementary functions. 550 ;; From A&S 10.1.1, 10.1.11-12 and 10.1.15: 551 ;; 552 ;; Y[1/2](z) = -J[-1/2](z) is a function of cos(z). 553 ;; Y[-1/2](z) = J[1/2](z) is a function of sin(z). 554 (bessel-y-half-order rat-order arg)) 555 556 ((and $bessel_reduce 557 (and (integerp order) 558 (plusp order) 559 (> order 1))) 560 ;; Reduce a bessel function of order > 2 to order 1 and 0. 561 ;; A&S 9.1.27: bessel_y(v,z) = 2*(v-1)/z*bessel_y(v-1,z)-bessel_y(v-2,z) 562 (sub (mul 2 563 (- order 1) 564 (inv arg) 565 (take '(%bessel_y) (- order 1) arg)) 566 (take '(%bessel_y) (- order 2) arg))) 567 568 ($hypergeometric_representation 569 (bessel-y-hypergeometric order arg)) 570 571 (t 572 (eqtest (list '(%bessel_y) order arg) expr))))) 573 574;; Returns the hypergeometric representation of bessel_y 575(defun bessel-y-hypergeometric (order arg) 576 ;; http://functions.wolfram.com/03.03.26.0002.01 577 ;; -hypergeometric([],[1-v],-z^2/4)*(2/z)^v*gamma(v)/%pi 578 ;; - hypergeometric([],[v+1],-z^2/4)*gamma(-v)*cos(%pi*v)*(z/2)^v/%pi 579 (add (mul -1 580 (power 2 order) 581 (inv (power arg order)) 582 (take '(%gamma) order) 583 (inv '$%pi) 584 (take '($hypergeometric) 585 (list '(mlist)) 586 (list '(mlist) (sub 1 order)) 587 (neg (div (mul arg arg) 4)))) 588 (mul -1 589 (inv (power 2 order)) 590 (power arg order) 591 (take '(%cos) (mul order '$%pi)) 592 (take '(%gamma) (neg order)) 593 (inv '$%pi) 594 (take '($hypergeometric) 595 (list '(mlist)) 596 (list '(mlist) (add 1 order)) 597 (neg (div (mul arg arg) 4)))))) 598 599;; Bessel function of the second kind, Y[n](z), for real or complex z 600(defun bessel-y (order arg) 601 (cond 602 ((zerop (imagpart arg)) 603 ;; We have numeric args and the first arg is purely 604 ;; real. Call the real-valued Bessel function. 605 ;; 606 ;; For negative values, use the analytic continuation formula 607 ;; A&S 9.1.36: 608 ;; 609 ;; Y[v](z*exp(m*%pi*%i)) = exp(-v*m*%pi*%i) * Y[v](z) 610 ;; + 2*%i*sin(m*v*%pi)*cot(v*%pi)*J[v](z) 611 ;; 612 ;; In particular for m = 1: 613 ;; Y[v](-z) = exp(-v*%pi*%i) * Y[v](z) + 2*%i*cos(v*%pi)*J[v](z) 614 (let ((arg (realpart arg))) 615 (cond 616 ((zerop order) 617 (cond ((>= arg 0) 618 (slatec:dbesy0 (float arg))) 619 (t 620 ;; For v = 0, this simplifies to 621 ;; Y[0](-z) = Y[0](z) + 2*%i*J[0](z) 622 ;; the return value has to be a CL number 623 (+ (slatec:dbesy0 (float (- arg))) 624 (complex 0 (* 2 (slatec:dbesj0 (float (- arg))))))))) 625 ((= order 1) 626 (cond ((>= arg 0) 627 (slatec:dbesy1 (float arg))) 628 (t 629 ;; For v = 1, this simplifies to 630 ;; Y[1](-z) = -Y[1](z) - 2*%i*J[1](v) 631 ;; the return value has to be a CL number 632 (+ (- (slatec:dbesy1 (float (- arg)))) 633 (complex 0 (* -2 (slatec:dbesj1 (float (- arg))))))))) 634 ((minusp order) 635 (cond ((zerop (nth-value 1 (truncate order))) 636 ;; Order is a negative integer or float representation. 637 ;; Use A&S 9.1.5: Y[-n](z) = (-1)^n*Y[n](z). 638 (if (evenp (floor order)) 639 (bessel-y (- order) arg) 640 (- (bessel-y (- order) arg)))) 641 (t 642 ;; Bessel function of negative order. We use the Hankel 643 ;; functions to compute this: 644 ;; Y[v](z) = (H1[v](z) - H2[v](z)) / 2 / %i 645 (let ((result (/ (- (hankel-1 order arg) 646 (hankel-2 order arg)) 647 (complex 0 2)))) 648 (cond ((= (nth-value 1 (floor order)) 1/2) 649 ;; ORDER is half-integral-value or a float 650 ;; representation thereof. 651 (if (minusp arg) 652 ;; arg is negative, the result is purely imaginary 653 (complex 0 (imagpart result)) 654 ;; arg is positive, the result is purely real 655 (realpart result))) 656 ;; in all other cases general complex result 657 (t result)))))) 658 (t 659 (multiple-value-bind (n alpha) (floor (float order)) 660 (let ((jvals (make-array (1+ n) :element-type 'flonum))) 661 ;; First we do the calculation for an positive argument. 662 (slatec:dbesy (abs (float arg)) alpha (1+ n) jvals) 663 664 ;; Now we look at the sign of the argument 665 (cond ((>= arg 0) 666 (aref jvals n)) 667 (t 668 (let* ((dpi (coerce pi 'flonum)) 669 (s1 (cis (- (* order dpi)))) 670 (s2 (* #c(0 2) (cos (* order dpi))))) 671 (let ((result (+ (* s1 (aref jvals n)) 672 (* s2 (bessel-j order (- arg)))))) 673 (cond ((zerop (nth-value 1 (truncate order))) 674 ;; ORDER is an integer or a float representation 675 ;; of an integer, and the arg is positive the 676 ;; result is general complex. 677 result) 678 ((= (nth-value 1 (floor order)) 1/2) 679 ;; ORDER is a half-integral-value or an float 680 ;; representation and we have arg < 0. The 681 ;; result is purely imaginary. 682 (complex 0 (imagpart result))) 683 ;; in all other cases general complex result 684 (t result)))))))))))) 685 (t 686 (cond ((minusp order) 687 ;; Bessel function of negative order. We use the Hankel function to 688 ;; compute this, because A&S 9.1.3 and 9.1.4 say 689 ;; H1[v](z) = J[v](z) + %i * Y[v](z) 690 ;; and 691 ;; H2[v](z) = J[v](z) - %i * Y[v](z). 692 ;; Thus 693 ;; Y[v](z) = (H1[v](z) - H2[v](z)) / 2 / %i 694 (/ (- (hankel-1 order arg) (hankel-2 order arg)) 695 (complex 0 2))) 696 (t 697 (multiple-value-bind (n alpha) (floor (float order)) 698 (let ((cyr (make-array (1+ n) :element-type 'flonum)) 699 (cyi (make-array (1+ n) :element-type 'flonum)) 700 (cwrkr (make-array (1+ n) :element-type 'flonum)) 701 (cwrki (make-array (1+ n) :element-type 'flonum))) 702 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi 703 v-nz v-cwrkr v-cwrki v-ierr) 704 (slatec::zbesy (float (realpart arg)) 705 (float (imagpart arg)) 706 alpha 1 (1+ n) cyr cyi 0 cwrkr cwrki 0) 707 (declare (ignore v-zr v-zi v-fnu v-kode v-n 708 v-cyr v-cyi v-cwrkr v-cwrki v-nz)) 709 710 ;; We should check for errors based on the value of v-ierr. 711 (when (plusp v-ierr) 712 (format t "zbesy ierr = ~A~%" v-ierr)) 713 714 (complex (aref cyr n) (aref cyi n)))))))))) 715 716;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 717;;; 718;;; Implementation of the Bessel I function 719;;; 720;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 721 722(defmfun $bessel_i (v z) 723 (simplify (list '(%bessel_i) v z))) 724 725(defprop $bessel_i %bessel_i alias) 726(defprop $bessel_i %bessel_i verb) 727(defprop %bessel_i $bessel_i reversealias) 728(defprop %bessel_i $bessel_i noun) 729 730(defprop %bessel_i simp-bessel-i operators) 731 732;; Bessel I distributes over lists, matrices, and equations 733 734(defprop %bessel_i (mlist $matrix mequal) distribute_over) 735 736(defprop %bessel_i 737 ((n x) 738 ;; Derivative wrt order n. A&S 9.6.42. 739 ;; 740 ;; I[n](x)*log(x/2) 741 ;; - (x/2)^n*sum(psi(n+k+1)/gamma(n+k+1)*(x^2/4)^k/k!,k,0,inf) 742 ((mplus) 743 ((mtimes) 744 ((%bessel_i) n x) 745 ((%log) ((mtimes) ((rat) 1 2) x))) 746 ((mtimes) -1 747 ((mexpt) ((mtimes) x ((rat) 1 2)) n) 748 ((%sum) 749 ((mtimes) 750 ((mexpt) ((mfactorial) $%k) -1) 751 ((mqapply) (($psi array) 0) ((mplus) 1 $%k n)) 752 ((mexpt) ((%gamma) ((mplus) 1 $%k n)) -1) 753 ((mexpt) ((mtimes) x x ((rat) 1 4)) $%k)) 754 $%k 0 $inf))) 755 ;; Derivative wrt to x. A&S 9.6.29. 756 ((mtimes) 757 ((mplus) ((%bessel_i) ((mplus) -1 n) x) 758 ((%bessel_i) ((mplus) 1 n) x)) 759 ((rat) 1 2))) 760 grad) 761 762;; Integral of the Bessel I function wrt z 763;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselI/21/01/01/ 764(defun bessel-i-integral-2 (v z) 765 (case v 766 (0 767 ;; integrate(bessel_i(0,z),z) 768 ;; = (1/2)*z*(bessel_i(0,z)*(%pi*struve_l(1,z)+2) 769 ;; -%pi*bessel_i(1,z)*struve_l(0,z)) 770 `((mtimes) ((rat) 1 2) ,z 771 ((mplus) 772 ((mtimes) -1 $%pi 773 ((%bessel_i) 1 ,z) 774 ((%struve_l) 0 ,z)) 775 ((mtimes) 776 ((%bessel_i) 0 ,z) 777 ((mplus) 2 778 ((mtimes) $%pi ((%struve_l) 1 ,z))))))) 779 (1 780 ;; integrate(bessel_i(1,z),z) = bessel_i(0,z) 781 `((%bessel_i) 0 ,z)) 782 (otherwise 783 ;; http://functions.wolfram.com/03.02.21.0002.01 784 ;; integrate(bessel_i(v,z),z) 785 ;; = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2) 786 ;; * hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],z^2/4) 787 ;; = 2^(-v)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],z^2/4) 788 ;; / gamma(v+2) 789 `((mtimes) 790 (($hypergeometric) 791 ((mlist) 792 ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) ,v))) 793 ((mlist) 794 ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) ,v)) 795 ((mplus) 1 ,v)) 796 ((mtimes) ((rat) 1 4) ((mexpt) ,z 2))) 797 ((mexpt) 2 ((mtimes) -1 ,v)) 798 ((mexpt) ((%gamma) ((mplus) 2 ,v)) -1) 799 ((mexpt) z ((mplus) 1 ,v)))))) 800 801(putprop '%bessel_i `((v z) nil ,#'bessel-i-integral-2) 'integral) 802 803;; Support a simplim%bessel_i function to handle specific limits 804 805(defprop %bessel_i simplim%bessel_i simplim%function) 806 807(defun simplim%bessel_i (expr var val) 808 ;; Look for the limit of the arguments. 809 (let ((v (limit (cadr expr) var val 'think)) 810 (z (limit (caddr expr) var val 'think))) 811 (cond 812 ;; Handle an argument 0 at this place. 813 ((or (zerop1 z) 814 (eq z '$zeroa) 815 (eq z '$zerob)) 816 (let ((sgn ($sign ($realpart v)))) 817 (cond ((and (eq sgn '$neg) 818 (not (maxima-integerp v))) 819 ;; bessel_i(v,0), Re(v)<0 and v not an integer 820 '$infinity) 821 ((and (eq sgn '$zero) 822 (not (zerop1 v))) 823 ;; bessel_i(v,0), Re(v)=0 and v #0 824 '$und) 825 ;; Call the simplifier of the function. 826 (t (simplify (list '(%bessel_i) v z)))))) 827 ((eq z '$inf) 828 ;; bessel_i(v,inf) 829 '$inf) 830 ((eq z '$minf) 831 ;; bessel_i(v,minf) 832 '$infinity) 833 (t 834 ;; All other cases are handled by the simplifier of the function. 835 (simplify (list '(%bessel_i) v z)))))) 836 837(defun simp-bessel-i (expr ignored z) 838 (declare (ignore ignored)) 839 (twoargcheck expr) 840 (let ((order (simpcheck (cadr expr) z)) 841 (arg (simpcheck (caddr expr) z)) 842 (rat-order nil)) 843 (cond 844 ((zerop1 arg) 845 ;; We handle the different case for zero arg carefully. 846 (let ((sgn ($sign ($realpart order)))) 847 (cond ((and (eq sgn '$zero) 848 (zerop1 ($imagpart order))) 849 ;; bessel_i(0,0) = 1 850 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 1)) 851 ((or (floatp order) (floatp arg)) 1.0) 852 (t 1))) 853 ((or (eq sgn '$pos) 854 (maxima-integerp order)) 855 ;; bessel_i(v,0) and Re(v)>0 or v an integer 856 (cond ((or ($bfloatp order) ($bfloatp arg)) ($bfloat 0)) 857 ((or (floatp order) (floatp arg)) 0.0) 858 (t 0))) 859 ((and (eq sgn '$neg) 860 (not (maxima-integerp order))) 861 ;; bessel_i(v,0) and Re(v)<0 and v not an integer 862 (simp-domain-error 863 (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.") 864 order arg)) 865 ((and (eq sgn '$zero) 866 (not (zerop1 ($imagpart order)))) 867 ;; bessel_i(v,0) and Re(v)=0 and v # 0 868 (simp-domain-error 869 (intl:gettext "bessel_i: bessel_i(~:M,~:M) is undefined.") 870 order arg)) 871 (t 872 ;; No information about the sign of the order 873 (eqtest (list '(%bessel_i) order arg) expr))))) 874 875 ((complex-float-numerical-eval-p order arg) 876 (cond ((= 0 ($imagpart order)) 877 ;; order is real, arg is real or complex 878 (complexify (bessel-i ($float order) 879 (complex ($float ($realpart arg)) 880 ($float ($imagpart arg)))))) 881 (t 882 ;; order is complex, arg is real or complex 883 (let (($numer t) 884 ($float t) 885 (order ($float order)) 886 (arg ($float arg))) 887 ($float 888 ($rectform 889 (bessel-i-hypergeometric order arg))))))) 890 891 ((and (integerp order) (minusp order)) 892 ;; Some special cases when the order is an integer 893 ;; A&S 9.6.6: I[-n](x) = I[n](x) 894 (take '(%bessel_i) (- order) arg)) 895 896 ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2))) 897 ;; When order is a fraction with a denominator of 2, we 898 ;; can express the result in terms of elementary functions. 899 ;; From A&S 10.2.13 and 10.2.14: 900 ;; 901 ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z) 902 ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z) 903 (bessel-i-half-order rat-order arg)) 904 905 ((and $bessel_reduce 906 (and (integerp order) 907 (plusp order) 908 (> order 1))) 909 ;; Reduce a bessel function of order > 2 to order 1 and 0. 910 ;; A&S 9.6.26: bessel_i(v,z) = -2*(v-1)/z*bessel_i(v-1,z)+bessel_i(v-2,z) 911 (add (mul -2 912 (- order 1) 913 (inv arg) 914 (take '(%bessel_i) (- order 1) arg)) 915 (take '(%bessel_i) (- order 2) arg))) 916 917 ((and $%iargs (multiplep arg '$%i)) 918 ;; bessel_i(v, %i*x) = (%i*x)^v/(x^v) * bessel_j(v, x) 919 ;; (From http://functions.wolfram.com/03.02.27.0002.01) 920 (let ((x (coeff arg '$%i 1))) 921 (mul (power (mul '$%i x) order) 922 (inv (power x order)) 923 (take '(%bessel_j) order x)))) 924 925 ($hypergeometric_representation 926 (bessel-i-hypergeometric order arg)) 927 928 (t 929 (eqtest (list '(%bessel_i) order arg) expr))))) 930 931;; Returns the hypergeometric representation of bessel_i 932(defun bessel-i-hypergeometric (order arg) 933 ;; A&S 9.6.47 / http://functions.wolfram.com/03.02.26.0002.01 934 ;; hypergeometric([],[v+1],z^2/4)*(z/2)^v/gamma(v+1) 935 (mul (inv (take '(%gamma) (add order 1))) 936 (power (div arg 2) order) 937 (take '($hypergeometric) 938 (list '(mlist)) 939 (list '(mlist) (add order 1)) 940 (div (mul arg arg) 4)))) 941 942;; Compute value of Modified Bessel function of the first kind of order ORDER 943(defun bessel-i (order arg) 944 (cond 945 ((zerop (imagpart arg)) 946 ;; We have numeric args and the first arg is purely 947 ;; real. Call the real-valued Bessel function. Use special 948 ;; routines for order 0 and 1, when possible 949 (let ((arg (realpart arg))) 950 (cond 951 ((zerop order) 952 (slatec:dbesi0 (float arg))) 953 ((= order 1) 954 (slatec:dbesi1 (float arg))) 955 ((or (minusp order) (< arg 0)) 956 (multiple-value-bind (order-int order-frac) (floor order) 957 (cond ((zerop order-frac) 958 ;; Order is an integer. From A&S 9.6.6 and 9.6.30 we have 959 ;; I[-n](z) = I[n](z) 960 ;; and 961 ;; I[n](-z) = (-1)^n*I[n](z) 962 (if (< arg 0) 963 (if (evenp order-int) 964 (bessel-i (abs order) (abs arg)) 965 (- (bessel-i (abs order) (abs arg)))) 966 (bessel-i (abs order) arg))) 967 (t 968 ;; Order or arg is negative and order is not an integer, use 969 ;; the bessel-j function for calculation. We know from 970 ;; the definition I[v](x) = x^v/(%i*x)^v * J[v](%i*x). 971 (let* ((arg (float arg)) 972 (result (* (expt arg order) 973 (expt (complex 0 arg) (- order)) 974 (bessel-j order (complex 0 arg))))) 975 ;; Try to clean up result if we know the result is 976 ;; purely real or purely imaginary. 977 (cond ((>= arg 0) 978 ;; Result is purely real for arg >= 0 979 (realpart result)) 980 ((zerop order-frac) 981 ;; Order is an integer or a float representation of 982 ;; an integer, the result is purely real. 983 (realpart result)) 984 ((= order-frac 1/2) 985 ;; Order is half-integral-value or a float 986 ;; representation and arg < 0, the result 987 ;; is purely imaginary. 988 (complex 0 (imagpart result))) 989 (t result))))))) 990 (t 991 ;; Now the case order > 0 and arg >= 0 992 (multiple-value-bind (n alpha) (floor (float order)) 993 (let ((jvals (make-array (1+ n) :element-type 'flonum))) 994 (slatec:dbesi (float (realpart arg)) alpha 1 (1+ n) jvals 0) 995 (aref jvals n))))))) 996 997 ((and (zerop (realpart arg)) 998 (zerop (rem order 1))) 999 ;; Handle the case for a pure imaginary arg and integer order. 1000 ;; In this case, the result is purely real or purely imaginary, 1001 ;; so we want to make sure that happens. 1002 ;; 1003 ;; bessel_i(n, %i*x) = (%i*x)^n/x^n * bessel_j(n,x) 1004 ;; = %i^n * bessel_j(n,x) 1005 (let* ((n (floor order)) 1006 (result (bessel-j (float n) (imagpart arg)))) 1007 (cond ((evenp n) 1008 ;; %i^(2*m) = (-1)^m, where n=2*m 1009 (if (evenp (/ n 2)) 1010 result 1011 (- result))) 1012 ((oddp n) 1013 ;; %i^(2*m+1) = %i*(-1)^m, where n = 2*m+1 1014 (if (evenp (floor n 2)) 1015 (complex 0 result) 1016 (complex 0 (- result))))))) 1017 1018 (t 1019 ;; The arg is complex. Use the complex-valued Bessel function. 1020 (multiple-value-bind (n alpha) (floor (abs (float order))) 1021 ;; Evaluate the function for positive order and fixup the result later. 1022 (let ((cyr (make-array (1+ n) :element-type 'flonum)) 1023 (cyi (make-array (1+ n) :element-type 'flonum))) 1024 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n 1025 v-cyr v-cyi v-nz v-ierr) 1026 (slatec::zbesi (float (realpart arg)) 1027 (float (imagpart arg)) 1028 alpha 1 (1+ n) cyr cyi 0 0) 1029 (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz)) 1030 ;; We should check for errors here based on the value of v-ierr. 1031 (when (plusp v-ierr) 1032 (format t "zbesi ierr = ~A~%" v-ierr)) 1033 1034 ;; We have evaluated I[|order|](arg), now we look at 1035 ;; the the sign of the order. 1036 (cond ((minusp order) 1037 ;; From A&S 9.6.2: 1038 ;; I[-a](z) = I[a](z) + (2/%pi)*sin(%pi*a)*K[a](z) 1039 (+ (complex (aref cyr n) (aref cyi n)) 1040 (let ((dpi (coerce pi 'flonum))) 1041 (* (/ 2.0 dpi) 1042 (sin (* dpi (- order))) 1043 (bessel-k (- order) arg))))) 1044 (t 1045 (complex (aref cyr n) (aref cyi n)))))))))) 1046 1047;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1048;;; 1049;;; Implementation of the Bessel K function 1050;;; 1051;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1052 1053(defmfun $bessel_k (v z) 1054 (simplify (list '(%bessel_k) v z))) 1055 1056(defprop $bessel_k %bessel_k alias) 1057(defprop $bessel_k %bessel_k verb) 1058(defprop %bessel_k $bessel_k reversealias) 1059(defprop %bessel_k $bessel_k noun) 1060 1061(defprop %bessel_k simp-bessel-k operators) 1062 1063;; Bessel K distributes over lists, matrices, and equations 1064 1065(defprop %bessel_k (mlist $matrix mequal) distribute_over) 1066 1067(defprop %bessel_k 1068 ((n x) 1069 ;; Derivativer wrt order n. A&S 9.6.43. 1070 ;; 1071 ;; %pi/2*csc(n*%pi)*['diff(bessel_i(-n,x),n)-'diff(bessel_i(n,x),n)] 1072 ;; - %pi*cot(n*%pi)*bessel_k(n,x) 1073 ((mplus) 1074 ((mtimes) -1 $%pi 1075 ((%bessel_k) n x) 1076 ((%cot) ((mtimes) $%pi n))) 1077 ((mtimes) 1078 ((rat) 1 2) 1079 $%pi 1080 ((%csc) ((mtimes) $%pi n)) 1081 ((mplus) 1082 ((%derivative) ((%bessel_i) ((mtimes) -1 n) x) n 1) 1083 ((mtimes) -1 1084 ((%derivative) ((%bessel_i) n x) n 1))))) 1085 ;; Derivative wrt to x. A&S 9.6.29. 1086 ((mtimes) 1087 -1 1088 ((mplus) ((%bessel_k) ((mplus) -1 n) x) 1089 ((%bessel_k) ((mplus) 1 n) x)) 1090 ((rat) 1 2))) 1091 grad) 1092 1093;; Integral of the Bessel K function wrt z 1094;; http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/21/01/01/ 1095(defun bessel-k-integral-2 (n z) 1096 (cond 1097 ((and ($integerp n) (<= 0 n)) 1098 (cond 1099 (($oddp n) 1100 ;; integrate(bessel_k(2*N+1,z),z) , N > 0 1101 ;; = -(-1)^((n-1)/2)*bessel_k(0,z) 1102 ;; + 2*sum((-1)^(k+(n-1)/2-1)*bessel_k(2*k,z),k,1,(n-1)/2) 1103 (let* ((k (gensym)) 1104 (answer `((mplus) 1105 ((mtimes) -1 ((%bessel_k) 0 ,z) 1106 ((mexpt) -1 1107 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))) 1108 ((mtimes) 2 1109 ((%sum) 1110 ((mtimes) ((%bessel_k) ((mtimes) 2 ,k) ,z) 1111 ((mexpt) -1 1112 ((mplus) -1 ,k 1113 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n))))) 1114 ,k 1 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n))))))) 1115 ;; Expand out the sum if n < 10. Otherwise fix up the indices 1116 (if (< n 10) 1117 (meval `(($ev) ,answer $sum)) ; Is there a better way? 1118 (simplify ($niceindices answer))))) 1119 (($evenp n) 1120 ;; integrate(bessel_k(2*N,z),z) , N > 0 1121 ;; = (1/2)*(-1)^(n/2)*%pi*z*(bessel_k(0,z)*struve_l(-1,z) 1122 ;; +bessel_k(1,z)*struve_l(0,z)) 1123 ;; + 2 * sum((-1)^(k+n/2)*bessel_k(2*k+1,z),k,0,n/2-1) 1124 (let* 1125 ((k (gensym)) 1126 (answer `((mplus) 1127 ((mtimes) 2 1128 ((%sum) 1129 ((mtimes) 1130 ((%bessel_k) ((mplus) 1 ((mtimes) 2 ,k)) ,z) 1131 ((mexpt) -1 1132 ((mplus) ,k ((mtimes) ((rat) 1 2) ,n)))) 1133 ,k 0 ((mplus) -1 ((mtimes) ((rat) 1 2) ,n)))) 1134 ((mtimes) 1135 ((rat) 1 2) 1136 $%pi 1137 ((mexpt) -1 ((mtimes) ((rat) 1 2) ,n)) 1138 z 1139 ((mplus) 1140 ((mtimes) 1141 ((%bessel_k) 0 ,z) 1142 ((%struve_l) -1 ,z)) 1143 ((mtimes) 1144 ((%bessel_k) 1 ,z) 1145 ((%struve_l) 0 ,z))))))) 1146 ;; expand out the sum if n < 10. Otherwise fix up the indices 1147 (if (< n 10) 1148 (meval `(($ev) ,answer $sum)) ; Is there a better way? 1149 (simplify ($niceindices answer))))))) 1150 (t nil))) 1151 1152(putprop '%bessel_k `((n z) nil ,#'bessel-k-integral-2) 'integral) 1153 1154;; Support a simplim%bessel_k function to handle specific limits 1155 1156(defprop %bessel_k simplim%bessel_k simplim%function) 1157 1158(defun simplim%bessel_k (expr var val) 1159 ;; Look for the limit of the arguments. 1160 (let ((v (limit (cadr expr) var val 'think)) 1161 (z (limit (caddr expr) var val 'think))) 1162 (cond 1163 ;; Handle an argument 0 at this place. 1164 ((or (zerop1 z) 1165 (eq z '$zeroa) 1166 (eq z '$zerob)) 1167 (cond ((zerop1 v) 1168 ;; bessel_k(0,0) 1169 '$inf) 1170 ((integerp v) 1171 ;; bessel_k(n,0), n an integer 1172 (cond ((evenp v) '$inf) 1173 (t (cond ((eq z '$zeroa) '$inf) 1174 ((eq z '$zerob) '$minf) 1175 (t '$infinity))))) 1176 ((not (zerop1 ($realpart v))) 1177 ;; bessel_k(v,0), Re(v)#0 1178 '$infinity) 1179 ((and (zerop1 ($realpart v)) 1180 (not (zerop1 v))) 1181 ;; bessel_k(v,0), Re(v)=0 and v#0 1182 '$und) 1183 ;; Call the simplifier of the function. 1184 (t (simplify (list '(%bessel_k) v z))))) 1185 ((eq z '$inf) 1186 ;; bessel_k(v,inf) 1187 0) 1188 ((eq z '$minf) 1189 ;; bessel_k(v,minf) 1190 '$infinity) 1191 (t 1192 ;; All other cases are handled by the simplifier of the function. 1193 (simplify (list '(%bessel_k) v z)))))) 1194 1195(defun simp-bessel-k (expr ignored z) 1196 (declare (ignore ignored)) 1197 (twoargcheck expr) 1198 (let ((order (simpcheck (cadr expr) z)) 1199 (arg (simpcheck (caddr expr) z)) 1200 (rat-order nil)) 1201 (cond 1202 ((zerop1 arg) 1203 ;; Domain error for a zero argument. 1204 (simp-domain-error 1205 (intl:gettext "bessel_k: bessel_k(~:M,~:M) is undefined.") order arg)) 1206 1207 ((complex-float-numerical-eval-p order arg) 1208 (cond ((= 0 ($imagpart order)) 1209 ;; order is real, arg is real or complex 1210 (complexify (bessel-k ($float order) 1211 (complex ($float ($realpart arg)) 1212 ($float ($imagpart arg)))))) 1213 (t 1214 ;; order is complex, arg is real or complex 1215 (let (($numer t) 1216 ($float t) 1217 (order ($float order)) 1218 (arg ($float arg))) 1219 ($float 1220 ($rectform 1221 (bessel-k-hypergeometric order arg))))))) 1222 1223 ((mminusp order) 1224 ;; A&S 9.6.6: K[-v](x) = K[v](x) 1225 (take '(%bessel_k) (mul -1 order) arg)) 1226 1227 ((and $besselexpand 1228 (setq rat-order (max-numeric-ratio-p order 2))) 1229 ;; When order is a fraction with a denominator of 2, we 1230 ;; can express the result in terms of elementary 1231 ;; functions. From A&S 10.2.16 and 10.2.17: 1232 ;; 1233 ;; K[1/2](z) = sqrt(%pi/2/z)*exp(-z) 1234 ;; = K[-1/2](z) 1235 (bessel-k-half-order rat-order arg)) 1236 1237 ((and $bessel_reduce 1238 (and (integerp order) 1239 (plusp order) 1240 (> order 1))) 1241 ;; Reduce a bessel function of order > 2 to order 1 and 0. 1242 ;; A&S 9.6.26: bessel_k(v,z) = 2*(v-1)/z*bessel_k(v-1,z)+bessel_k(v-2,z) 1243 (add (mul 2 1244 (- order 1) 1245 (inv arg) 1246 (take '(%bessel_k) (- order 1) arg)) 1247 (take '(%bessel_k) (- order 2) arg))) 1248 1249 ($hypergeometric_representation 1250 (bessel-k-hypergeometric order arg)) 1251 1252 (t 1253 (eqtest (list '(%bessel_k) order arg) expr))))) 1254 1255;; Returns the hypergeometric representation of bessel_k 1256(defun bessel-k-hypergeometric (order arg) 1257 ;; http://functions.wolfram.com/03.04.26.0002.01 1258 ;; hypergeometric([],[1-v],z^2/4)*2^(v-1)*gamma(v)/z^v 1259 ;; + hypergeometric([],[v+1],z^2/4)*2^(-v-1)*gamma(-v)*z^v 1260 (add (mul (power 2 (sub order 1)) 1261 (take '(%gamma) order) 1262 (inv (power arg order)) 1263 (take '($hypergeometric) 1264 (list '(mlist)) 1265 (list '(mlist) (sub 1 order)) 1266 (div (mul arg arg) 4))) 1267 (mul (inv (power 2 (add order 1))) 1268 (power arg order) 1269 (take '(%gamma) (neg order)) 1270 (take '($hypergeometric) 1271 (list '(mlist)) 1272 (list '(mlist) (add order 1)) 1273 (div (mul arg arg) 4))))) 1274 1275;; Compute value of Modified Bessel function of the second kind of order ORDER 1276(defun bessel-k (order arg) 1277 (cond 1278 ((zerop (imagpart arg)) 1279 ;; We have numeric args and the first arg is purely real. Call the 1280 ;; real-valued Bessel function. Handle orders 0 and 1 specially, when 1281 ;; possible. 1282 (let ((arg (realpart arg))) 1283 (cond 1284 ((< arg 0) 1285 ;; This is the extension for negative arg. 1286 ;; We use A&S 9.6.6 and 9.6.31 for evaluation: 1287 ;; K[v](-z) = K[|v|](-z) 1288 ;; = exp(-%i*%pi*|v|)*K[|v|](z) - %i*%pi*I[|v|](z) 1289 (let* ((dpi (coerce pi 'flonum)) 1290 (s1 (cis (* dpi (- (abs order))))) 1291 (s2 (* (complex 0 -1) dpi)) 1292 (result (+ (* s1 (bessel-k (abs order) (- arg))) 1293 (* s2 (bessel-i (abs order) (- arg))))) 1294 (rem (nth-value 1 (floor order)))) 1295 (cond ((zerop rem) 1296 ;; order is an integer or a float representation of an 1297 ;; integer, the result is a general complex 1298 result) 1299 ((= rem 1/2) 1300 ;; order is half-integral-value or an float representation 1301 ;; and arg < 0, the result is pure imaginary 1302 (complex 0 (imagpart result))) 1303 ;; in all other cases general complex result 1304 (t result)))) 1305 ((= order 0) 1306 (slatec:dbesk0 (float arg))) 1307 ((= order 1) 1308 (slatec:dbesk1 (float arg))) 1309 (t 1310 ;; From A&S 9.6.6, K[-v](z) = K[v](z), so take the 1311 ;; absolute value of the order. 1312 (multiple-value-bind (n alpha) (floor (abs (float order))) 1313 (let ((jvals (make-array (1+ n) :element-type 'flonum))) 1314 (slatec:dbesk (float arg) alpha 1 (1+ n) jvals 0) 1315 (aref jvals n))))))) 1316 (t 1317 ;; The arg is complex. Use the complex-valued Bessel function. From 1318 ;; A&S 9.6.6, K[-v](z) = K[v](z), so take the absolute value of the order. 1319 (multiple-value-bind (n alpha) (floor (abs (float order))) 1320 (let ((cyr (make-array (1+ n) :element-type 'flonum)) 1321 (cyi (make-array (1+ n) :element-type 'flonum))) 1322 (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n 1323 v-cyr v-cyi v-nz v-ierr) 1324 (slatec::zbesk (float (realpart arg)) 1325 (float (imagpart arg)) 1326 alpha 1 (1+ n) cyr cyi 0 0) 1327 (declare (ignore v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz)) 1328 1329 ;; We should check for errors here based on the value of v-ierr. 1330 (when (plusp v-ierr) 1331 (format t "zbesk ierr = ~A~%" v-ierr)) 1332 (complex (aref cyr n) (aref cyi n)))))))) 1333 1334;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1335;;; 1336;;; Expansion of Bessel J and Y functions of half-integral order. 1337;;; 1338;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1339;; 1340;; From A&S 10.1.1, we have 1341;; 1342;; J[n+1/2](z) = sqrt(2*z/%pi)*j[n](z) 1343;; Y[n+1/2](z) = sqrt(2*z/%pi)*y[n](z) 1344;; 1345;; where j[n](z) is the spherical bessel function of the first kind 1346;; and y[n](z) is the spherical bessel function of the second kind. 1347;; 1348;; A&S 10.1.8 and 10.1.9 give 1349;; 1350;; 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)] 1351;; 1352;; 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)] 1353;; 1354 1355;; A&S 10.1.10 1356;; 1357;; j[n](z) = f[n](z)*sin(z) + (-1)^(n+1)*f[-n-1](z)*cos(z) 1358;; 1359;; f[0](z) = 1/z, f[1](z) = 1/z^2 1360;; 1361;; f[n-1](z) + f[n+1](z) = (2*n+1)/z*f[n](z) 1362;; 1363(defun f-fun (n z) 1364 (cond ((= n 0) 1365 (div 1 z)) 1366 ((= n 1) 1367 (div 1 (mul z z))) 1368 ((= n -1) 1369 0) 1370 ((>= n 2) 1371 ;; f[n+1](z) = (2*n+1)/z*f[n](z) - f[n-1](z) or 1372 ;; f[n](z) = (2*n-1)/z*f[n-1](z) - f[n-2](z) 1373 (sub (mul (div (+ n n -1) z) 1374 (f-fun (1- n) z)) 1375 (f-fun (- n 2) z))) 1376 (t 1377 ;; Negative n 1378 ;; 1379 ;; f[n-1](z) = (2*n+1)/z*f[n](z) - f[n+1](z) or 1380 ;; f[n](z) = (2*n+3)/z*f[n+1](z) - f[n+2](z) 1381 (sub (mul (div (+ n n 3) z) 1382 (f-fun (1+ n) z)) 1383 (f-fun (+ n 2) z))))) 1384 1385;; Compute sqrt(2*z/%pi) 1386(defun root-2z/pi (z) 1387 (power (mul 2 z (inv '$%pi)) '((rat simp) 1 2))) 1388 1389(defun bessel-j-half-order (order arg) 1390 "Compute J[n+1/2](z)" 1391 (let* ((n (floor order)) 1392 (sign (if (oddp n) -1 1)) 1393 (jn (sub (mul ($expand (f-fun n arg)) 1394 (take '(%sin) arg)) 1395 (mul sign 1396 ($expand (f-fun (- (- n) 1) arg)) 1397 (take '(%cos) arg))))) 1398 (mul (root-2z/pi arg) 1399 jn))) 1400 1401(defun bessel-y-half-order (order arg) 1402 "Compute Y[n+1/2](z)" 1403 ;; A&S 10.1.1: 1404 ;; Y[n+1/2](z) = sqrt(2*z/%pi)*y[n](z) 1405 ;; 1406 ;; A&S 10.1.15: 1407 ;; y[n](z) = (-1)^(n+1)*j[-n-1](z) 1408 ;; 1409 ;; So 1410 ;; Y[n+1/2](z) = sqrt(2*z/%pi)*(-1)^(n+1)*j[-n-1](z) 1411 ;; = (-1)^(n+1)*sqrt(2*z/%pi)*j[-n-1](z) 1412 ;; = (-1)^(n+1)*J[-n-1/2](z) 1413 (let* ((n (floor order)) 1414 (jn (bessel-j-half-order (- (- order) 1/2) arg))) 1415 (if (evenp n) 1416 (mul -1 jn) 1417 jn))) 1418 1419;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1420;;; 1421;;; Expansion of Bessel I and K functions of half-integral order. 1422;;; 1423;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1424 1425;; See A&S 10.2.12 1426;; 1427;; I[n+1/2](z) = sqrt(2*z/%pi)*[g[n](z)*sinh(z) + g[-n-1](z)*cosh(z)] 1428;; 1429;; g[0](z) = 1/z, g[1](z) = -1/z^2 1430;; 1431;; g[n-1](z) - g[n+1](z) = (2*n+1)/z*g[n](z) 1432;; 1433;; 1434(defun g-fun (n z) 1435 (declare (type integer n)) 1436 (cond ((= n 0) 1437 (div 1 z)) 1438 ((= n 1) 1439 (div -1 (mul z z))) 1440 ((>= n 2) 1441 ;; g[n](z) = g[n-2](z) - (2*n-1)/z*g[n-1](z) 1442 (sub (g-fun (- n 2) z) 1443 (mul (div (+ n n -1) z) 1444 (g-fun (- n 1) z)))) 1445 (t 1446 ;; n is negative 1447 ;; 1448 ;; g[n](z) = (2*n+3)/z*g[n+1](z) + g[n+2](z) 1449 (add (mul (div (+ n n 3) z) 1450 (g-fun (+ n 1) z)) 1451 (g-fun (+ n 2) z))))) 1452 1453;; See A&S 10.2.12 1454;; 1455;; I[n+1/2](z) = sqrt(2*z/%pi)*[g[n](z)*sinh(z) + g[-n-1](z)*cosh(z)] 1456 1457(defun bessel-i-half-order (order arg) 1458 (let ((order (floor order))) 1459 (mul (root-2z/pi arg) 1460 (add (mul ($expand (g-fun order arg)) 1461 (take '(%sinh) arg)) 1462 (mul ($expand (g-fun (- (+ order 1)) arg)) 1463 (take '(%cosh) arg)))))) 1464 1465;; See A&S 10.2.15 1466;; 1467;; 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) 1468;; 1469;; or 1470;; 1471;; K[n+1/2](z) = sqrt(%pi/2)/sqrt(z)*exp(-z)*sum((n+1/2,k)/(2*z)^k,k,0,n) 1472;; 1473;; where (A&S 10.1.9) 1474;; 1475;; (n+1/2,k) = (n+k)!/k!/(n-k)! 1476 1477(defun k-fun (n z) 1478 (declare (type unsigned-byte n)) 1479 ;; Computes the sum above 1480 (let ((sum 1) 1481 (term 1)) 1482 (loop for k from 0 upto n do 1483 (setf term (mul term 1484 (div (div (* (- n k) (+ n k 1)) 1485 (+ k 1)) 1486 (mul 2 z)))) 1487 (setf sum (add sum term))) 1488 sum)) 1489 1490(defun bessel-k-half-order (order arg) 1491 (let ((order (truncate (abs order)))) 1492 (mul (mul (power (div '$%pi 2) '((rat simp) 1 2)) 1493 (power arg '((rat simp) -1 2)) 1494 (power '$%e (neg arg))) 1495 (k-fun (abs order) arg)))) 1496 1497;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1498;;; 1499;;; Realpart und imagpart of Bessel functions 1500;;; 1501;;; We handle the special cases when we know that the Bessel function 1502;;; is pure real or pure imaginary. In all other cases Maxima generates 1503;;; a general noun form as result. 1504;;; 1505;;; To get the complex sign of the order an argument of the Bessel function 1506;;; the function $csign is used which calls $sign in a complex mode. 1507;;; 1508;;; This is an extension of of the algorithm of the function risplit in 1509;;; the file rpart.lisp. risplit looks for a risplit-function on the 1510;;; property list and call it if available. 1511;;; 1512;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1513 1514(defvar *debug-bessel* nil) 1515 1516;;; Put the risplit-function for Bessel J and Bessel I on the property list 1517 1518(defprop %bessel_j risplit-bessel-j-or-i risplit-function) 1519(defprop %bessel_i risplit-bessel-j-or-i risplit-function) 1520 1521;;; realpart und imagpart for Bessel J and Bessel I function 1522 1523(defun risplit-bessel-j-or-i (expr) 1524 (when *debug-bessel* 1525 (format t "~&RISPLIT-BESSEL-J with ~A~%" expr)) 1526 (let ((order (cadr expr)) 1527 (arg (caddr expr)) 1528 (sign-order ($csign (cadr expr))) 1529 (sign-arg ($csign (caddr expr)))) 1530 1531 (when *debug-bessel* 1532 (format t "~& : order = ~A~%" order) 1533 (format t "~& : arg = ~A~%" arg) 1534 (format t "~& : sign-order = ~A~%" sign-order) 1535 (format t "~& : sign-arg = ~A~%" sign-arg)) 1536 1537 (cond 1538 ((or (member sign-order '($complex $imaginary)) 1539 (eq sign-arg '$complex)) 1540 ;; order or arg are complex, return general noun-form 1541 (risplit-noun expr)) 1542 ((eq sign-arg '$imaginary) 1543 ;; arg is pure imaginary 1544 (cond 1545 ((or ($oddp order) 1546 ($featurep order '$odd)) 1547 ;; order is an odd integer, pure imaginary noun-form 1548 (cons 0 (mul -1 '$%i expr))) 1549 ((or ($evenp order) 1550 ($featurep order '$even)) 1551 ;; order is an even integer, real noun-form 1552 (cons expr 0)) 1553 (t 1554 ;; order is not an odd or even integer, or Maxima can not 1555 ;; determine it, return general noun-form 1556 (risplit-noun expr)))) 1557 ;; At this point order and arg are real. 1558 ;; We have to look for some special cases involing a negative arg 1559 ((or (maxima-integerp order) 1560 (member sign-arg '($pos $pz))) 1561 ;; arg is positive or order an integer, real noun-form 1562 (cons expr 0)) 1563 ;; At this point we know that arg is negative or the sign is not known 1564 ((zerop1 (sub ($truncate ($multthru 2 order)) ($multthru 2 order))) 1565 ;; order is half integral 1566 (cond 1567 ((eq sign-arg '$neg) 1568 ;; order is half integral and arg negative, imaginary noun-form 1569 (cons 0 (mul -1 '$%i expr))) 1570 (t 1571 ;; the sign of arg or order is not known 1572 (risplit-noun expr)))) 1573 (t 1574 ;; the sign of arg or order is not known 1575 (risplit-noun expr))))) 1576 1577;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1578 1579;;; Put the risplit-function for Bessel K and Bessel Y on the property list 1580 1581(defprop %bessel_k risplit-bessel-k-or-y risplit-function) 1582(defprop %bessel_y risplit-bessel-k-or-y risplit-function) 1583 1584;;; realpart und imagpart for Bessel K and Bessel Y function 1585 1586(defun risplit-bessel-k-or-y (expr) 1587 (when *debug-bessel* 1588 (format t "~&RISPLIT-BESSEL-K with ~A~%" expr)) 1589 (let ((order (cadr expr)) 1590 (arg (caddr expr)) 1591 (sign-order ($csign (cadr expr))) 1592 (sign-arg ($csign (caddr expr)))) 1593 1594 (when *debug-bessel* 1595 (format t "~& : order = ~A~%" order) 1596 (format t "~& : arg = ~A~%" arg) 1597 (format t "~& : sign-order = ~A~%" sign-order) 1598 (format t "~& : sign-arg = ~A~%" sign-arg)) 1599 1600 (cond 1601 ((or (member sign-order '($complex $imaginary)) 1602 (member sign-arg '($complex '$imaginary))) 1603 (risplit-noun expr)) 1604 ;; At this point order and arg are real valued. 1605 ;; We have to look for some special cases involing a negative arg 1606 ((member sign-arg '($pos $pz)) 1607 ;; arg is positive 1608 (cons expr 0)) 1609 ;; At this point we know that arg is negative or the sign is not known 1610 ((and (not (maxima-integerp order)) 1611 (zerop1 (sub ($truncate ($multthru 2 order)) ($multthru 2 order)))) 1612 ;; order is half integral 1613 (cond 1614 ((eq sign-arg '$neg) 1615 (cons 0 (mul -1 '$%i expr))) 1616 (t 1617 ;; the sign of arg is not known 1618 (risplit-noun expr)))) 1619 (t 1620 ;; the sign of arg is not known 1621 (risplit-noun expr))))) 1622 1623;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1624;;; 1625;;; Implementation of scaled Bessel I functions 1626;;; 1627;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1628 1629;; FIXME: The following scaled functions need work. They should be 1630;; extended to match bessel_i, but still carefully compute the scaled 1631;; value. 1632 1633;; I think g0(x) = exp(-x)*I[0](x), g1(x) = exp(-x)*I[1](x), and 1634;; gn(x,n) = exp(-x)*I[n](x), based on some simple numerical 1635;; evaluations. 1636 1637(defmfun $scaled_bessel_i0 ($x) 1638 (cond ((mnump $x) 1639 ;; XXX Should we return noun forms if $x is rational? 1640 (slatec:dbsi0e ($float $x))) 1641 (t 1642 (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil))) 1643 `((%bessel_i) 0 ,$x))))) 1644 1645(defmfun $scaled_bessel_i1 ($x) 1646 (cond ((mnump $x) 1647 ;; XXX Should we return noun forms if $x is rational? 1648 (slatec:dbsi1e ($float $x))) 1649 (t 1650 (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil))) 1651 `((%bessel_i) 1 ,$x))))) 1652 1653(defmfun $scaled_bessel_i ($n $x) 1654 (cond ((and (mnump $x) (mnump $n)) 1655 ;; XXX Should we return noun forms if $n and $x are rational? 1656 (multiple-value-bind (n alpha) (floor ($float $n)) 1657 (let ((iarray (make-array (1+ n) :element-type 'flonum))) 1658 (slatec:dbesi ($float $x) alpha 2 (1+ n) iarray 0) 1659 (aref iarray n)))) 1660 (t 1661 (mul (power '$%e (neg (simplifya `((mabs) ,$x) nil))) 1662 ($bessel_i $n $x))))) 1663 1664;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1665;;; 1666;;; Implementation of the Hankel 1 function 1667;;; 1668;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1669 1670(defmfun $hankel_1 (v z) 1671 (simplify (list '(%hankel_1) v z))) 1672 1673(defprop $hankel_1 %hankel_1 alias) 1674(defprop $hankel_1 %hankel_1 verb) 1675(defprop %hankel_1 $hankel_1 reversealias) 1676(defprop %hankel_1 $hankel_1 noun) 1677 1678;; hankel_1 distributes over lists, matrices, and equations 1679 1680(defprop %hankel_1 (mlist $matrix mequal) distribute_over) 1681 1682(defprop %hankel_1 simp-hankel-1 operators) 1683 1684; Derivatives of the Hankel 1 function 1685(defprop %hankel_1 1686 ((n x) 1687 nil 1688 1689 ;; Derivative wrt arg x. A&S 9.1.27. 1690 ((mtimes) 1691 ((mplus) 1692 ((%hankel_1) ((mplus) -1 n) x) 1693 ((mtimes) -1 ((%hankel_1) ((mplus) 1 n) x))) 1694 ((rat) 1 2))) 1695 grad) 1696 1697(defun simp-hankel-1 (expr ignored z) 1698 (declare (ignore ignored)) 1699 (twoargcheck expr) 1700 (let ((order (simpcheck (cadr expr) z)) 1701 (arg (simpcheck (caddr expr) z)) 1702 rat-order) 1703 (cond 1704 ((zerop1 arg) 1705 (simp-domain-error 1706 (intl:gettext "hankel_1: hankel_1(~:M,~:M) is undefined.") 1707 order arg)) 1708 ((complex-float-numerical-eval-p order arg) 1709 (cond ((= 0 ($imagpart order)) 1710 ;; order is real, arg is real or complex 1711 (complexify (hankel-1 ($float order) 1712 (complex ($float ($realpart arg)) 1713 ($float ($imagpart arg)))))) 1714 (t 1715 ;; The order is complex. Use 1716 ;; hankel_1(v,z) = bessel_j(v,z) + %i*bessel_y(v,z) 1717 ;; and evaluate using the hypergeometric function 1718 (let (($numer t) 1719 ($float t) 1720 (order ($float order)) 1721 (arg ($float arg))) 1722 ($float 1723 ($rectform 1724 (add (bessel-j-hypergeometric order arg) 1725 (mul '$%i 1726 (bessel-y-hypergeometric order arg))))))))) 1727 ((and $besselexpand 1728 (setq rat-order (max-numeric-ratio-p order 2))) 1729 ;; When order is a fraction with a denominator of 2, we can express 1730 ;; the result in terms of elementary functions. 1731 ;; Use the definition hankel_1(v,z) = bessel_j(v,z)+%i*bessel_y(v,z) 1732 (sratsimp 1733 (add (bessel-j-half-order rat-order arg) 1734 (mul '$%i 1735 (bessel-y-half-order rat-order arg))))) 1736 (t (eqtest (list '(%hankel_1) order arg) expr))))) 1737 1738;; Numerically compute H1[v](z). 1739;; 1740;; A&S 9.1.3 says H1[v](z) = J[v](z) + %i * Y[v](z) 1741;; 1742(defun hankel-1 (v z) 1743 (let ((v (float v)) 1744 (z (coerce z '(complex flonum)))) 1745 (cond ((minusp v) 1746 ;; A&S 9.1.6: 1747 ;; 1748 ;; H1[-v](z) = exp(v*%pi*%i)*H1[v](z) 1749 ;; 1750 ;; or 1751 ;; 1752 ;; H1[v](z) = exp(-v*%pi*%i)*H1[-v](z) 1753 1754 (* (cis (* (float pi) (- v))) (hankel-1 (- v) z))) 1755 (t 1756 (multiple-value-bind (n fnu) (floor v) 1757 (let ((zr (realpart z)) 1758 (zi (imagpart z)) 1759 (cyr (make-array (1+ n) :element-type 'flonum)) 1760 (cyi (make-array (1+ n) :element-type 'flonum))) 1761 (multiple-value-bind (dzr dzi df dk dm dn dcyr dcyi nz ierr) 1762 (slatec::zbesh zr zi fnu 1 1 (1+ n) cyr cyi 0 0) 1763 (declare (ignore dzr dzi df dk dm dn dcyr dcyi nz ierr)) 1764 (complex (aref cyr n) (aref cyi n))))))))) 1765 1766;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1767;;; 1768;;; Implementation of the Hankel 2 function 1769;;; 1770;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1771 1772(defmfun $hankel_2 (v z) 1773 (simplify (list '(%hankel_2) v z))) 1774 1775(defprop $hankel_2 %hankel_2 alias) 1776(defprop $hankel_2 %hankel_2 verb) 1777(defprop %hankel_2 $hankel_2 reversealias) 1778(defprop %hankel_2 $hankel_2 noun) 1779 1780;; hankel_2 distributes over lists, matrices, and equations 1781 1782(defprop %hankel_2 (mlist $matrix mequal) distribute_over) 1783 1784(defprop %hankel_2 simp-hankel-2 operators) 1785 1786; Derivatives of the Hankel 2 function 1787(defprop %hankel_2 1788 ((n x) 1789 nil 1790 1791 ;; Derivative wrt arg x. A&S 9.1.27. 1792 ((mtimes) 1793 ((mplus) 1794 ((%hankel_2) ((mplus) -1 n) x) 1795 ((mtimes) -1 ((%hankel_2) ((mplus) 1 n) x))) 1796 ((rat) 1 2))) 1797 grad) 1798 1799(defun simp-hankel-2 (expr ignored z) 1800 (declare (ignore ignored)) 1801 (twoargcheck expr) 1802 (let ((order (simpcheck (cadr expr) z)) 1803 (arg (simpcheck (caddr expr) z)) 1804 rat-order) 1805 (cond 1806 ((zerop1 arg) 1807 (simp-domain-error 1808 (intl:gettext "hankel_2: hankel_2(~:M,~:M) is undefined.") 1809 order arg)) 1810 ((complex-float-numerical-eval-p order arg) 1811 (cond ((= 0 ($imagpart order)) 1812 ;; order is real, arg is real or complex 1813 (complexify (hankel-2 ($float order) 1814 (complex ($float ($realpart arg)) 1815 ($float ($imagpart arg)))))) 1816 (t 1817 ;; The order is complex. Use 1818 ;; hankel_2(v,z) = bessel_j(v,z) - %i*bessel_y(v,z) 1819 ;; and evaluate using the hypergeometric function 1820 (let (($numer t) 1821 ($float t) 1822 (order ($float order)) 1823 (arg ($float arg))) 1824 ($float 1825 ($rectform 1826 (sub (bessel-j-hypergeometric order arg) 1827 (mul '$%i 1828 (bessel-y-hypergeometric order arg))))))))) 1829 ((and $besselexpand 1830 (setq rat-order (max-numeric-ratio-p order 2))) 1831 ;; When order is a fraction with a denominator of 2, we can express 1832 ;; the result in terms of elementary functions. 1833 ;; Use the definition hankel_2(v,z) = bessel_j(v,z)-%i*bessel_y(v,z) 1834 (sratsimp 1835 (sub (bessel-j-half-order rat-order arg) 1836 (mul '$%i 1837 (bessel-y-half-order rat-order arg))))) 1838 (t (eqtest (list '(%hankel_2) order arg) expr))))) 1839 1840;; Numerically compute H2[v](z). 1841;; 1842;; A&S 9.1.4 says H2[v](z) = J[v](z) - %i * Y[v](z) 1843;; 1844(defun hankel-2 (v z) 1845 (let ((v (float v)) 1846 (z (coerce z '(complex flonum)))) 1847 (cond ((minusp v) 1848 ;; A&S 9.1.6: 1849 ;; 1850 ;; H2[-v](z) = exp(-v*%pi*%i)*H2[v](z) 1851 ;; 1852 ;; or 1853 ;; 1854 ;; H2[v](z) = exp(v*%pi*%i)*H2[-v](z) 1855 1856 (* (cis (* (float pi) v)) (hankel-2 (- v) z))) 1857 (t 1858 (multiple-value-bind (n fnu) (floor v) 1859 (let ((zr (realpart z)) 1860 (zi (imagpart z)) 1861 (cyr (make-array (1+ n) :element-type 'flonum)) 1862 (cyi (make-array (1+ n) :element-type 'flonum))) 1863 (multiple-value-bind (dzr dzi df dk dm dn dcyr dcyi nz ierr) 1864 (slatec::zbesh zr zi fnu 1 2 (1+ n) cyr cyi 0 0) 1865 (declare (ignore dzr dzi df dk dm dn dcyr dcyi nz ierr)) 1866 (complex (aref cyr n) (aref cyi n))))))))) 1867 1868;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1869;;; 1870;;; Implementation of Struve H function 1871;;; 1872;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1873 1874(defmfun $struve_h (v z) 1875 (simplify (list '(%struve_h) v z))) 1876 1877(defprop $struve_h %struve_h alias) 1878(defprop $struve_h %struve_h verb) 1879(defprop %struve_h $struve_h reversealias) 1880(defprop %struve_h $struve_h noun) 1881 1882;; struve_h distributes over lists, matrices, and equations 1883 1884(defprop %struve_h (mlist $matrix mequal) distribute_over) 1885 1886(defprop %struve_h simp-struve-h operators) 1887 1888; Derivatives of the Struve H function 1889(defprop %struve_h 1890 ((v z) 1891 nil 1892 1893 ;; Derivative wrt arg z. A&S 12.1.10. 1894 ((mtimes) 1895 ((rat) 1 2) 1896 ((mplus) 1897 ((%struve_h) ((mplus) -1 v) z) 1898 ((mtimes) -1 ((%struve_h) ((mplus) 1 v) z)) 1899 ((mtimes) 1900 ((mexpt) $%pi ((rat) -1 2)) 1901 ((mexpt) 2 ((mtimes) -1 v)) 1902 ((mexpt) ((%gamma) ((mplus) ((rat) 3 2) v)) -1) 1903 ((mexpt) z v))))) 1904 grad) 1905 1906;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1907 1908(defun simp-struve-h (expr ignored z) 1909 (declare (ignore ignored)) 1910 (twoargcheck expr) 1911 (let ((order (simpcheck (cadr expr) z)) 1912 (arg (simpcheck (caddr expr) z))) 1913 (cond 1914 1915 ;; Check for special values 1916 1917 ((zerop1 arg) 1918 (cond ((eq ($csign (add order 1)) '$zero) 1919 ; http://functions.wolfram.com/03.09.03.0018.01 1920 (cond ((or ($bfloatp order) 1921 ($bfloatp arg)) 1922 ($bfloat (div 2 '$%pi))) 1923 ((or (floatp order) 1924 (floatp arg)) 1925 ($float (div 2 '$%pi))) 1926 (t 1927 (div 2 '$%pi)))) 1928 ((eq ($sign (add ($realpart order) 1)) '$pos) 1929 ; http://functions.wolfram.com/03.09.03.0001.01 1930 arg) 1931 ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz)) 1932 (simp-domain-error 1933 (intl:gettext "struve_h: struve_h(~:M,~:M) is undefined.") 1934 order arg)) 1935 (t 1936 (eqtest (list '(%struve_h) order arg) expr)))) 1937 1938 ;; Check for numerical evaluation 1939 1940 ((complex-float-numerical-eval-p order arg) 1941 ; A&S 12.1.21 1942 (let (($numer t) ($float t)) 1943 ($rectform 1944 ($float 1945 (mul 1946 ($rectform (power arg (add order 1.0))) 1947 ($rectform (inv (power 2.0 order))) 1948 (inv (power ($float '$%pi) 0.5)) 1949 (inv (simplify (list '(%gamma) (add order 1.5)))) 1950 (simplify (list '($hypergeometric) 1951 (list '(mlist) 1) 1952 (list '(mlist) '((rat simp) 3 2) 1953 (add order '((rat simp) 3 2))) 1954 (div (mul arg arg) -4.0)))))))) 1955 1956 ((complex-bigfloat-numerical-eval-p order arg) 1957 ; A&S 12.1.21 1958 (let (($ratprint nil) 1959 (arg ($bfloat arg)) 1960 (order ($bfloat order))) 1961 ($rectform 1962 ($bfloat 1963 (mul 1964 ($rectform (power arg (add order 1))) 1965 ($rectform (inv (power 2 order))) 1966 (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2)))) 1967 (inv (simplify (list '(%gamma) 1968 (add order ($bfloat '((rat simp) 3 2)))))) 1969 (simplify (list '($hypergeometric) 1970 (list '(mlist) 1) 1971 (list '(mlist) '((rat simp) 3 2) 1972 (add order '((rat simp) 3 2))) 1973 (div (mul arg arg) ($bfloat -4))))))))) 1974 1975 ;; Transformations and argument simplifications 1976 1977 ((and $besselexpand 1978 (ratnump order) 1979 (integerp (mul 2 order))) 1980 (cond 1981 ((eq ($sign order) '$pos) 1982 ;; Expansion of Struve H for a positive half integral order. 1983 (sratsimp 1984 (add 1985 (mul 1986 (inv (simplify (list '(mfactorial) (sub order 1987 '((rat simp) 1 2))))) 1988 (inv (power '$%pi '((rat simp) 1 2 ))) 1989 (power (div arg 2) (add order -1)) 1990 (let ((index (gensumindex))) 1991 (dosum 1992 (mul 1993 (simplify (list '($pochhammer) '((rat simp) 1 2) index)) 1994 (simplify (list '($pochhammer) 1995 (sub '((rat simp) 1 2) order) 1996 index)) 1997 (power (mul -1 arg arg (inv 4)) (mul -1 index))) 1998 index 0 (sub order '((rat simp) 1 2)) t))) 1999 (mul 2000 (power (div 2 '$%pi) '((rat simp) 1 2)) 2001 (power -1 (add order '((rat simp) 1 2))) 2002 (inv (power arg '((rat simp) 1 2))) 2003 (add 2004 (mul 2005 (simplify 2006 (list '(%sin) 2007 (add (mul '((rat simp) 1 2) 2008 '$%pi 2009 (add order '((rat simp) 1 2))) 2010 arg))) 2011 (let ((index (gensumindex))) 2012 (dosum 2013 (mul 2014 (power -1 index) 2015 (simplify (list '(mfactorial) 2016 (add (mul 2 index) 2017 order 2018 '((rat simp) -1 2)))) 2019 (inv (simplify (list '(mfactorial) (mul 2 index)))) 2020 (inv (simplify (list '(mfactorial) 2021 (add (mul -2 index) 2022 order 2023 '((rat simp) -1 2))))) 2024 (inv (power (mul 2 arg) (mul 2 index)))) 2025 index 0 2026 (simplify (list '($floor) 2027 (div (sub (mul 2 order) 1) 4))) 2028 t))) 2029 (mul 2030 (simplify (list '(%cos) 2031 (add (mul '((rat simp) 1 2) 2032 '$%pi 2033 (add order '((rat simp) 1 2))) 2034 arg))) 2035 (let ((index (gensumindex))) 2036 (dosum 2037 (mul 2038 (power -1 index) 2039 (simplify (list '(mfactorial) 2040 (add (mul 2 index) 2041 order 2042 '((rat simp) 1 2)))) 2043 (power (mul 2 arg) (mul -1 (add (mul 2 index) 1))) 2044 (inv (simplify (list '(mfactorial) 2045 (add (mul 2 index) 1)))) 2046 (inv (simplify (list '(mfactorial) 2047 (add (mul -2 index) 2048 order 2049 '((rat simp) -3 2)))))) 2050 index 0 2051 (simplify (list '($floor) 2052 (div (sub (mul 2 order) 3) 4))) 2053 t)))))))) 2054 2055 ((eq ($sign order) '$neg) 2056 ;; Expansion of Struve H for a negative half integral order. 2057 (sratsimp 2058 (add 2059 (mul 2060 (power (div 2 '$%pi) '((rat simp) 1 2)) 2061 (power -1 (add order '((rat simp) 1 2))) 2062 (inv (power arg '((rat simp) 1 2))) 2063 (add 2064 (mul 2065 (simplify (list '(%sin) 2066 (add 2067 (mul 2068 '((rat simp) 1 2) 2069 '$%pi 2070 (add order '((rat simp) 1 2))) 2071 arg))) 2072 (let ((index (gensumindex))) 2073 (dosum 2074 (mul 2075 (power -1 index) 2076 (simplify (list '(mfactorial) 2077 (add (mul 2 index) 2078 (neg order) 2079 '((rat simp) -1 2)))) 2080 (inv (simplify (list '(mfactorial) (mul 2 index)))) 2081 (inv (simplify (list '(mfactorial) 2082 (add (mul -2 index) 2083 (neg order) 2084 '((rat simp) -1 2))))) 2085 (inv (power (mul 2 arg) (mul 2 index)))) 2086 index 0 2087 (simplify (list '($floor) 2088 (div (add (mul 2 order) 1) -4))) 2089 t))) 2090 (mul 2091 (simplify (list '(%cos) 2092 (add 2093 (mul 2094 '((rat simp) 1 2) 2095 '$%pi 2096 (add order '((rat simp) 1 2))) 2097 arg))) 2098 (let ((index (gensumindex))) 2099 (dosum 2100 (mul 2101 (power -1 index) 2102 (simplify (list '(mfactorial) 2103 (add (mul 2 index) 2104 (neg order) 2105 '((rat simp) 1 2)))) 2106 (power (mul 2 arg) (mul -1 (add (mul 2 index) 1))) 2107 (inv (simplify (list '(mfactorial) 2108 (add (mul 2 index) 1)))) 2109 (inv (simplify (list '(mfactorial) 2110 (add (mul -2 index) 2111 (neg order) 2112 '((rat simp) -3 2)))))) 2113 index 0 2114 (simplify (list '($floor) 2115 (div (add (mul 2 order) 3) -4))) 2116 t)))))))))) 2117 (t 2118 (eqtest (list '(%struve_h) order arg) expr))))) 2119 2120;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2121;;; 2122;;; Implementation of Struve L function 2123;;; 2124;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2125 2126(defmfun $struve_l (v z) 2127 (simplify (list '(%struve_l) v z))) 2128 2129(defprop $struve_l %struve_l alias) 2130(defprop $struve_l %struve_l verb) 2131(defprop %struve_l $struve_l reversealias) 2132(defprop %struve_l $struve_l noun) 2133 2134(defprop %struve_l simp-struve-l operators) 2135 2136;; struve_l distributes over lists, matrices, and equations 2137 2138(defprop %struve_l (mlist $matrix mequal) distribute_over) 2139 2140; Derivatives of the Struve L function 2141(defprop %struve_l 2142 ((v z) 2143 nil 2144 2145 ;; Derivative wrt arg z. A&S 12.2.5. 2146 ((mtimes) 2147 ((rat) 1 2) 2148 ((mplus) 2149 ((%struve_l) ((mplus) -1 v) z) 2150 ((%struve_l) ((mplus) 1 v) z) 2151 ((mtimes) 2152 ((mexpt) $%pi ((rat) -1 2)) 2153 ((mexpt) 2 ((mtimes) -1 v)) 2154 ((mexpt) ((%gamma) ((mplus) ((rat) 3 2) v)) -1) 2155 ((mexpt) z v))))) 2156 grad) 2157 2158;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2159 2160(defun simp-struve-l (expr ignored z) 2161 (declare (ignore ignored)) 2162 (twoargcheck expr) 2163 (let ((order (simpcheck (cadr expr) z)) 2164 (arg (simpcheck (caddr expr) z))) 2165 (cond 2166 2167 ;; Check for special values 2168 2169 ((zerop1 arg) 2170 (cond ((eq ($csign (add order 1)) '$zero) 2171 ; http://functions.wolfram.com/03.10.03.0018.01 2172 (cond ((or ($bfloatp order) 2173 ($bfloatp arg)) 2174 ($bfloat (div 2 '$%pi))) 2175 ((or (floatp order) 2176 (floatp arg)) 2177 ($float (div 2 '$%pi))) 2178 (t 2179 (div 2 '$%pi)))) 2180 ((eq ($sign (add ($realpart order) 1)) '$pos) 2181 ; http://functions.wolfram.com/03.10.03.0001.01 2182 arg) 2183 ((member ($sign (add ($realpart order) 1)) '($zero $neg $nz)) 2184 (simp-domain-error 2185 (intl:gettext "struve_l: struve_l(~:M,~:M) is undefined.") 2186 order arg)) 2187 (t 2188 (eqtest (list '(%struve_l) order arg) expr)))) 2189 2190 ;; Check for numerical evaluation 2191 2192 ((complex-float-numerical-eval-p order arg) 2193 ; http://functions.wolfram.com/03.10.26.0002.01 2194 (let (($numer t) ($float t)) 2195 ($rectform 2196 ($float 2197 (mul 2198 ($rectform (power arg (add order 1.0))) 2199 ($rectform (inv (power 2.0 order))) 2200 (inv (power ($float '$%pi) 0.5)) 2201 (inv (simplify (list '(%gamma) (add order 1.5)))) 2202 (simplify (list '($hypergeometric) 2203 (list '(mlist) 1) 2204 (list '(mlist) '((rat simp) 3 2) 2205 (add order '((rat simp) 3 2))) 2206 (div (mul arg arg) 4.0)))))))) 2207 2208 ((complex-bigfloat-numerical-eval-p order arg) 2209 ; http://functions.wolfram.com/03.10.26.0002.01 2210 (let (($ratprint nil)) 2211 ($rectform 2212 ($bfloat 2213 (mul 2214 ($rectform (power arg (add order 1))) 2215 ($rectform (inv (power 2 order))) 2216 (inv (power ($bfloat '$%pi) ($bfloat '((rat simp) 1 2)))) 2217 (inv (simplify (list '(%gamma) (add order '((rat simp) 3 2))))) 2218 (simplify (list '($hypergeometric) 2219 (list '(mlist) 1) 2220 (list '(mlist) '((rat simp) 3 2) 2221 (add order '((rat simp) 3 2))) 2222 (div (mul arg arg) 4)))))))) 2223 2224 ;; Transformations and argument simplifications 2225 2226 ((and $besselexpand 2227 (ratnump order) 2228 (integerp (mul 2 order))) 2229 (cond 2230 ((eq ($sign order) '$pos) 2231 ;; Expansion of Struve L for a positive half integral order. 2232 (sratsimp 2233 (add 2234 (mul -1 2235 (power 2 (sub 1 order)) 2236 (power arg (sub order 1)) 2237 (inv (power '$%pi '((rat simp) 1 2))) 2238 (inv (simplify (list '(mfactorial) (sub order 2239 '((rat simp) 1 2))))) 2240 (let ((index (gensumindex))) 2241 (dosum 2242 (mul 2243 (simplify (list '($pochhammer) '((rat simp) 1 2) index)) 2244 (simplify (list '($pochhammer) 2245 (sub '((rat simp) 1 2) order) 2246 index)) 2247 (power (mul arg arg (inv 4)) (mul -1 index))) 2248 index 0 (sub order '((rat simp) 1 2)) t))) 2249 (mul -1 2250 (power (div 2 '$%pi) '((rat simp) 1 2)) 2251 (inv (power arg '((rat simp) 1 2))) 2252 ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2)) 2253 (add 2254 (mul 2255 (let (($trigexpand t)) 2256 (simplify (list '(%sinh) 2257 (sub (mul '((rat simp) 1 2) 2258 '$%pi 2259 '$%i 2260 (add order '((rat simp) 1 2))) 2261 arg)))) 2262 (let ((index (gensumindex))) 2263 (dosum 2264 (mul 2265 (simplify (list '(mfactorial) 2266 (add (mul 2 index) 2267 (simplify (list '(mabs) order)) 2268 '((rat simp) -1 2)))) 2269 (inv (simplify (list '(mfactorial) (mul 2 index)))) 2270 (inv (simplify (list '(mfactorial) 2271 (add (simplify (list '(mabs) 2272 order)) 2273 (mul -2 index) 2274 '((rat simp) -1 2))))) 2275 (inv (power (mul 2 arg) (mul 2 index)))) 2276 index 0 2277 (simplify (list '($floor) 2278 (div (sub (mul 2 2279 (simplify (list '(mabs) 2280 order))) 2281 1) 2282 4))) 2283 t))) 2284 (mul 2285 (let (($trigexpand t)) 2286 (simplify (list '(%cosh) 2287 (sub (mul '((rat simp) 1 2) 2288 '$%pi 2289 '$%i 2290 (add order '((rat simp) 1 2))) 2291 arg)))) 2292 (let ((index (gensumindex))) 2293 (dosum 2294 (mul 2295 (simplify (list '(mfactorial) 2296 (add (mul 2 index) 2297 (simplify (list '(mabs) order)) 2298 '((rat simp) 1 2)))) 2299 (power (mul 2 arg) (neg (add (mul 2 index) 1))) 2300 (inv (simplify (list '(mfactorial) 2301 (add (mul 2 index) 1)))) 2302 (inv (simplify (list '(mfactorial) 2303 (add (simplify (list '(mabs) 2304 order)) 2305 (mul -2 index) 2306 '((rat simp) -3 2)))))) 2307 index 0 2308 (simplify (list '($floor) 2309 (div (sub (mul 2 2310 (simplify (list '(mabs) 2311 order))) 2312 3) 2313 4))) 2314 t)))))))) 2315 ((eq ($sign order) '$neg) 2316 ;; Expansion of Struve L for a negative half integral order. 2317 (sratsimp 2318 (add 2319 (mul -1 2320 (power (div 2 '$%pi) '((rat simp) 1 2)) 2321 (inv (power arg '((rat simp) 1 2))) 2322 ($exp (div (mul '$%pi '$%i (add order '((rat simp) 1 2))) 2)) 2323 (add 2324 (mul 2325 (let (($trigexpand t)) 2326 (simplify (list '(%sinh) 2327 (sub (mul '((rat simp) 1 2) 2328 '$%pi 2329 '$%i 2330 (add order '((rat simp) 1 2))) 2331 arg)))) 2332 (let ((index (gensumindex))) 2333 (dosum 2334 (mul 2335 (simplify (list '(mfactorial) 2336 (add (mul 2 index) 2337 (neg order) 2338 '((rat simp) -1 2)))) 2339 (inv (simplify (list '(mfactorial) (mul 2 index)))) 2340 (inv (simplify (list '(mfactorial) 2341 (add (neg order) 2342 (mul -2 index) 2343 '((rat simp) -1 2))))) 2344 (inv (power (mul 2 arg) (mul 2 index)))) 2345 index 0 2346 (simplify (list '($floor) 2347 (div (add (mul 2 order) 1) -4))) 2348 t))) 2349 (mul 2350 (let (($trigexpand t)) 2351 (simplify (list '(%cosh) 2352 (sub (mul '((rat simp) 1 2) 2353 '$%pi 2354 '$%i 2355 (add order '((rat simp) 1 2))) 2356 arg)))) 2357 (let ((index (gensumindex))) 2358 (dosum 2359 (mul 2360 (simplify (list '(mfactorial) 2361 (add (mul 2 index) 2362 (neg order) 2363 '((rat simp) 1 2)))) 2364 (power (mul 2 arg) (neg (add (mul 2 index) 1))) 2365 (inv (simplify (list '(mfactorial) 2366 (add (mul 2 index) 1)))) 2367 (inv (simplify (list '(mfactorial) 2368 (add (neg order) 2369 (mul -2 index) 2370 '((rat simp) -3 2)))))) 2371 index 0 2372 (simplify (list '($floor) 2373 (div (add (mul 2 order) 3) -4))) 2374 t)))))))))) 2375 (t 2376 (eqtest (list '(%struve_l) order arg) expr))))) 2377 2378;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2379 2380(defmspec $gauss (form) 2381 (format t 2382"NOTE: The gauss function is superseded by random_normal in the `distrib' package. 2383Perhaps you meant to enter `~a'.~%" 2384 (print-invert-case (implode (mstring `(($random_normal) ,@ (cdr form)))))) 2385 '$done) 2386