Commit 6ab5226ca3e32b443d87934ec138ff0efc8aaecc

Authored by Raymond Toy
1 parent d795ba71

First cut at Bessel Y. Not working yet.

Showing 1 changed file with 36 additions and 0 deletions
qd-bessel.lisp
... ... @@ -376,6 +376,7 @@
376 376 (format t " term = ~S~%" term)
377 377 (format t " sum = ~S~%" sum))))))
378 378  
  379 +;;
379 380 ;; TODO:
380 381 ;; o For |z| <= 1 use the series.
381 382 ;; o Currently accuracy is not good for large z and half-integer
... ... @@ -402,6 +403,41 @@
402 403 (+ (/ -1 z)
403 404 (sum-ab big-n v z)
404 405 (sum-big-ia big-n v z)))))))))
  406 +
  407 +;; Bessel Y
  408 +;;
  409 +;; bessel_y(v, z) = 1/(2*%pi*%i)*(exp(-%i*v*%pi/2)*I(%i*v,z) - exp(%i*v*%pi/2)*I(-%i*z, v))
  410 +;; + z/v/%pi*((1-cos(v*%pi)/z) + S(N,z,v)*cos(v*%pi)-S(N,z,-v))
  411 +;;
  412 +;; where
  413 +;;
  414 +;; S(N,z,v) = sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum(exp(-k*z)*a[n](k,v),k,1,N),n,0,inf)
  415 +;; + sum(A[n](v)*I[n](N+1/2,z,v),n,0,inf)
  416 +;;
  417 +(defun bessel-y (v z)
  418 + (flet ((ipart (v z)
  419 + (let* ((iz (* #c(0 1) z))
  420 + (c+ (exp (* v (float-pi z) 1/2)))
  421 + (c- (exp (* v (float-pi z) -1/2)))
  422 + (i+ (exp-arc-i-2 iz v))
  423 + (i- (exp-arc-i-2 (- iz) v)))
  424 + (/ (- (* c- i+) (* c+ i-))
  425 + (* #c(0 2) (float-pi z)))))
  426 + (s (big-n z v)
  427 + (+ (sum-ab big-n v z)
  428 + (sum-big-ia big-n v z))))
  429 + (let* ((big-n 100)
  430 + (vpi (* v (float-pi z)))
  431 + (c (cos vpi)))
  432 + (+ (ipart v z)
  433 + (* (/ z vpi)
  434 + (+ (/ (- 1 c)
  435 + z)
  436 + (* c
  437 + (s big-n z v))
  438 + (- (s big-n z (- v)))))))))
  439 +
  440 +
405 441  
406 442 (defun paris-series (v z n)
407 443 (labels ((pochhammer (a k)
... ...