Commit 7c5a3186070096ee93e16f2ddf51b2c84e7c5895

Authored by Raymond Toy
1 parent bba9f894

Correct some comments, remove unused code.

Showing 1 changed file with 14 additions and 7 deletions
qd-bessel.lisp
... ... @@ -37,7 +37,7 @@
37 37  
38 38 ;; B[k](p) = 1/2^(k+3/2)*integrate(exp(-p*u)*u^(k-1/2),u,0,1)
39 39 ;; = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p)
40   -;; = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)
  40 +;; = 1/2^(k+3/2)/p^(k+1/2) * G(k+1/2, p)
41 41 ;;
42 42 ;; where G(a,z) is the lower incomplete gamma function.
43 43 ;;
... ... @@ -109,6 +109,12 @@
109 109 ;; = conj(B[k](%i*z).
110 110 ;;
111 111 ;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real.
  112 +;;
  113 +;; Also note that when v is an integer of the form (2*m+1)/2, then
  114 +;; r[2*k+1](-2*%i*v) = r[2*k+1](-%i*(2*m+1))
  115 +;; = -%i*(2*m+1)*product(-(2*m+1)^2+(2*j-1)^2, j, 1, k)
  116 +;; so the product is zero when k >= m and the series I(p, q) is
  117 +;; finite.
112 118 (defun exp-arc-i (p q)
113 119 (let* ((sqrt2 (sqrt (float 2 (realpart p))))
114 120 (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
... ... @@ -144,14 +150,9 @@
144 150 (format t " sum - ~S~%" sum)))))
145 151  
146 152 (defun exp-arc-i-2 (p q)
147   - (let* ((sqrt2 (sqrt (float 2 (realpart p))))
148   - (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
149   - (v (* #c(0 -2) q))
  153 + (let* ((v (* #c(0 -2) q))
150 154 (v2 (expt v 2))
151 155 (eps (epsilon (realpart p))))
152   - (when *debug-exparc*
153   - (format t "sqrt2 = ~S~%" sqrt2)
154   - (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))
155 156 (do* ((k 0 (1+ k))
156 157 (bk (bk 0 p)
157 158 (bk k p))
... ... @@ -162,6 +163,12 @@
162 163 (* ratio bk))
163 164 (sum term (+ sum term)))
164 165 ((< (abs term) (* (abs sum) eps))
  166 + (when *debug-exparc*
  167 + (format t "Final k= ~D~%" k)
  168 + (format t " bk = ~S~%" bk)
  169 + (format t " ratio = ~S~%" ratio)
  170 + (format t " term = ~S~%" term)
  171 + (format t " sum - ~S~%" sum))
165 172 (* sum #c(0 2) (/ (exp p) q)))
166 173 (when *debug-exparc*
167 174 (format t "k = ~D~%" k)
... ...