;;; Lisplab, level0-specunc.lisp
;;; Special functions for numeric arguments. Using Slatec.
;;; 

;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License along
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

(in-package :lisplab)

(defun neg-integer-p (x)
  (and (< x 0) (= x (truncate x))))

(defmethod .besj (n (x number))
  "f2cl slatec based implementation"
  ;; Bessel J function, for n >=0, real and complex numbers. 
  ;; TODO: what about negaive n and complex n?
  (if (neg-integer-p n)
      (* (expt -1 n) (.besj (- n) x))
      (typecase x 
	(complex (let ((rx (to-df (realpart x)))
		       (cx (to-df (imagpart x)))
		       (ry (make-dvec 1))
		       (cy (make-dvec 1)))
		   (slatec:zbesj rx cx (to-df n) 1 1 ry cy 0 0)
		   (complex (aref ry 0) (aref cy 0))))
	(t (let ((x (to-df x)))
	     (case n 
	       (0 (slatec:dbesj0 x))
	       (1 (slatec:dbesj1 x))
	       (t (let ((y (make-dvec 1)))
		    (slatec:dbesj x (to-df n) 1 y 0)
		    (aref y 0)))))))))

(defmethod .besy (n (x number))
  "f2cl slatec based implementation"
  ;; Bessel Y function (Neumann function), for n >=0, x>0, real and complex numbers. 
  ;; TODO: what about negaive n, negative x and complex n?
   (typecase x 
    (complex (let ((rx (to-df (realpart x)))
		   (cx (to-df (imagpart x)))
		   (ry (make-dvec 1))
		   (cy (make-dvec 1))
		   (rw (make-dvec 1))
		   (cw (make-dvec 1)))
	       (slatec:zbesy rx cx (to-df n) 1 1 ry cy 0 rw cw 0)
	       (complex (aref ry 0) (aref cy 0))))
    (t (let ((x (to-df x)))
	 (case n 
	   (0 (slatec:dbesy0 x))
	   (1 (slatec:dbesy1 x))
	   (t (let ((y (make-dvec 1)))
		(slatec:dbesy x (to-df n) 1 y)
		(aref y 0))))))))

(defmethod .besi (n (x number))
  "f2cl slatec based implementation"
  ;; Bessel I function, for n >=0, x>0, real and complex numbers. 
  ;; TODO: what about negaive n, negative x and complex n?
  (typecase x 
    (complex (let ((rx (to-df (realpart x)))
		   (cx (to-df (imagpart x)))
		   (ry (make-dvec 1))
		   (cy (make-dvec 1)))
	       (slatec:zbesi  rx cx (to-df n) 1 1 ry cy 0 0)
	       (complex (aref ry 0) (aref cy 0))))
    (t (let ((x (to-df x)))
	 (case n 
	   (0 (slatec:dbesi0 x))
	   (1 (slatec:dbesi1 x))
	   (t (let ((y (make-dvec 1)))
		(slatec:dbesi x (to-df n) 1 1 y 0)
		(aref y 0))))))))

(defmethod .besk (n (x number))
  "f2cl slatec based implementation"
  ;; Bessel K function, for n >=0, x>0, real and complex numbers. 
  ;; TODO: what about negaive n, negative x and complex n?
  (typecase x 
    (complex (let ((rx (to-df (realpart x)))
		   (cx (to-df (imagpart x)))
		   (ry (make-dvec 1))
		   (cy (make-dvec 1)))
	       (slatec:zbesk rx cx (to-df n) 1 1 ry cy 0 0)
	       (complex (aref ry 0) (aref cy 0))))
    (t (let ((x (to-df x)))
	 (case n 
	   (0 (slatec:dbesk0 x))
	   (1 (slatec:dbesk1 x))
	   (t (let ((y (make-dvec 1)))
		(slatec:dbesk x (to-df n) 1 1 y 0)
		(aref y 0))))))))

(defmethod .besh1 (n (x number))
  "f2cl slatec based implementation"
  ;; Bessel H1 function, for n >=0, x>0, real and complex numbers. 
  ;; TODO: what about negaive n, negative x and complex n?
  (let ((rx (to-df (realpart x)))
	(cx (to-df (imagpart x)))
	(ry (make-dvec 1))
	(cy (make-dvec 1)))
    (slatec:zbesh rx cx (to-df n) 1 1 1 ry cy 0 0 )
    (complex (aref ry 0) (aref cy 0))))

(defmethod .besh2 (n (x number))
  "f2cl slatec based implementation"
  ;; Bessel H2 function, for n >=0, x>0, real and complex numbers. 
  ;; TODO: what about negaive n, negative x and complex n?
  (let ((rx (to-df (realpart x)))
	(cx (to-df (imagpart x)))
	(ry (make-dvec 1))
	(cy (make-dvec 1)))
    (slatec:zbesh rx cx (to-df n) 1 2 1 ry cy 0 0 )
    (complex (aref ry 0) (aref cy 0))))

(defmethod .erf ((x number))
  "f2cl slatec based implementation"
  (let ((x (to-df x)))
    (slatec:derf x)))

(defmethod .erfc ((x number))
  "f2cl slatec based implementation"
  (let ((x (to-df x)))
    (slatec:derfc x)))

(defmethod .gamma ((x number))
  "f2cl slatec based implementation"
  (let ((x (to-df x)))
    (slatec:dgamma x)))



#|




(defgeneric .erfc (x)) 

(defgeneric .gamma (x)) 

(defmethod .Ai (x) 
  (slatec:dai x))


|#