;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                                BESSEL-J
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun BESSEL-J (order x)
  "
;;;; BESSEL FUNCTIONS
;;; The Bessel functions J0(x) J1(x) . . .

;; Examples:

 (title \"Plot of Bessel Functions\")
 (plot
   (lambda(x) (bessel-j 0 x))
   (lambda(x) (bessel-j 1 x))
   (lambda(x) (bessel-j 2 x))
   (lambda(x) (bessel-j 3 x))
   0 12)
"
  (declare (float z2))
  (cond ((= order 0)				;NBS 9.4.1 eps < 5e-8
	 (cond ((<= (abs x) 3.0)
		(let ((z2 (expt (/ x 3.0) 2)))
		  (clmath-poly z2
			1.0000000  -2.2499997
			1.2656208  -0.3163866
			0.0444479  -0.0039444
			0.0002100)))
	       ((>= x 3.0)			;NBS 9.4.3
		(let* ((z      (/ 3.0 x))
		       (f0     (clmath-poly z		; eps 1.6E-8
				     0.79788456 -.00000077 -.00552740
				     -.00009512 0.00137237 -.00072805
				     0.00014476))
		       (theta0 (+ x
				  (clmath-poly z	; eps 7.0E-8
					-.78539816 -.04166397
					-.00003954 0.00262537
					-.00054125 -.00029333
					0.00013558))))
		  (declare (float z f0 theta0))
		  ;; Y0 is ... (/$ (*$ f0 (sin theta0)) (sqrt x))
		  (/ (* f0 (cos theta0)) (sqrt x))))
	       (t
		(BESSEL-J ORDER (- X)))))
	((= order 1)
	 (cond ((<= (abs x) 3.0)
		(let ((z2 (expt (/ x 3.0) 2)))
		  (* (* x
			(clmath-poly z2		;eps < x * 1.3e-8
			      0.50000000 -0.56249985
			      0.21093573 -0.03954289
			      0.00443319 -0.00031761
			      0.00001109)
			))))
	       ((>= x 3.0)			;NBS 9.4.3
		(let* ((z      (/ 3.0 x))
		       (f1     (clmath-poly z		; eps 4.0E-8
				     0.79788456 0.00000156 0.01659667
				     0.00017105 -.00249511 0.00113653
				     -.00020033))
		       (theta1 (+ x
				  (clmath-poly z	; eps 9.0E-8
					-2.35619449 +0.12499612
					+0.00005650 -0.00637879
					+0.00074348 +0.00079824
					-0.00029166))))
		  (declare (float z f1 theta1))
		  ;; Y1 is ... (/$ (*$ f1 (sin theta1)) (sqrt x))
		  (/ (* f1 (cos theta1)) (sqrt x))))
	       (t
		(- (bessel-j order (- x))))))
	((< order 0)
	 (setq order (- order))
	 (cond ((oddp order)
		(- (BESSEL-J order x)))
	       (t (BESSEL-J order x))))
	(t
	 (if (= x 0) 0
	   (- (* (/ (* 2 (1- order)) x)
		 (BESSEL-J (1- order) x))
	      (BESSEL-J (- order 2) x))
	   )
	 )
	)
  )