Changeset 6ab5226ca3e32b443d87934ec138ff0efc8aaecc

Timestamp:

04/11/12 20:32:23 (14 months ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

9fd2ebcbeed3ecae899f732b15fd279f6fb0f14f

Parents:

d795ba718dc53f591c82994811f50250aceec1d7

git-committer:

Raymond Toy <toy.raymond@…> (04/11/12 20:32:23)

Message:

First cut at Bessel Y. Not working yet.

Files:

• qd-bessel.lisp (modified) (2 diffs)

qd-bessel.lisp

@@ -377 +377 @@
           (format t " sum  = ~S~%" sum))))))
 
+;; 
 ;; TODO:
 ;;  o For |z| <= 1 use the series.
@@ -403 +404 @@
                       (sum-ab big-n v z)
                       (sum-big-ia big-n v z)))))))))
+
+;; Bessel Y
+;;
+;; 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))
+;;                   + z/v/%pi*((1-cos(v*%pi)/z) + S(N,z,v)*cos(v*%pi)-S(N,z,-v))
+;;
+;; where
+;;
+;;   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)
+;;               + sum(A[n](v)*I[n](N+1/2,z,v),n,0,inf)
+;;
+(defun bessel-y (v z)
+  (flet ((ipart (v z)
+           (let* ((iz (* #c(0 1) z))
+                  (c+ (exp (* v (float-pi z) 1/2)))
+                  (c- (exp (* v (float-pi z) -1/2)))
+                  (i+ (exp-arc-i-2 iz v))
+                  (i- (exp-arc-i-2 (- iz) v)))
+             (/ (- (* c- i+) (* c+ i-))
+                (* #c(0 2) (float-pi z)))))
+         (s (big-n z v)
+           (+ (sum-ab big-n v z)
+              (sum-big-ia big-n v z))))
+    (let* ((big-n 100)
+           (vpi (* v (float-pi z)))
+           (c (cos vpi)))
+      (+ (ipart v z)
+         (* (/ z vpi)
+            (+ (/ (- 1 c)
+                  z)
+               (* c
+                  (s big-n z v))
+               (- (s big-n z (- v)))))))))
+           
+  
 
 (defun paris-series (v z n)