Changeset 78801d


Ignore:
Timestamp:
04/15/12 17:50:16 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
8c5195, b1d9be
Parents:
5566bc
Message:
  • Add comments for integer-bessel-j-exp-arc.
  • Simplify sum-an so we stop the sum when the terms no longer contribute to the sum.
  • Change big-n. This still needs work.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • qd-bessel.lisp

    r235ac2 r78801d  
    180180
    181181
    182 ;;
     182;; Not really just for Bessel J for integer orders, but in that case,
     183;; this is all that's needed to compute Bessel J.  For other values,
     184;; this is just part of the computation needed.
     185;;
     186;; Compute
     187;;
     188;;  1/(2*%pi) * (exp(-%i*v*%pi/2) * I(%i*z, v) + exp(%i*v*%pi/2) * I(-%i*z, v))
    183189(defun integer-bessel-j-exp-arc (v z)
    184190  (let* ((iz (* #c(0 1) z))
     
    258264;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
    259265;;
     266#+nil
    260267(defun sum-an (big-n n v z)
    261268  (let ((sum 0))
     
    265272                          (an n k v))))
    266273    sum))
     274
     275;; Like above, but we just stop when the terms no longer contribute to
     276;; the sum.
     277(defun sum-an (big-n n v z)
     278  (let ((eps (epsilon (realpart z))))
     279    (do* ((k 1 (+ 1 k))
     280          (term (* (exp (- (* k z)))
     281                   (an n k v))
     282                (* (exp (- (* k z)))
     283                   (an n k v)))
     284          (sum term (+ sum term)))
     285         ((or (<= (abs term) (* eps (abs sum)))
     286              (> k big-n))
     287          sum))))
    267288
    268289;; SUM-AB computes the series
     
    400421          (t
    401422           ;; Need to fine-tune the value of big-n.
    402            (let ((big-n 100)
     423           (let ((big-n 10)
    403424                 (vpi (* v (float-pi (realpart z)))))
    404425             (+ (integer-bessel-j-exp-arc v z)
Note: See TracChangeset for help on using the changeset viewer.