;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;                                                                    ;;;;;
;;;  Updated and modified by Raymond Toy, 2001                         ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "MAXIMA")
(macsyma-module ellipt)

;;;
;;; We use the definitions from Abramowitz and Stegun where our
;;; sn(u,m) is sn(u|m).  That is, the second arg is the parameter,
;;; instead of the modulus k or modular angle alpha.
;;;
;;; Note that m = k^2 and k = sin(alpha).
;;;

;;
;; Routines for computing the basic elliptic functions sn, cn, and dn.
;;
;;
;; A&S gives several methods for computing elliptic functions
;; including the AGM method (16.4) and ascending and descending Landen
;; transformations (16.12 and 16.14).  We use these latter because
;; they are actually quite fast, only requiring simple arithmetic and
;; square roots for the transformation until the last step.  The AGM
;; requires evaluation of several trignometric functions at each
;; stage.
;;
;; In addition, the Landen transformations are valid for all u and
;; Thus, we can compute the elliptic functions for complex u and
;; m. (The code below supports this, but we could make it run much
;; faster if we specialized it to for double-floats.  However, if we
;; do that, we won't be able to handle the cases where m < 0 or m > 1.
;; We'll have to handle these specially via A&S 16.10 and 16.11.)
;;
;; See A&S 16.12 and 16.14.

(flet ((ascending-transform (u m)
	 ;; A&S 16.14.1
	 ;;
	 ;; Take care in computing this transform.  For the case where
	 ;; m is complex, we should compute sqrt(mu1) first as
	 ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
	 ;; If not, we may choose the wrong branch when computing
	 ;; sqrt(mu1).
	 (let* ((root-m (sqrt m))
		(mu (/ (* 4 root-m)
		       (expt (1+ root-m) 2)))
		(root-mu1 (/ (- 1 root-m) (+ 1 root-m)))
		(mu1 (* root-mu1 root-mu1))
		(v (/ u (1+ root-mu1))))
	   (values v mu root-mu1)))
       (descending-transform (u m)
	 ;; Note: Don't calculate mu first, as given in 16.12.1.  We
	 ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
	 ;; then compute mu = sqrt(mu)^2.  If we calculate mu first,
	 ;; sqrt(mu) loses information when m or m1 is complex.
	 (let* ((root-m1 (sqrt (- 1 m)))
		(root-mu (/ (- 1 root-m1) (+ 1 root-m1)))
		(mu (* root-mu root-mu))
		(v (/ u (1+ root-mu))))
	   (values v mu root-mu))))
  (declaim (inline descending-transform ascending-transform))


  ;; Could use the descending transform, but some of my tests show
  ;; that it has problems without roundoff errors.
  (defun elliptic-dn-ascending (u m)
    (if (< (abs (- 1 m)) (* 4 long-float-epsilon))
	;; A&S 16.6.3
	(/ (cosh u))
	(multiple-value-bind (v mu root-mu1)
	    (ascending-transform u m)
	  ;; A&S 16.14.4
	  (let* ((new-dn (elliptic-dn-ascending v mu)))
	    (* (/ (- 1 root-mu1) mu)
	       (/ (+ root-mu1 (* new-dn new-dn))
		  new-dn))))))

  ;; Don't use the descending version because it requires cn, dn, and
  ;; sn.
  (defun elliptic-cn-ascending (u m)
    (if (< (abs (- 1 m)) (* 4 long-float-epsilon))
	;; A&S 16.6.2
	(/ (cosh u))
	(multiple-value-bind (v mu root-mu1)
	    (ascending-transform u m)
	  ;; A&S 16.14.3
	  (let* ((new-dn (elliptic-dn-ascending v mu)))
	    (* (/ (+ 1 root-mu1) mu)
	       (/ (- (* new-dn new-dn) root-mu1)
		  new-dn))))))

  ;; We don't use the ascending transform here because it requires
  ;; evaluating sn, cn, and dn.  The ascending transform only needs
  ;; sn.
  (defun elliptic-sn-descending (u m)
    ;; A&S 16.12.2
    (if (< (abs m) long-float-epsilon)
	(sin u)
	(multiple-value-bind (v mu root-mu)
	    (descending-transform u m)
	  (let* ((new-sn (elliptic-sn-descending v mu)))
	    (/ (* (1+ root-mu) new-sn)
	       (1+ (* root-mu new-sn new-sn)))))))
  #+nil
  (defun elliptic-sn-ascending (u m)
    (if (< (abs (- 1 m)) (* 4 long-float-epsilon))
	;; A&S 16.6.1
	(tanh u)
	(multiple-value-bind (v mu root-mu1)
	    (ascending-transform u m)
	  ;; A&S 16.14.2
	  (let* ((new-cn (elliptic-cn-ascending v mu))
		 (new-dn (elliptic-dn-ascending v mu))
		 (new-sn (elliptic-sn-ascending v mu)))
	    (/ (* (+ 1 root-mu1) new-sn new-cn)
	       new-dn)))))
)

(defun sn (u m)
  (if (and (realp u) (realp m))
      (realpart (elliptic-sn-descending u m))
      (elliptic-sn-descending u m)))

(defun cn (u m)
  (if (and (realp u) (realp m))
      (realpart (elliptic-cn-ascending u m))
      (elliptic-cn-ascending u m)))

(defun dn (u m)
  (if (and (realp u) (realp m))
      (realpart (elliptic-dn-ascending u m))
      (elliptic-dn-ascending u m)))



;;
;; How this works, I think.
;;
;; $sn is the user visible function.  This gets mapped to the symbol
;; %sn.  We put properties on this symbol so maxima can figure out
;; what to do with it.

;; Tell maxima how to simplify the functions %sn, etc.  This borrows
;; heavily from trigi.lisp.
(defprop %sn simp-%sn operators)
(defprop %cn simp-%cn operators)
(defprop %dn simp-%dn operators)
(defprop %asn simp-%asn operators)
(defprop %acn simp-%acn operators)
(defprop %adn simp-%adn operators)

;; Tell maxima what the derivatives are.  We currently only know how
;; to take derivatives of the first argument, not the second.  There
;; is a formula for this, but it's quite complicated so I'm not sure I
;; want to add that.
(defprop %sn ((u m)
	      ((mtimes) ((%cn) u m) ((%dn) u m))
	      ((%derivative) ((%sn simp) u m) m 1))
  grad)

(defprop %cn ((u m)
	      ((mtimes) -1 ((%sn) u m) ((%dn) u m))
	      ((%derivative) ((%cn simp) u m) m 1))
  grad)

(defprop %dn ((u m)
	      ((mtimes) -1 m ((%sn) u m) ((%cn) u m))
	      ((%derivative) ((%dn simp) u m) m 1))
  grad)

;; The inverse elliptic functions.
(defprop %asn ((x m)
	       ;; 1/sqrt(1-x^2)/sqrt(1-m*x^2)
	       ((MTIMES SIMP)
		((MEXPT SIMP) ((MPLUS SIMP) 1 ((MTIMES SIMP) -1 ((MEXPT SIMP) x 2)))
		 ((RAT SIMP) -1 2))
		((MEXPT SIMP) ((MPLUS SIMP) 1 ((MTIMES SIMP) -1 m ((MEXPT SIMP) x 2)))
		 ((RAT SIMP) -1 2)))
	       ((%derivative ((%asn) x m) m 1)))
  grad)

(defprop %acn ((x m)
	       ;; -1/sqrt(1-x^2)/sqrt(1-m*x^2)
	       ((mtimes simp) -1
		((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
		 ((rat simp) -1 2))
		((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
		 ((rat simp) -1 2)))
	       ((%derivative ((%acn) x m) m 1)))
  grad)

(defprop %adn ((x m)
	       ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
	       ((mtimes simp)
		((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
		 ((rat simp) -1 2))
		((mexpt simp) ((mplus simp) -1 m ((mexpt simp) x 2)) ((rat simp) -1 2)))
	       ((%derivative ((%adn) x m) m 1)))
  grad)

;; Define the actual functions for the user
(defmfun $sn (u m)
  (simplify (list '(%sn) (resimplify u) (resimplify m))))
(defmfun $cn (u m)
  (simplify (list '(%cn) (resimplify u) (resimplify m))))
(defmfun $dn (u m)
  (simplify (list '(%dn) (resimplify u) (resimplify m))))

(defmfun $asn (u m)
  (simplify (list '(%asn) (resimplify u) (resimplify m))))

(defmfun $acn (u m)
  (simplify (list '(%acn) (resimplify u) (resimplify m))))

(defmfun $adn (u m)
  (simplify (list '(%adn) (resimplify u) (resimplify m))))


;; A complex number looks like
;;
;; ((MPLUS SIMP) 0.70710678118654757 ((MTIMES SIMP) 0.70710678118654757 $%I))
;;
;; or
;;
;; ((MPLUS SIMP) 1 $%I))
;;
(defun complex-number-p (u)
  ;; Return non-NIL if U is a complex number (or number)
  (or (numberp u)
      (and (consp u)
	   (numberp (second u))
	   (or (and (consp (third u))
		    (numberp (second (third u)))
		    (eq (third (third u)) '$%i))
	       (and (eq (third u) '$%i))))))

(defun complexify (x)
  ;; Convert a Lisp number to a maxima number
  (if (realp x)
      x
      `((mplus simp) ,(realpart x)
	       ((mtimes simp) ,(imagpart x) $%i))))
;; Tell maxima how to simplify the functions
;;
;; FORM is list containing the actual expression.  I don't really know
;; what Y and Z contain.  Most of this modeled after SIMP-%SIN.
(defmfun simp-%sn (form y z)
  (declare (ignore y))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate sn
	   (sn (float u 1d0) (float m 1d0)))
	  #+nil
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   ;; For complex values.  Should we really do this?
	   (format t "real u = ~A~%" ($realpart u))
	   (format t "imag u = ~A~%" ($imagpart u))
	   (let ((result (sn (complex ($realpart u) ($imagpart u))
			     (complex ($realpart m) ($imagpart m)))))
	     (complexify result)))
	  ((zerop1 u)
	   ;; A&S 16.5.1
	   0)
	  ((zerop1 m)
	   ;; A&S 16.6.1
	   `(($sin) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.1
	   `(($tanh) ,u))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%sn) u m) form)))))

(defmfun simp-%cn (form y z)
  (declare (ignore y))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (cn (float u 1d0) (float m 1d0)))
	  #+nil
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (format t "real u = ~A~%" ($realpart u))
	   (format t "imag u = ~A~%" ($imagpart u))
	   (let ((result (cn (complex ($realpart u) ($imagpart u))
			     (complex ($realpart m) ($imagpart m)))))
	     (complexify result)))
	  ((zerop1 u)
	   ;; A&S 16.5.1
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.2
	   `(($cos) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.2
	   `(($sech) ,u))
	  (t (eqtest (list '(%cn) u m) form)))))

(defmfun simp-%dn (form y z)
  (declare (ignore y))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (dn (float u 1d0) (float m 1d0)))
	  #+nil
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (format t "real u = ~A~%" ($realpart u))
	   (format t "imag u = ~A~%" ($imagpart u))
	   (let ((result (dn (complex ($realpart u) ($imagpart u))
			     (complex ($realpart m) ($imagpart m)))))
	     (complexify result)))
	  ((zerop1 u)
	   ;; A&S 16.5.1
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.3
	   1)
	  ((onep1 m)
	   ;; A&S 16.6.3
	   `(($sech) ,u))
	  (t (eqtest (list '(%dn) u m) form)))))

;; Should we simplify the inverse elliptic functions into the
;; appropriate incomplete elliptic integral?  I think we should leave
;; it, but perhaps allow some way to do that transformation if
;; desired.

(defmfun simp-%asn (form y z)
  (declare (ignore y))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate asn
	   ;;
	   ;; asn(x,m) = F(asin(x),m)
	   (elliptic-f (asin x) m))
	  ((zerop1 m)
	   ;; asn(x,0) = F(asin(x),0) = asin(x)
	   `((%ellipt_f) ((%asin) ,u) 0))
	  ((onep1 m)
	   ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
	   `((%ellipt_f) ((%asin) ,u) 1))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%asn) u m) form)))))

(defmfun simp-%acn (form y z)
  (declare (ignore y))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate acn
	   ;;
	   ;; acn(x,m) = F(acos(x),m)
	   (elliptic-f (acos x) m))
	  ((zerop1 m)
	   ;; asn(x,0) = F(acos(x),0) = acos(x)
	   `((%ellipt_f) ((%acos) ,u) 0))
	  ((onep1 m)
	   ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
	   `((%ellipt_f) ((%acos) ,u) 1))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%acn) u m) form)))))

;;
;; adn(u) = x.
;; u = dn(x) = sqrt(1-m*sn(x)^2)
;; sn(x)^2 = (1-u^2)/m
;; sn(x) = sqrt((1-u^2)/m)
;; x = asn(sqrt((1-u^2)/m))
;; x = adn(u) = asn(sqrt((1-u^2)/m))

(defmfun simp-%adn (form y z)
  (declare (ignore y))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate asn
	   ;;
	   (elliptic-f (asin (/ (* (sqrt (- 1 u)) (sqrt (+ 1 u)))
				(sqrt m)))
		       m))
	  ((onep1 m)
	   ;; x = dn(u,1) = sech(u).  so u = asech(x)
	   `((%asech) ,u))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%adn) u m) form)))))


;;;; Elliptic integrals

;; Carlson's elliptic integral of the first kind.
;;
;; ***PURPOSE  Compute the incomplete or complete elliptic integral of the
;;             1st kind.  For X, Y, and Z non-negative and at most one of
;;             them zero, RF(X,Y,Z) = Integral from zero to infinity of
;;                                 -1/2     -1/2     -1/2
;;                       (1/2)(t+X)    (t+Y)    (t+Z)    dt.
;;             If X, Y or Z is zero, the integral is complete.
;;
;;          Value of IER assigned by the DRF routine
;;
;;                   Value assigned         Error Message Printed
;;                   IER = 1                MIN(X,Y,Z) .LT. 0.0D0
;;                       = 2                MIN(X+Y,X+Z,Y+Z) .LT. LOLIM
;;                       = 3                MAX(X,Y,Z) .GT. UPLIM
;;
;;    Special double precision functions via DRF
;;
;;
;;
;;
;;                   Legendre form of ELLIPTIC INTEGRAL of 1st kind
;;
;;                   -----------------------------------------
;;
;;
;;
;;                                              2         2   2
;;                   F(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1)
;;
;;
;;                                   2
;;                   K(K) = DRF(0,1-K ,1)
;;
;;
;;                          PI/2     2   2      -1/2
;;                        = INT  (1-K SIN (PHI) )   D PHI
;;                           0
;;
;;
;;
;;                   Bulirsch form of ELLIPTIC INTEGRAL of 1st kind
;;
;;                   -----------------------------------------
;;
;;
;;                                           2 2    2
;;                   EL1(X,KC) = X DRF(1,1+KC X ,1+X )
;;
;;
;;                   Lemniscate constant A
;;
;;                   -----------------------------------------
;;
;;
;;                        1      4 -1/2
;;                   A = INT (1-S )    DS = DRF(0,1,2) = DRF(0,2,1)
;;                        0
;;
;;
;;
;;     -------------------------------------------------------------------
;;
;; ***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
;;                  elliptic integrals, ACM Transactions on Mathematical
;;                  Software 7, 3 (September 1981), pp. 398-403.
;;                B. C. Carlson, Computing elliptic integrals by
;;                  duplication, Numerische Mathematik 33, (1979),
;;                  pp. 1-16.
;;                B. C. Carlson, Elliptic integrals of the first kind,
;;                  SIAM Journal of Mathematical Analysis 8, (1977),
;;                 pp. 231-242.

(let ((errtol (expt (* 4 long-float-epsilon) 1/6))
      (uplim (/ most-positive-long-float 5))
      (lolim (* least-positive-normalized-long-float 5))
      (c1 (float 1/24 1L0))
      (c2 (float 3/44 1L0))
      (c3 (float 1/14 1L0)))
  (declare (long-float errtol c1 c2 c3))
  (defun drf (x y z)
    "Compute Carlson's incomplete or complete elliptic integral of the
first kind:

                   INF
                  /
                  [                     1
  RF(x, y, z) =   I    ----------------------------------- dt
                  ]    SQRT(x + t) SQRT(y + t) SQRT(z + t)
                  /
                   0


where x >= 0, y >= 0, z >=0, and at most one of x, y, z is zero.
"
    (declare (long-float x y z))
    ;; Check validity of input
    (assert (and (>= x 0) (>= y 0) (>= z 0)
		 (plusp (+ x y)) (plusp (+ x z)) (plusp (+ y z))))
    (assert (<= (max x y z) uplim))
    (assert (>= (min (+ x y) (+ x z) (+ y z)) lolim))
    (locally 
	(declare (type (long-float 0L0) x y)
		 (type (long-float (0L0)) z)
		 (optimize (speed 3)))
      (loop
	  (let* ((mu (/ (+ x y z) 3))
		 (x-dev (- 2 (/ (+ mu x) mu)))
		 (y-dev (- 2 (/ (+ mu y) mu)))
		 (z-dev (- 2 (/ (+ mu z) mu))))
	    (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
	      (let ((e2 (- (* x-dev y-dev) (* z-dev z-dev)))
		    (e3 (* x-dev y-dev z-dev)))
		(return (/ (+ 1
			      (* e2 (- (* c1 e2)
				       1/10
				       (* c2 e3)))
			      (* c3 e3))
			   (sqrt mu)))))
	    (let* ((x-root (sqrt x))
		   (y-root (sqrt y))
		   (z-root (sqrt z))
		   (lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
	      (setf x (* (+ x lam) 1/4))
	      (setf y (* (+ y lam) 1/4))
	      (setf z (* (+ z lam) 1/4))))))))

;; Elliptic integral of the first kind (Legendre's form):
;;
;;
;;      phi
;;     /
;;     [             1
;;     I    ------------------- ds
;;     ]                  2
;;     /    SQRT(1 - m SIN (s))
;;     0

(defun elliptic-f (phi m)
  (declare (double-float phi m))
  (cond ((> m 1)
	 ;; A&S 17.4.15
	 (/ (elliptic-f (asin (* (sqrt m) (sin phi))) (/ m))))
	((< m 0)
	 ;; A&S 17.4.17
	 (let* ((m (- m))
		(m+1 (+ 1 m))
		(root (sqrt m+1))
		(m/m+1 (/ m m+1)))
	   (- (/ (elliptic-f (/ pi 2) m/m+1)
		 root)
	      (/ (elliptic-f (- (/ pi 2) phi) m/m+1)
		 root))))
	((= m 0)
	 ;; A&S 17.4.19
	 phi)
	((= m 1)
	 ;; A&S 17.4.21
	 (log (tan (+ (/ phi 2) (/ pi 2)))))
	(t
	 (let ((sin-phi (sin phi))
	       (cos-phi (cos phi))
	       (k (sqrt m)))
	   (* sin-phi
	      (drf (* cos-phi cos-phi)
		   (* (- 1 (* k sin-phi))
		      (+ 1 (* k sin-phi)))
		   1d0))))))
;;
;; Carlsons' elliptic integral of the second kind.
;;
;;   1.     DRD
;;          Evaluate an INCOMPLETE (or COMPLETE) ELLIPTIC INTEGRAL
;;          of the second kind
;;          Standard FORTRAN function routine
;;          Double precision version
;;          The routine calculates an approximation result to
;;          DRD(X,Y,Z) = Integral from zero to infinity of
;;                              -1/2     -1/2     -3/2
;;                    (3/2)(t+X)    (t+Y)    (t+Z)    dt,
;;          where X and Y are nonnegative, X + Y is positive, and Z is
;;          positive.  If X or Y is zero, the integral is COMPLETE.

;;    -------------------------------------------------------------------
;;
;;
;;   Special double precision functions via DRD and DRF
;;
;;
;;                  Legendre form of ELLIPTIC INTEGRAL of 2nd kind
;;
;;                  -----------------------------------------
;;
;;
;;                                             2         2   2
;;                  E(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1) -
;;
;;                     2      3             2         2   2
;;                  -(K/3) SIN (PHI) DRD(COS (PHI),1-K SIN (PHI),1)
;;
;;
;;                                  2        2            2
;;                  E(K) = DRF(0,1-K ,1) - (K/3) DRD(0,1-K ,1)
;;
;;                         PI/2     2   2      1/2
;;                       = INT  (1-K SIN (PHI) )  D PHI
;;                          0
;;
;;                  Bulirsch form of ELLIPTIC INTEGRAL of 2nd kind
;;
;;                  -----------------------------------------
;;
;;                                               2 2    2
;;                  EL2(X,KC,A,B) = AX DRF(1,1+KC X ,1+X ) +
;;
;;                                              3          2 2    2
;;                                 +(1/3)(B-A) X DRD(1,1+KC X ,1+X )
;;
;;
;;
;;
;;                  Legendre form of alternative ELLIPTIC INTEGRAL
;;                  of 2nd kind
;;
;;                  -----------------------------------------
;;
;;
;;
;;                            Q     2       2   2  -1/2
;;                  D(Q,K) = INT SIN P  (1-K SIN P)     DP
;;                            0
;;
;;
;;
;;                                     3          2     2   2
;;                  D(Q,K) = (1/3) (SIN Q) DRD(COS Q,1-K SIN Q,1)
;;
;;
;;
;;
;;                  Lemniscate constant  B
;;
;;                  -----------------------------------------
;;
;;
;;
;;
;;                       1    2    4 -1/2
;;                  B = INT  S (1-S )    DS
;;                       0
;;
;;
;;                  B = (1/3) DRD (0,2,1)
;;
;;
;;                  Heuman's LAMBDA function
;;
;;                  -----------------------------------------
;;
;;
;;
;;                  (PI/2) LAMBDA0(A,B) =
;;
;;                                    2                2
;;                 = SIN(B) (DRF(0,COS (A),1)-(1/3) SIN (A) *
;;
;;                            2               2         2       2
;;                  *DRD(0,COS (A),1)) DRF(COS (B),1-COS (A) SIN (B),1)
;;
;;                            2       3             2
;;                  -(1/3) COS (A) SIN (B) DRF(0,COS (A),1) *
;;
;;                           2         2       2
;;                   *DRD(COS (B),1-COS (A) SIN (B),1)
;;
;;
;;
;;                  Jacobi ZETA function
;;
;;                  -----------------------------------------
;;
;;                             2                 2       2   2
;;                  Z(B,K) = (K/3) SIN(B) DRF(COS (B),1-K SIN (B),1)
;;
;;
;;                                       2             2
;;                             *DRD(0,1-K ,1)/DRF(0,1-K ,1)
;;
;;                               2       3           2       2   2
;;                            -(K /3) SIN (B) DRD(COS (B),1-K SIN (B),1)
;;
;;

(let ((errtol (expt (/ long-float-epsilon 3) 1/6))
      (c1 (float 3/14 1L0))
      (c2 (float 1/6 1L0))
      (c3 (float 9/22 1L0))
      (c4 (float 3/26 1L0)))
  (declare (long-float errtol c1 c2 c3 c4))
  (defun drd (x y z)
    (declare (real x y z))
    ;; Check validity of input

    (assert (and (>= x 0) (>= y 0) (>= z 0) (plusp (+ x y)))
	    (x y z))
	    
    (let ((x (float x 1L0))
	  (y (float y 1L0))
	  (z (float z 1L0))
	  (sigma 0L0)
	  (power4 1L0))
      (declare (type (long-float 0L0) x y power4 sigma)
	       (type (long-float (0L0)) z)
	       (optimize (speed 3)))
      (loop
	  (let* ((mu (* 1/5 (+ x y (* 3 z))))
		 (x-dev (/ (- mu x) mu))
		 (y-dev (/ (- mu y) mu))
		 (z-dev (/ (- mu z) mu)))
	    (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
	      (let* ((ea (* x-dev y-dev))
		     (eb (* z-dev z-dev))
		     (ec (- ea eb))
		     (ed (- ea (* 6 eb)))
		     (ef (+ ed ec ec))
		     (s1 (* ed (+ (- c1)
				  (* 1/4 c3 ed)
				  (* -3/2 c4 z-dev ef))))
		     (s2 (* z-dev (+ (* c2 ef)
				     (* z-dev (+ (* (- c3) ec)
						 (* z-dev c4 ea)))))))
		(return (+ (* 3 sigma)
			   (/ (* power4 (+ 1 s1 s2))
			      (* mu (sqrt mu)))))))
	    (let* ((x-root (sqrt x))
		   (y-root (sqrt y))
		   (z-root (sqrt z))
		   (lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
	      (incf sigma (/ power4 z-root (+ z lam)))
	      (setf power4 (/ power4 4))
	      (setf x (/ (+ x lam) 4))
	      (setf y (/ (+ y lam) 4))
	      (setf z (/ (+ z lam) 4))))))))

;; Elliptic integral of the second kind (Legendre's form):
;;
;;
;;      phi
;;     /
;;     [                  2
;;     I    SQRT(1 - m SIN (s)) ds
;;     ]
;;     /
;;      0

(defun elliptic-e (phi m)
  (declare (double-float phi m))
  (cond ((= m 0)
	 ;; A&S 17.4.23
	 phi)
	((= m 1)
	 ;; A&S 17.4.25
	 (sin phi))
	(t
	 (let* ((sin-phi (sin phi))
		(cos-phi (cos phi))
		(k (sqrt m))
		(y (* (- 1 (* k sin-phi))
		      (+ 1 (* k sin-phi)))))
	   (- (* sin-phi
		 (drf (* cos-phi cos-phi) y 1d0))
	      (* (/ m 3)
		 (expt sin-phi 3)
		 (drd (* cos-phi cos-phi) y 1d0)))))))


;; Define the elliptic integrals for maxima
;;
;; We use the definitions given in A&S 17.2.6 and 17.2.8.  In particular:
;;
;;                 phi
;;                /
;;                [             1
;; F(phi|m)  =    I    ------------------- ds
;;                ]                  2
;;                /    SQRT(1 - m SIN (s))
;;                 0
;;
;; and
;;
;;              phi
;;             /
;;             [                  2
;; E(phi|m) =  I    SQRT(1 - m SIN (s)) ds
;;             ]
;;             /
;;              0
;;
;; That is, we do not use the modular angle, alpha, as the second arg;
;; the parameter m = sin(alpha)^2 is used.
;;


(defprop %ellipt_f simp-%ellipt_f operators)
(defprop %ellipt_e simp-%ellipt_e operators)

;; The derivatives wrt to k (m) are complicated, so I'm not going to
;; do this.
(defprop %ellipt_f ((phi m)
		    ;; 1/sqrt(1-m*sin(phi)^2)
		    ((mexpt simp)
		     ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
		     ((rat simp) -1 2))
		    ((%derivative) ((%ellipt_f simp) phi m) m 1))
  grad)

(defprop %ellipt_e ((phi m)
		    ;; (1-m*sin(phi)^2)
		    ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
		    ((%derivative) ((%ellipt_e simp) phi m) m 1))
  grad)
		    
(defmfun $ellipt_f (phi m)
  (simplify (list '(%ellipt_f) (resimplify phi) (resimplify m))))
(defmfun $ellipt_e (phi m)
  (simplify (list '(%ellipt_e) (resimplify phi) (resimplify m))))

(defmfun simp-%ellipt_f (form y z)
  (twoargcheck form)
  (let ((phi (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp phi) (floatp m))
	       (and $numer (numberp phi) (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-f (float phi 1d0) (float m 1d0)))
	  ((zerop1 m)
	   ;; A&S 17.4.19
	   phi)
	  ((onep1 m)
	   ;; A&S 17.4.21.  Let's pick the log tan form.
	   `((%log) ((%tan)
		     ((mplus) ((mtimes simp) $%pi ((rat simp) 1 4))
		      ((mtimes simp) ((rat simp) 1 2) ,phi)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%ellipt_f) phi m) form)))))

(defmfun simp-%ellipt_e (form y z)
  (twoargcheck form)
  (let ((phi (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp phi) (floatp m))
	       (and $numer (numberp phi) (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-e (float phi 1d0) (float m 1d0)))
	  ((zerop1 m)
	   ;; A&S 17.4.23
	   phi)
	  ((onep1 m)
	   ;; A&S 17.4.25
	   `((%sin) ,phi))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%ellipt_e) phi m) form)))))