Commit 581a0a08 by Raymond Toy

### Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct

parents 76ea9b13 78801d63
 ... ... @@ -179,7 +179,13 @@ (format t " sum - ~S~%" sum))))) ;; Not really just for Bessel J for integer orders, but in that case, ;; this is all that's needed to compute Bessel J. For other values, ;; this is just part of the computation needed. ;; ;; Compute ;; ;; 1/(2*%pi) * (exp(-%i*v*%pi/2) * I(%i*z, v) + exp(%i*v*%pi/2) * I(-%i*z, v)) (defun integer-bessel-j-exp-arc (v z) (let* ((iz (* #c(0 1) z)) (i+ (exp-arc-i-2 iz v))) ... ... @@ -257,6 +263,7 @@ ;; ;; sum(exp(-k*z)*a[n](k,v), k, 1, N) ;; #+nil (defun sum-an (big-n n v z) (let ((sum 0)) (loop for k from 1 upto big-n ... ... @@ -265,6 +272,20 @@ (an n k v)))) sum)) ;; Like above, but we just stop when the terms no longer contribute to ;; the sum. (defun sum-an (big-n n v z) (let ((eps (epsilon (realpart z)))) (do* ((k 1 (+ 1 k)) (term (* (exp (- (* k z))) (an n k v)) (* (exp (- (* k z))) (an n k v))) (sum term (+ sum term))) ((or (<= (abs term) (* eps (abs sum))) (> k big-n)) sum)))) ;; SUM-AB computes the series ;; ;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf) ... ... @@ -399,7 +420,7 @@ (integer-bessel-j-exp-arc v z)) (t ;; Need to fine-tune the value of big-n. (let ((big-n 100) (let ((big-n 10) (vpi (* v (float-pi (realpart z))))) (+ (integer-bessel-j-exp-arc v z) (if (= vv v) ... ...
 ... ... @@ -1673,9 +1673,54 @@ nil) ;; Bessel J for complex args (rt:deftest bessel-j-complex.pos-order.d.1 (rt:deftest bessel-j-complex-arg.d.1 (let ((b (bessel-j 0d0 #c(1d0 1))) (true #c(0.9376084768060293d0 -0.4965299476091221d0))) (check-accuracy 50.73 b true)) nil) (rt:deftest bessel-j-complex-arg.d.2 (let ((b (bessel-j 1d0 #c(1d0 1))) (true #c(0.6141603349229036d0 0.3650280288270878d0))) (check-accuracy 52.51 b true)) nil) (rt:deftest bessel-j-complex-arg.d.3 (let ((b (bessel-j 2d0 #c(1d0 1))) (true #c(0.0415798869439621d0 0.2473976415133063d0))) (check-accuracy 50.41 b true)) nil) (rt:deftest bessel-j-complex-arg.d.4 (let ((b (bessel-j 2.3d0 #c(1d0 1))) (true #c(-0.0141615213034667d0 0.1677798241687935d0))) (check-accuracy 48.56 b true)) nil) (rt:deftest bessel-j-complex-arg.d.5 (let ((b (bessel-j -2.3d0 #c(1d0 1))) (true #c(0.1920598664138632d0 -0.5158676904105332d0))) (check-accuracy 50.97 b true)) nil) (rt:deftest bessel-j-1/2-complex.d.1 (loop for k from 0 below 10 for x = (complex (random (/ pi 2)) (random (/ pi 2))) for b = (bessel-j 0.5d0 x) for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))) for result = (check-accuracy 49.8 b true) when result append (list (list (list k x) result))) nil) (rt:deftest bessel-j-1/2-complex.q.1 (loop for k from 0 below 10 for x = (complex (random (/ (float-pi #q1) 2)) (random (/ (float-pi #q1) 2))) for b = (bessel-j #q0.5 x) for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) for result = (check-accuracy 212 b true) when result append (list (list (list k x) result))) nil)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!