/*
Copyright (C) 2000, 2001 Barton Willis

Brief Description: 

Maxima code for evaluating  orthogonal polynomials listed in
Chapter 22 of Abramowitz and Stegun. 

Author:  

Barton Willis, willisb@unk.edu
Department of Mathematics and Statistics
University of Nebraska at Kearney
Kearney NE 68849

License: 

GPL;  see COPYING for details.  You  may redistribute unmodified
copies of the specfun package; you may not redistribute
modified versions of specfun.

Maxima version required: 

specfun has been tested under gcl / maxima 5.4  for Linux (Intel), 
gcl / maxima 5.5 for  Windows NT 4.0,  Macsyma 2.2 for Windows 
NT 4.0,  and Macsyma 422 for Linux (Intel).

*/

/*  Ported to maxima 5.4 on 11 Sept 2000 and 17 March .  I changed:

  1. added a floor function,
  
  2. changed some (...)!! into genfact statements; Macsyma applies 
      ceil to the second argument of genfact, Maxima doesn't.  
      Thus (-1)!! = 1 in Macsyma and it is an error in Maxima. Thus
      (2*n-1)!! is replaced by genfact(2*n-1,n,2) and 
      (2*n+1)!! is replaced by genfact(2*n+1,n+1,2).
*/

/*  Code for Jacobi polynomials, ultraspherical polynomials, 
    associated Legendre polynomials of the first kind,
    Chebyshev polynomials of first and second kinds,
    Laguerre polynomials, generalized Laguerre polynomials,  
    Legendre polynomials, (both p and q) Hermite polynomials, 
    spherical Bessel functions of first and second kinds, spherical 
    Hankel functions of the first and second kinds, and spherical 
    harmonics.

    With  few exceptions,  I use the notation and conventions of 
    Abramowitz and Stegun (A&S);  page numbers refer to the tenth 
    printing, December 1972.   I also refer to Gradshteyn  and 
    Ryzhik (G&S), the 1980 corrected and enlarged edition and
    Eugen Mertzbacher "Quantum Mechanics"  2nd edition, 1970.
    
    This code was developed and tested using Macsyma 422 under
    (Intel) RedHat Linux 6.1. It was then ported to Maxima 5.4 also
    under (Intel) RedHat Linux 6.1.  
    
    When convenient, functions are computed as special cases of
    the Jacobi polynomials. For the most part, formulae are written 
    as they are in  A&S; tinkering with the formulae can increase speed 
    at the expense  of greater risk  of errors.    
    
    Special values for the arguments are not computed any differently 
    than generic arguments; I don't want to have worry whether 
    normalizations  are consistent with the function definitions.

   This code is primarily intended for symbolic computation. It is hoped
   that it gives accurate floating point results as well; however, I don't
   make any claims that the algorithms are well suited for floating point
   evaluation.   Most polynomials are returned in the form
          b (a0 + a1 * f + ... + an f^n),  
    where f is a  polynomial of degree two or less and b, a0, ..., an  are 
    constants.   Some day I may  experiment with a new version of 
    specfun that is intended for floating point numbers.

    Domain errors are sometimes, but not always, trapped. Generally,
    instead of a domain error, an unevaluated function is returned. 
    Often, something else (usually the gamma function)  issues an error.
    
    When one or more arguments are floating point numbers, the return
    type is supposed to be given by the widest float in the  argument 
    list.  Presumably, all conversions are handled automatically; however,
    When  an argument is a complex float, this doesn't work.  Using
    coerce_return_type should fix this problem.
    
    I'd like to locally  set  float2bf to true to eliminate warnings about 
    conversion to  bfloat's; however, the code
    
   ... :=  block([..., float2bf : true], ... 
   
   sometimes gives a "Signal 11 caught."  (Macsyma 422) error 
   message when compiled. The translated Lisp code
   seems to work okay when float2bf is locally set to true.
   For now, I'll globally set float2bf to true. 
   
   The file "test_specfun.mc" tests this code.
*/

put('specfun,110,'version)$

/* Macsyma, but not maxima, has a floor function.  Here is a
floor function that is sufficient.  

This function returns floor(m,n).
*/

floor_rat(m,n) := block([ ],
  mode_declare([m, n, function(floor_rat)], fixnum),
  first(divide(m,n))
)$

eval_when([load, batch, batchload],
     float2bf : true
);

/* True iff x is a float or a big float.
*/

float_bfloatp(x) := block([ ],
     floatnump(x) or bfloatp(x)
);

/* True iff x is a complex float. The real and imaginary parts can
    either be floats or big floats.
*/

complex_floatp(x) := block([ ],
    if  atom(x) then (
        false
    ) else (
        constantp(x) and float_bfloatp(realpart(x)) and float_bfloatp(imagpart(x))
    )       
);

/* True iff x is a float, a big float, or a complex float.
*/

floating_pointp(x) := block([  ],
     float_bfloatp(x) or complex_floatp(x)
);
	      
/*  Coerce the return type y to match the input type x.
*/

coerce_return_type(x, y) := block([ratprint : false],
    if floatnump(x)  then (
         float( y )
    ) else if bfloatp(x) then (
         bfloat( y )
    ) else if complex_floatp (x) then (
         expand (y)   
   ) else (
         y
   )
);    	

	 
/* Compute the Jacobi polynomial using the table on page 789 of
A&S.  When a <= -1 or b <= -1, the polynomials are defined, but
they don't form an orthogonal set.  If you need the jacobi polynomials
for a <= -1 or b <= -1, compute rat(jacobi_p(n,a,b,x)) and substitute
for a and b.  It works.  

When n is an integer and a, b, and x are floats (but not bfloats),
jacobi_p calls float modedeclared function jacobi_pf.  This speeds
the evaluation of the Jacobi polynomials for floats.

If featurep(n, 'integer)  = true, return a sum representation 
(See G&R 8.960 page 1035).  To avoid problems caused 0^0 terms 
in a sum, I removed the m=0 and m=n terms  from the summation.  This 
is inelegant, but it avoids the bogus results like the following:

(c1) e : sum(x^k, k,0,n)$
(c2) ev(e, x=0);
(d2)              0

*/
     
jacobi_p(n,a,b,x) := block([sofar, w, abn, m, s_args],
   mode_declare(m,fixnum, [n, a, b, x, abn, sofar, w, function(jacobi_p)],any, s_args, list),
   if integerp(n) and n > -1 and (not integerp(a) or a > -1) and (not integerp(b) or b > -1) then (
      if apply ("and", map ('constantp, [a,b,x])) and apply ("and",  map ('floatnump,[a,b,x])) then (
  	     jacobi_pf (n, float (a), float (b), float(x))
       ) else (	     
           sofar : 1,
           w : 1,
           x : (x - 1) / 2,
           abn : a + b + n,
           for m : 1 thru n do (
               w : w * (n + 1 - m) * (abn + m) * x / (m * (a + m)),
	       sofar : sofar + w
           ),
           coerce_return_type (x, binomial(n + a, n) * sofar) 
      )	   
   )  else (
       if featurep(n, 'integer) then (
             s_args : [binomial(n+a,'m) * binomial(n+b,n-'m) * (x-1)^(n-'m) * (x + 1)^'m,'m,1,n-1],
             (  binomial(n+a,0) * binomial(n+b,n) * (x-1)^n
	    + binomial(n+a,n) * binomial(n+b,0) * (x + 1)^n
	    + apply('sum, s_args)) / 2^n
       ) else (
            funmake ('jacobi_p,[n,a,b,x])
       )	    
   )
)$

/* jacobi_p: A&S 22.8.1  page 783.*/

gradef(jacobi_p(n,a,b,x),
       'diff('jacobi_p(n,a,b,x),n),
        'diff('jacobi_p(n,a,b,x),a),
        'diff('jacobi_p(n,a,b,x),b),
        (n*(a-b-(2*n+a+b)*x)*'jacobi_p(n,a,b,x)+2*(n+a)*(n+b)
            * 'jacobi_p(n-1, a,b,x)) / ((2*n+a+b)*(1-x^2)));

/* jacobi_pf is the essentially the same as jacobi_p except  a,b, and 
x are declared to be floats and n is declared to be a fixnum.  This
isn't intended to be a user function; a user should use jacobi_p instead.
*/

jacobi_pf(n,a,b,x) := block([sofar, w, abn, np1, m],
    mode_declare([n,m],fixnum, [a, b, x, abn, np1, sofar, w, function(jacobi_pf)],float),       
   
    sofar : 1.0,
    w : 1.0,
    x : (x - 1.0) / 2.0,
    abn : a + b + n,
    np1 : 1.0 + n,
    for m : 1 thru n do (
         w : w * (np1 - m) * (abn + m) * x / (m * (a + m)),
	 sofar : sofar + w
    ),
    binomial(n + a, n) * sofar	   
)$

/* Compute the ultraspherical polynomials using the jacobi polynomials;
    see A&S 22.5.27  page  778.
*/
    
ultraspherical(n,a,x) := block([ ],
   mode_declare([n,a,x, function(ultraspherical)], any),
   gamma(a + 1/2) * gamma(2*a + n) * jacobi_p(n, a - 1/2, a - 1/2,x) / ( gamma(2 * a) * gamma(a + n + 1/2))
);


/* For now, n & m must be integers.  See A & S 22.5.37 page 779, 
   A & S 8.6.6 (second equation) page 334, and A & S 8.2.5 page 333.
*/

assoc_legendre_p(n,m,x) := block([ c ],
   mode_declare([n, m, function(assoc_legendre_p),x, c], any),
   
   if integerp ( n ) and integerp ( m ) and n > -1  then (
       if abs ( m ) > n then (
            0
       ) else if m < 0 then (
            gamma(n - m + 1) * assoc_legendre_p(n,-m,x) / gamma(n + m + 1)
       ) else (
           if m = 0 then (
	       c : 1
	   ) else (
	       c : genfact (2 * m - 1, m - 1, 2)
	   ), 
           c * (-1)^m * (1 - x^2)^(m / 2) * ultraspherical (n - m, m + 1 / 2, x) 
       )	   
   ) else (
      funmake('assoc_legendre_p, [n,m,x])
   )
);

/* See A&S 8.5.4 page 334.
*/

gradef(assoc_legendre_p(n, m, x),
        'diff('assoc_legendre_p(m,m,x),n),
        'diff('alegendre_p(n,m,x),m),
         (n * x * 'assoc_legendre_p(n, m, x)-(n+m) * 'assoc_legendre_p(n-1, m, x)) / (x^2 - 1) 
); 

/* See A&S 8.6.19 page 334.
*/

legendre_q(n, x) := block([m, leg_p, sofar, ratprint : false],
    mode_declare(m, fixnum, leg_p, list, [n, x, function(legendre_q), sofar], any),
    if integerp(n) and n > -1  then (
	leg_p : [ ],
	for m : 0 thru n do (
	    leg_p : endcons(legendre_p(m,x), leg_p)
	),
	
	sofar : 0,
	for m : 1 thru n do (
	   sofar : sofar + part(leg_p, m) * part(leg_p, n-m+1) / m
	),
	
        part(leg_p, n+1) * log((1+x)/(1-x)) / 2 - sofar
    ) else (
       funmake('legendre_q, [n,x])
    )    	
)$  

/* See A&S 8.6.7 and 8.2.6 pages 333 and 334.   I chose the 
definition that is real valued on (-1,1).  

For negative m, A&S 8.2.6 page 333 and  G&R 8.706 page 1000
disagree; the factor of exp(i m pi) in A&S 8.1.6 suggests to me that
A&S  8.2.6 is bogus.   As further evidence, Macsyma 2.4 seems to 
agree with G&R 8.706.   I'll use G&R.
*/

assoc_legendre_q(n,m,x) := block([ x%, f, ratprint : false],
   mode_declare([n,m,assoc_legendre_q, x], any),
   
   if integerp(n) and n > -1 and integerp(m) then (
      if m < 0 then (
           (-1)^m * (n+m)! * assoc_legendre_q(n,-m,x) / (n-m)!
      ) else (      
           f : ratsimp((-1)^m * (1 - x^2)^(m/2) * diff(legendre_q(n,x%),x%,m)),
	   subst(x, x%, f)
      )
  ) else (
      funmake('assoc_legendre_q, [n,m, x])
  )    
)$


/*  See A&S 22.5.31 page 778 and A&S 6.1.22 page 256.
*/
         
chebyshev_t(n,x) := block([ ],
   mode_declare([n, x, function(chebyshev_t)], any),
   
   if (integerp(n) and n > -1) or featurep(n, 'integer) then (
       n ! * 2^n * jacobi_p(n, -1/2, -1/2, x) / genfact(2*n - 1, n, 2)
   ) else (
       funmake('chebyshev_t,[n,x])
   )
)$      

/* See  A&S, 22.8.3 page 783. 
*/

gradef(chebyshev_t(n, x),
      'diff('chebyshev_t(n, x),n),
      (-n * x * 'chebyshev_t(n, x) + n * 'chebyshev_t(n-1, x))/(1-x^2)
);

/*  See A&S 22.5.32 page 778 and A&S 6.1.22 page 256.
*/

chebyshev_u(n,x) := block([ ],
   mode_declare([n,x, function(chebyshev_u)], any),
   
   if integerp(n) and n > -1 or featurep(n, 'integer) then (
       (n + 1) ! * 2^n * jacobi_p(n, 1/2, 1/2, x) / genfact(2*n + 1, n+1,2) 
   ) else (
       funmake('chebyshev_u,[n,x])
   )
)$      

/* See  A&S, 22.8.4 page 783. 
*/

gradef(chebyshev_u(n, x),
               'diff('chebyshev_u(n, x),n),
               (-n * x * 'chebyshev_u(n, x) + (n + 1) * 'chebyshev_u(n-1, x))/(1-x^2)
);
      
/* See A&S 22.5.16, page 778.
*/

laguerre(n, x) := block([ratprint : false],
   mode_declare([n, x, function(laguerre)], any),
   gen_laguerre(n, 0, x)
);

/* See table in A&S on page 789.
*/

gen_laguerre(n, a, x) := block([sofar, w, m, ratprint : false],
   mode_declare(m, fixnum, [n, x, a, sofar, w, function(gen_laguerre)], any),
   if integerp(n) and n > -1 then (
     sofar : 1,
     w : 1,
     for m : 1 thru n do (
         w : -w * x * (n - m + 1) / (m * (m + a)),
	 sofar : sofar + w
     ),
     coerce_return_type(x, binomial(n + a, n) * sofar)
  ) else (
     funmake('gen_laguerre,[n,a, x])   
  )
)$

/* See A&S 22.8.6 page 783.
*/

gradef(gen_laguerre(n, a, x),
      'diff('gen_laguerre(n, a, x),n),
      'diff(gen_laguerre(n, a, x),a),
      (n * gen_laguerre(n, a, x) - (n + a) * gen_laguerre(n-1, a, x)) / x
);
      

/*  See A&S 22.5.35 page 779.
*/

legendre_p(n,x) := block([   ],
   mode_declare([n, x, function(legendre_p)], any),
   if integerp(n) and n > -1 or featurep(n, 'integer) then (
         jacobi_p(n,0,0,x)
   ) else (
       funmake('legendre_p, [n,x])
   )
);

gradef(legendre_p (n,x),
     'diff('legendre_p (n,x),n),
      (-n * x * legendre_p (n,x) + n * legendre_p (n-1,x) ) / (1 - x^2)
);
   
 
/* See A&S 22.5.40 and 22.5.41, page 779.
*/

hermite(n,x) := block([ ],
   mode_declare([n, x, function(hermite)], any),
   
   if integerp(n) then (
      if evenp(n) then (
	 n : floor_rat(n,2),
	 coerce_return_type(x, (-1)^n * 4^n * n! * gen_laguerre(n, -1/2, x^2))
      ) else (
	 n : floor_rat(n, 2),
	 coerce_return_type(x, (-1)^n * 2^(2*n+1) * n! * x * gen_laguerre(n, 1/2, x^2))
      )
   ) else (
      funmake('hermite,[n,x])
   )
)$            	  	 

/* See A&S 22.8.7 page 783.
*/

gradef(hermite(n,x),
    'diff('hermite(n,x),n),
    2 * n * hermite(n-1,x)
);
    
/* See A&S 10.1.17 page 439.
*/

spherical_hankel2(n,x) := block([sofar, w, m],
   mode_declare(m, fixnum, [n, x, sofar, w, function(spherical_hankel2)], any),
   if integerp(n)  then (
      if n > -1 then (
          sofar : 1,
          w : 1,
          for m : 0 thru n do (
               w : w * (n + m + 1) * (n - m) / (2 * %i * x *  (m + 1)),
	       sofar : sofar + w
          ),
          sofar : %i^(n+1) * exp(-%i * x) * sofar / x,
          coerce_return_type(x, sofar)
     ) else (
        %i * (-1)^n * spherical_hankel2(-n-1,x)          
     ) 	  
  ) else (
     funmake('spherical_hankel2, [n, x])   
  )
)$
   
/* See A&S 10.1.36 page 439. (Let m = 0 in 10.1.36.)
*/

spherical_hankel1(n,x) := (-1)^n * spherical_hankel2(n,-x)$

/* See A&S 10.1.9 page 437.
*/

p_fun(n,x) := block([sofar, w, n1, m ],
   mode_declare([n, n1, m], fixnum, [sofar, w, x, function(p_fun)], any),
   sofar : 1,
   w : 1,
   n1 : floor_rat(n,2) - 1,
   for m : 0 thru n1 do (
       w : -w * (n + 2 * m + 2) * (n + 2 * m + 1) * (n - 2 * m) * (n - 2 * m - 1) 
               / (4 * (2 * m + 1) * (2 * m + 2) *  x^2),
       sofar : sofar + w
     ),
     sofar
)$     

q_fun(n,x) := block([sofar, w, n1, m],
   mode_declare([n, n1, m], fixnum, [sofar, w, x, function(q_fun)], any),
   
   if n = 0 then (
      0
   ) else (  
      sofar : 1,
      w : 1,
      n1 : floor_rat(n-1,2) - 1,
      for m : 0 thru n1 do (
          w : -w * (n + 2 * m + 3) * (n + 2 * m + 2) * (n - 2*m - 1) * (n - 2*m - 2)
	       / (4 * (2 * m + 3) * (2 * m + 2) * x^2),
         sofar : sofar + w
      ),
      (n + 1) ! * sofar /(2 * x * (n-1)!)
   )         
)$   

/* See A&S 10.1.8 page 437 and A&S 10.1.15 page 439.
*/
   
spherical_bessel_j(n,x) := block([  ],
   mode_declare([n, x, function(spherical_bessel_j)], any),
   if integerp(n)  then (
      if n > -1 then (
          (p_fun(n,x) * sin(x - n * %pi /2) + q_fun(n, x)  * cos(x - n * %pi / 2)) / x
      ) else (
         (-1)^n * spherical_bessel_y(-n-1,x)
      )	  
  ) else (
     funmake('spherical_bessel_j, [n, x])   
  )
)$


/* See A&S 10.1.9 page 437 and 10.1.15 page 439.
*/

spherical_bessel_y(n,x) := block([ ],
   mode_declare([n, x, function(spherical_bessel_y)], any),
   if integerp(n)  then (
       if n > -1 then (
            (-1)^(n + 1) * (p_fun(n,x) * cos(x + n * %pi /2) - q_fun(n, x)  * sin(x + n * %pi / 2)) / x
       ) else (
             (-1)^(n+1) * spherical_bessel_j(-n-1,x) 
	)	    
  ) else (
     funmake('spherical_bessel_y, [n, x])   
  )
)$

/* See A&S 10.1.19 page 439.
*/

gradef('spherical_hankel2(n,x),
       'diff('spherical_hankel2(n,x),n),
       (n * spherical_hankel2(n-1,x) - (n+1) * spherical_hankel2(n+1,x))/(2*n+1)
);
       
gradef('spherical_hankel1(n,x),
       'diff('spherical_hankel1(n,x),n),
       (n * spherical_hankel1(n-1,x) - (n+1) * spherical_hankel1(n+1,x))/(2*n+1)
);

gradef('spherical_bessel_j(n,x),
       'diff('spherical_bessel_j(n,x),n),
       (n * spherical_bessel_j(n-1,x) - (n+1) * spherical_bessel_j(n+1,x))/(2*n+1)
);

gradef('spherical_bessel_y(n,x),
       'diff('spherical_bessel_y(n,x),n),
       (n * spherical_bessel_y(n-1,x) - (n+1) * spherical_bessel_y(n+1,x))/(2*n+1)
);


/*  Compute P_n^m(cos(theta)).  See Mertzbacher, 9.59 page 184
and 9.64 page 185, and A & S 22.5.37 page 779.  assoc_leg_cos
lacks error checking; its should only be called by spherical_harmonic.
*/

assoc_leg_cos(n,m,x) := block([ c ],
   mode_declare([n, m], fixnum, [c, function(assoc_leg_cos),x], any),
   if m = 0 then (
	c : 1
    ) else (
	 c : (2 * m - 1)!!
    ),	       
    c * sin ( x )^m  * ultraspherical (n - m, m + 1 / 2, cos (x))
);

/* See Mertzbacher, "Quantum Mechanics"  9.64 page 185.  We 
need to be careful; the A&S associated Legendre function differs
from Mertzbacher.  Compare Mertzbacher 9.59 to A&S 8.6.6. 
*/

spherical_harmonic(n, m, theta, phi) := block([ ],
   mode_declare([n,m, theta, phi, function(spherical_harmonic)], any),

   if integerp(n) and integerp(m) and n > -1 and m <= n then (
       if m < 0 then (
            (-1)^m * spherical_harmonic(n,-m,theta,-phi)
        ) else (
            rat(exp(%i * m * phi)   * sqrt((2 * n + 1)  * (n - m)! / (4 *  %pi *  (n + m)!))
                   * assoc_leg_cos(n, m, theta))	  
        ) 		   
   ) else (
        funmake('spherical_harmonic,[n, m, theta, phi])
  )
)$  	

/* See Mertzbacher 9.66 page 185 and  9.72a-9.72b page 187.
*/

gradef(spherical_harmonic(n, m, theta, phi),
    'diff(spherical_harmonic(n, m, theta, phi),n),
    'diff(spherical_harmonic(n, m, theta, phi),m),
     -sqrt((n-m)*(n+m+1)) * exp(-%i * phi)  * spherical_harmonic(n, m+1, theta, phi)/2
     + sqrt((n+m)*(n-m+1)) * exp(%i * phi)  * spherical_harmonic(n, m-1, theta, phi)/2,
    %i * m * spherical_harmonic(n, m, theta, phi)
);