327 | | ) |

328 | | |

| 327 | (let ((big-n-1/2 (+ big-n 1/2)) |

| 328 | (eps (epsilon z))) |

| 329 | (do* ((n 0 (1+ n)) |

| 330 | (term (* (big-a 0 v) |

| 331 | (big-i 0 big-n-1/2 z v)) |

| 332 | (* (big-a n v) |

| 333 | (big-i n big-n-1/2 z v))) |

| 334 | (sum term (+ sum term))) |

| 335 | ((<= (abs term) (* eps (abs sum))) |

| 336 | sum) |

| 337 | #+nil |

| 338 | (progn |

| 339 | (format t "n = ~D~%" n) |

| 340 | (format t " term = ~S~%" term) |

| 341 | (format t " sum = ~S~%" sum))))) |

| 342 | |

| 343 | ;; Series for bessel J: |

| 344 | ;; |

| 345 | ;; (z/2)^v*sum((-1)^k/Gamma(k+v+1)/k!*(z^2//4)^k, k, 0, inf) |

| 346 | (defun s-bessel-j (v z) |

| 347 | (with-floating-point-contagion (v z) |

| 348 | (let ((z2/4 (* z z 1/4)) |

| 349 | (eps (epsilon z))) |

| 350 | (do* ((k 0 (+ 1 k)) |

| 351 | (f (gamma (+ v 1)) |

| 352 | (* f (* k (+ v k)))) |

| 353 | (term (/ f) |

| 354 | (/ (* (- term) z2/4) f)) |

| 355 | (sum term (+ sum term))) |

| 356 | ((<= (abs term) (* eps (abs sum))) |

| 357 | (* sum (expt (* z 1/2) v))) |

| 358 | (format t "k = ~D~%" k) |

| 359 | (format t " term = ~S~%" term) |

| 360 | (format t " sum = ~S~%" sum))))) |

| 361 | |