The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
pp_addpm({At=>Top},<<'EOD');
=head1 NAME

PDL::GSLSF::LEGENDRE - PDL interface to GSL Special Functions

=head1 DESCRIPTION

This is an interface to the Special Function package present in the GNU Scientific Library. 

=head1 SYNOPSIS

=cut


EOD

# PP interface to GSL

pp_addhdr('
#include <gsl/gsl_sf.h>
#include "../gslerr.h"
');

pp_def('gsl_sf_legendre_Pl',
       GenericTypes => [D],
       OtherPars =>'int l',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_legendre_Pl_e,($COMP(l),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'P_l(x)'
      );

pp_def('gsl_sf_legendre_Pl_array',
       GenericTypes => [D],
       OtherPars =>'int l=>num',
       Pars=>'double x(); double [o]y(num)',
       Code =>'
GSLERR(gsl_sf_legendre_Pl_array,($COMP(l)-1,$x(),$P(y)))
',
       Doc =>'P_l(x) from 0 to n-1.'
      );

pp_def('gsl_sf_legendre_Ql',
       GenericTypes => [D],
       OtherPars =>'int l',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_legendre_Ql_e,($COMP(l),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Q_l(x)'
      );

pp_def('gsl_sf_legendre_Plm',
       GenericTypes => [D],
       OtherPars =>'int l; int m',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_legendre_Plm_e,($COMP(l),$COMP(m),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'P_lm(x)'
      );

my $v = `gsl-config --version`;
if (defined($v) && $v>=2.0){

    pp_def('gsl_sf_legendre_array',
	   GenericTypes => [D],
	   OtherPars => 'char norm;  int lmax; int csphase',
	   Pars => 'double x(); double [o]y(n); double [t]work(wn)',
	   RedoDimsCode => '
$SIZE(wn)=gsl_sf_legendre_array_n($COMP(lmax));
$SIZE(n)=$COMP(lmax)*($COMP(lmax)+1)/2+$COMP(lmax)+1;
',
	   Code => <<'EOC',
           int i;
           if($x()<-1||$x()>1) barf("The input to gsl_sf_legendre_array must be abs(x)<=1, and you input %f. Try normalizing your input.",$x());
	   switch ($COMP(norm)){
	       case 'P' :
		   GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_NONE, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
	       break;
	       case 'S' :
		   GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_SCHMIDT, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
	       break;
	       case 'Y' :
		   GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_SPHARM, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
	       break;
	       case 'N' :
		   GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_FULL, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
	       break;
               default :
                ;
           }
           for(i=0; i<$SIZE(n); i++) {
              $y(n=>i) = $work(wn=>i);
           }

EOC

	   HandleBad => 1,
	   BadCode => <<'EOBC',
	   int i;
	   if ( $ISBAD( x() ) ) {
              loop(n) %{ $SETBAD ( y() ); %}
	   } else {
             if($x()<-1||$x()>1) barf("The input to gsl_sf_legendre_array must be abs(x)<=1, and you input %f",$x());
	     switch ($COMP(norm)) {
		   case 'P' :
		       GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_NONE, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
		   break;
		   case 'S' :
		       GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_SCHMIDT, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
		   break;
		   case 'Y' :
		       GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_SPHARM, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
		   break;
		   case 'N' :
		       GSLERR(gsl_sf_legendre_array_e,(GSL_SF_LEGENDRE_FULL, $COMP(lmax), $x(), $COMP(csphase), $P(work)));
		   break;
		   default :
		       ;
	       }
             for(i=0; i<$SIZE(n); i++) {
                $y(n=>i) = $work(wn=>i);
             }
	   }

EOBC



	   Doc => <<'EOD',
=for ref

Calculate all normalized associated Legendre polynomials.

=for usage

$Plm = gsl_sf_legendre_array($x,'P',4,-1);

The calculation is done for degree 0 <= l <= lmax and order 0 <= m <= l on the range abs(x)<=1.

The parameter norm should be:

=over 3

=item 'P' for unnormalized associated Legendre polynomials P_l^m(x),

=item 'S' for Schmidt semi-normalized associated Legendre polynomials S_l^m(x),

=item 'Y' for spherical harmonic associated Legendre polynomials Y_l^m(x), or

=item 'N' for fully normalized associated Legendre polynomials N_l^m(x).

=back

lmax is the maximum degree l.
csphase should be (-1) to INCLUDE the Condon-Shortley phase factor (-1)^m, or (+1) to EXCLUDE it.

See L<gsl_sf_legendre_array_index> to get the value of C<l> and C<m> in the returned vector.

EOD

	);

    pp_def('gsl_sf_legendre_array_index',
	   OtherPars => 'int lmax',
	   Pars => 'int [o]l(n); int [o]m(n)',
	   RedoDimsCode => '$SIZE(n)=$COMP(lmax)*($COMP(lmax)+1)/2+$COMP(lmax)+1;',
	   Code => q/
           int ell, em, index;
	   for (ell=0; ell<=$COMP(lmax); ell++){
	       for (em=0; em<=ell; em++){
		   index = gsl_sf_legendre_array_index(ell,em);
		   $l(n=>index)=ell;
		   $m(n=>index)=em;
	       }
	   }/,
	   Doc =>'=for ref

Calculate the relation between gsl_sf_legendre_arrays index and l and m values.

=for usage
($l,$m) = gsl_sf_legendre_array_index($lmax);

Note that this function is called differently than the corresponding GSL function, to make it more useful for PDL: here you just input the maximum l (lmax) that was used in C<gsl_sf_legendre_array> and it calculates all l and m values.'
	);


} elsif (defined($v) && $v<2.0) {

pp_def('gsl_sf_legendre_Plm_array',
       GenericTypes => [D],
       OtherPars =>'int l=>num; int m',
       Pars=>'double x(); double [o]y(num)',
       Code =>'
GSLERR(gsl_sf_legendre_Plm_array,($COMP(l)-2+$COMP(m),$COMP(m),$x(),$P(y)))
',
       Doc =>'P_lm(x) for l from 0 to n-2+m.
gsl_sf_legendre_Plm_array has been deprecated in GSL version 2.0. It is included here for backwards compatability and may be removed in a future release.  New code should use L<gsl_sf_legendre_array> instead.'
      );

pp_def('gsl_sf_legendre_sphPlm_array',
       GenericTypes => [D],
       OtherPars =>'int n=>num; int m',
       Pars=>'double x(); double [o]y(num)',
       Code =>'
GSLERR(gsl_sf_legendre_sphPlm_array,($COMP(n)-2+$COMP(m),$COMP(m),$x(),$P(y)))
',
       Doc =>'P_lm(x), normalized properly for use in spherical harmonics for l from 0 to n-2+m.
gsl_sf_legendre_sphPlm_array has been deprecated in GSL version 2.0. It is included here for backwards compatability and may be removed in a future release.  New code should use L<gsl_sf_legendre_array> instead.'
      );

} else {
die("Could not determine GSL version from gsl-config, so can not determine which legendre array functions to define.");
}

pp_def('gsl_sf_legendre_sphPlm',
       GenericTypes => [D],
       OtherPars =>'int l; int m',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_legendre_sphPlm_e,($COMP(l),$COMP(m),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'P_lm(x), normalized properly for use in spherical harmonics'
      );

pp_def('gsl_sf_conicalP_half',
       GenericTypes => [D],
       OtherPars =>'double lambda',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_conicalP_half_e,($COMP(lambda),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Irregular Spherical Conical Function P^{1/2}_{-1/2 + I lambda}(x)'
      );

pp_def('gsl_sf_conicalP_mhalf',
       GenericTypes => [D],
       OtherPars =>'double lambda',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_conicalP_mhalf_e,($COMP(lambda),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Regular Spherical Conical Function P^{-1/2}_{-1/2 + I lambda}(x)'
      );

pp_def('gsl_sf_conicalP_0',
       GenericTypes => [D],
       OtherPars =>'double lambda',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_conicalP_0_e,($COMP(lambda),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Conical Function P^{0}_{-1/2 + I lambda}(x)'
      );

pp_def('gsl_sf_conicalP_1',
       GenericTypes => [D],
       OtherPars =>'double lambda',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_conicalP_1_e,($COMP(lambda),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Conical Function P^{1}_{-1/2 + I lambda}(x)'
      );

pp_def('gsl_sf_conicalP_sph_reg',
       GenericTypes => [D],
       OtherPars =>'int l; double lambda',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_conicalP_sph_reg_e,($COMP(l),$COMP(lambda),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Regular Spherical Conical Function P^{-1/2-l}_{-1/2 + I lambda}(x)'
      );

pp_def('gsl_sf_conicalP_cyl_reg_e',
       GenericTypes => [D],
       OtherPars =>'int m; double lambda',
       Pars=>'double x(); double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_conicalP_cyl_reg_e,($COMP(m),$COMP(lambda),$x(),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'Regular Cylindrical Conical Function P^{-m}_{-1/2 + I lambda}(x)'
      );

pp_def('gsl_sf_legendre_H3d',
       GenericTypes => [D],
       OtherPars =>'int l; double lambda; double eta',
       Pars=>'double [o]y(); double [o]e()',
       Code =>'
gsl_sf_result r;
GSLERR(gsl_sf_legendre_H3d_e,($COMP(l),$COMP(lambda),$COMP(eta),&r))
$y() = r.val;
$e() = r.err; 
',
       Doc =>'lth radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space.'
      );

pp_def('gsl_sf_legendre_H3d_array',
       GenericTypes => [D],
       OtherPars =>'int l=>num; double lambda; double eta',
       Pars=>'double [o]y(num)',
       Code =>'
GSLERR(gsl_sf_legendre_H3d_array,($COMP(l)-1,$COMP(lambda),$COMP(eta),$P(y)))
',
       Doc =>'Array of H3d(ell), for l from 0 to n-1.'
      );

pp_addpm({At=>Bot},<<'EOD');

=head1 AUTHOR

This file copyright (C) 1999 Christian Pellegrin <chri@infis.univ.trieste.it>
All rights reserved. There
is no warranty. You are allowed to redistribute this software /
documentation under certain conditions. For details, see the file
COPYING in the PDL distribution. If this file is separated from the
PDL distribution, the copyright notice should be included in the file.

The GSL SF modules were written by G. Jungman.

=cut


EOD

pp_add_boot('gsl_set_error_handler_off();
');


pp_done();