Function: calcFunc-besJ

calcFunc-besJ is an autoloaded and byte-compiled function defined in calc-funcs.el.gz.

Signature

(calcFunc-besJ V X)

Source Code

;; Defined in /usr/src/emacs/lisp/calc/calc-funcs.el.gz
;;; Bessel functions.

;;; Should generalize this to handle arbitrary precision!

(defun calcFunc-besJ (v x)
  (or (math-numberp v) (math-reject-arg v 'numberp))
  (or (math-numberp x) (math-reject-arg x 'numberp))
  (let ((calc-internal-prec (min 8 calc-internal-prec)))
    (math-with-extra-prec 3
      (setq x (math-float (math-normalize x)))
      (setq v (math-float (math-normalize v)))
      (cond ((math-zerop x)
	     (if (math-zerop v)
		 '(float 1 0)
	       '(float 0 0)))
	    ((math-inexact-result))
	    ((not (math-num-integerp v))
	     (let ((start (math-div 1 (calcFunc-fact v))))
	       (math-mul (math-besJ-series start start
					   0
					   (math-mul '(float -25 -2)
						     (math-sqr x))
					   v)
			 (math-pow (math-div x 2) v))))
	    ((math-negp (setq v (math-trunc v)))
	     (if (math-oddp v)
		 (math-neg (calcFunc-besJ (math-neg v) x))
	       (calcFunc-besJ (math-neg v) x)))
	    ((eq v 0)
	     (math-besJ0 x))
	    ((eq v 1)
	     (math-besJ1 x))
	    ((Math-lessp v (math-abs-approx x))
	     (let ((j 0)
		   (bjm (math-besJ0 x))
		   (bj (math-besJ1 x))
		   (two-over-x (math-div 2 x))
		   bjp)
	       (while (< (setq j (1+ j)) v)
		 (setq bjp (math-sub (math-mul (math-mul j two-over-x) bj)
				     bjm)
		       bjm bj
		       bj bjp))
	       bj))
	    (t
	     (if (Math-lessp 100 v) (math-reject-arg v 'range))
	     (let* ((j (logior (+ v (cl-isqrt (* 40 v))) 1))
		    (two-over-x (math-div 2 x))
		    (jsum nil)
		    (bjp '(float 0 0))
		    (sum '(float 0 0))
		    (bj '(float 1 0))
		    bjm ans)
	       (while (> (setq j (1- j)) 0)
		 (setq bjm (math-sub (math-mul (math-mul j two-over-x) bj)
				     bjp)
		       bjp bj
		       bj bjm)
		 (if (> (nth 2 (math-abs-approx bj)) 10)
		     (setq bj (math-mul bj '(float 1 -10))
			   bjp (math-mul bjp '(float 1 -10))
			   ans (and ans (math-mul ans '(float 1 -10)))
			   sum (math-mul sum '(float 1 -10))))
		 (or (setq jsum (not jsum))
		     (setq sum (math-add sum bj)))
		 (if (= j v)
		     (setq ans bjp)))
	       (math-div ans (math-sub (math-mul 2 sum) bj))))))))