Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[nu,mu1,mu2,darg,mode,dpn,ipn,isig,ierror]=dxnrmp(nu,mu1,mu2,darg,mode,dpn,ipn,isig,ierror);
```function [nu,mu1,mu2,darg,mode,dpn,ipn,isig,ierror]=dxnrmp(nu,mu1,mu2,darg,mode,dpn,ipn,isig,ierror);
persistent c1 c2 dk gt i ip ip1 ip2 j k mu p p1 p2 p3 s sx t tx x ;

if isempty(i), i=0; end;
if isempty(ip), ip=0; end;
if isempty(ip1), ip1=0; end;
if isempty(ip2), ip2=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(mu), mu=0; end;
%***BEGIN PROLOGUE  DXNRMP
%***PURPOSE  Compute normalized Legendre polynomials.
%***LIBRARY   SLATEC
%***CATEGORY  C3A2, C9
%***TYPE      doubleprecision (XNRMP-S, DXNRMP-D)
%***KEYWORDS  LEGENDRE FUNCTIONS
%***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
%           Smith, John M., (NBS and George Mason University)
%***DESCRIPTION
%
%        subroutine TO CALCULATE NORMALIZED LEGENDRE POLYNOMIALS
%        (XNRMP is single-precision version)
%        DXNRMP calculates normalized Legendre polynomials of varying
%        order and fixed argument and degree. The order MU and degree
%        NU are non-negative integers and the argument is real. Because
%        the algorithm requires the use of numbers outside the normal
%        machine range, this subroutine employs a special arithmetic
%        called extended-range arithmetic. See J.M. Smith, F.W.J. Olver,
%        and D.W. Lozier, Extended-Range Arithmetic and Normalized
%        Legendre Polynomials, ACM Transactions on Mathematical Soft-
%        ware, 93-105, March 1981, for a complete description of the
%        algorithm and special arithmetic. Also see program comments
%        in DXSET.
%
%        The normalized Legendre polynomials are multiples of the
%        associated Legendre polynomials of the first kind where the
%        normalizing coefficients are chosen so as to make the integral
%        from -1 to 1 of the square of each function equal to 1. See
%        E. Jahnke, F. Emde and F. Losch, Tables of Higher Functions,
%        McGraw-Hill, New York, 1960, p. 121.
%
%        The input values to DXNRMP are NU, MU1, MU2, DARG, and MODE.
%        These must satisfy
%          1. NU .GE. 0 specifies the degree of the normalized Legendre
%             polynomial that is wanted.
%          2. MU1 .GE. 0 specifies the lowest-order normalized Legendre
%             polynomial that is wanted.
%          3. MU2 .GE. MU1 specifies the highest-order normalized Leg-
%             endre polynomial that is wanted.
%         4a. MODE = 1 and -1.0D0 .LE. DARG .LE. 1.0D0 specifies that
%             Normalized Legendre(NU, MU, DARG) is wanted for MU = MU1,
%             MU1 + 1, ..., MU2.
%         4b. MODE = 2 and -3.14159... .LT. DARG .LT. 3.14159... spec-
%             ifies that Normalized Legendre(NU, MU, COS(DARG)) is
%             wanted for MU = MU1, MU1 + 1, ..., MU2.
%
%        The output of DXNRMP consists of the two vectors DPN and IPN
%        and the error estimate ISIG. The computed values are stored as
%        extended-range numbers such that
%             (DPN(1),IPN(1))=NORMALIZED LEGENDRE(NU,MU1,DX)
%             (DPN(2),IPN(2))=NORMALIZED LEGENDRE(NU,MU1+1,DX)
%                .
%                .
%             (DPN(K),IPN(K))=NORMALIZED LEGENDRE(NU,MU2,DX)
%        where K = MU2 - MU1 + 1 and DX = DARG or COS(DARG) according
%        to whether MODE = 1 or 2. Finally, ISIG is an estimate of the
%        number of decimal digits lost through rounding errors in the
%        computation. For example if DARG is accurate to 12 significant
%        decimals, then the computed function values are accurate to
%        12 - ISIG significant decimals (except in neighborhoods of
%        zeros).
%
%        The interpretation of (DPN(I),IPN(I)) is DPN(I)*(IR**IPN(I))
%        where IR is the internal radix of the computer arithmetic. When
%        IPN(I) = 0 the value of the normalized Legendre polynomial is
%        contained entirely in DPN(I) and subsequent double-precision
%        computations can be performed without further consideration of
%        extended-range arithmetic. However, if IPN(I) .NE. 0 the corre-
%        sponding value of the normalized Legendre polynomial cannot be
%        represented in double-precision because of overflow or under-
%        flow. THE USER MUST TEST IPN(I) IN HIS/HER PROGRAM. In the case
%        that IPN(I) is nonzero, the user could rewrite his/her program
%        to use extended range arithmetic.
%
%
%
%        The interpretation of (DPN(I),IPN(I)) can be changed to
%        DPN(I)*(10**IPN(I)) by calling the extended-range subroutine
%        DXCON. This should be done before printing the computed values.
%        As an example of usage, the Fortran coding
%              J = K
%              DO 20 I = 1, K
%              CALL DXCON(DPN(I), IPN(I),IERROR)
%              IF (IERROR.NE.0) RETURN
%              PRINT 10, DPN(I), IPN(I)
%           10 FORMAT(1X, D30.18 , I15)
%              IF ((IPN(I) .EQ. 0) .OR. (J .LT. K)) GO TO 20
%              J = I - 1
%           20 CONTINUE
%        will print all computed values and determine the largest J
%        such that IPN(1) = IPN(2) = ... = IPN(J) = 0. Because of the
%        change of representation caused by calling DXCON, (DPN(I),
%        IPN(I)) for I = J+1, J+2, ... cannot be used in subsequent
%        extended-range computations.
%
%        IERROR is an error indicator. If no errors are detected,
%        IERROR=0 when control returns to the calling routine. If
%        an error is detected, IERROR is returned as nonzero. The
%        calling routine must check the value of IERROR.
%
%        If IERROR=212 or 213, invalid input was provided to DXNRMP.
%        If IERROR=201,202,203, or 204, invalid input was provided
%        to DXSET.
%        If IERROR=205 or 206, an internal consistency error occurred
%        in DXSET (probably due to a software malfunction in the
%        library routine I1MACH).
%        If IERROR=207, an overflow or underflow of an extended-range
%        number was detected in DXADJ.
%        If IERROR=208, an overflow or underflow of an extended-range
%        number was detected in DXC210.
%
%***REFERENCES  Smith, Olver and Lozier, Extended-Range Arithmetic and
%                 Normalized Legendre Polynomials, ACM Trans on Math
%                 Softw, v 7, n 1, March 1981, pp 93--105.
%***REVISION HISTORY  (YYMMDD)
%   820712  DATE WRITTEN
%   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
%   901019  Revisions to prologue.  (DWL and WRB)
%   901106  Changed all specific intrinsics to generic.  (WRB)
%           Corrected order of sections in prologue and added TYPE
%           section.  (WRB)
%           CALLs to XERROR changed to CALLs to XERMSG.  (WRB)
%   920127  Revised PURPOSE section of prologue.  (DWL)
%***end PROLOGUE  DXNRMP
if isempty(gt), gt=0; end;
dpn_shape=size(dpn);dpn=reshape(dpn,1,[]);
ipn_shape=size(ipn);ipn=reshape(ipn,1,[]);
if isempty(c1), c1=0; end;
if isempty(c2), c2=0; end;
if isempty(p), p=0; end;
if isempty(p1), p1=0; end;
if isempty(p2), p2=0; end;
if isempty(p3), p3=0; end;
if isempty(s), s=0; end;
if isempty(sx), sx=0; end;
if isempty(t), t=0; end;
if isempty(tx), tx=0; end;
if isempty(x), x=0; end;
if isempty(dk), dk=0; end;
% CALL DXSET TO INITIALIZE EXTENDED-RANGE ARITHMETIC (SEE DXSET
% LISTING FOR DETAILS)
%***FIRST EXECUTABLE STATEMENT  DXNRMP
ierror = 0;
[dumvar1,dumvar2,dumvar3,dumvar4,ierror]=dxset(0,0,0.0d0,0,ierror);
if( ierror==0 )
%
%        TEST FOR PROPER INPUT VALUES.
%
if( nu>=0 )
if( mu1>=0 )
if( mu1<=mu2 )
while (1);
gt=0;
if( nu~=0 )
if( mode<1 || mode>2 )
xermsg('SLATEC','DXNRMP','NU, MU1, MU2 or MODE not valid',212,1);
ierror = 212;
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
if( mode==2 )
if( abs(darg)<=4.0d0.*atan(1.0d0) )
if( darg==0.0d0 )
break;
end;
x = cos(darg);
sx = abs(sin(darg));
tx = x./sx;
isig = fix(log10(2.0d0.*nu.*(5.0d0+abs(darg.*tx))));
gt=1;
end;
elseif( abs(darg)<=1.0d0 ) ;
if( abs(darg)==1.0d0 )
break;
end;
x = darg;
sx = sqrt((1.0d0+abs(x)).*((0.5d0-abs(x))+0.5d0));
tx = x./sx;
isig = fix(log10(2.0d0.*nu.*(5.0d0+tx.^2)));
gt=1;
end;
if(gt==0)
xermsg('SLATEC','DXNRMP','DARG out of range',213,1);
ierror = 213;
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
%
%        BEGIN CALCULATION
%
mu = fix(mu2);
i = fix(mu2 - mu1 + 1);
%
%        IF MU.GT.NU, NORMALIZED LEGENDRE(NU,MU,X)=0.
%
while( mu>nu );
dpn(i) = 0.0d0;
ipn(i) = 0;
i = fix(i - 1);
mu = fix(mu - 1);
if( i<=0 )
isig = 0;
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
end;
mu = fix(nu);
%
%        P1 = 0. = NORMALIZED LEGENDRE(NU,NU+1,X)
%
p1 = 0.0d0;
ip1 = 0;
%
%        CALCULATE P2 = NORMALIZED LEGENDRE(NU,NU,X)
%
p2 = 1.0d0;
ip2 = 0;
p3 = 0.5d0;
dk = 2.0d0;
for j = 1 : nu;
p3 =((dk+1.0d0)./dk).*p3;
p2 = p2.*sx;
if( ierror~=0 )
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
dk = dk + 2.0d0;
end; j = fix(nu+1);
p2 = p2.*sqrt(p3);
if( ierror~=0 )
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
s = 2.0d0.*tx;
t = 1.0d0./nu;
gt=0;
if( mu2>=nu )
dpn(i) = p2;
ipn(i) = fix(ip2);
i = fix(i - 1);
if( i==0 )
gt=1;
end;
end;
%
%        RECURRENCE PROCESS
%
if(gt==0)
while( true );
p = mu.*t;
c1 = 1.0d0./sqrt((1.0d0-p+t).*(1.0d0+p));
c2 = s.*p.*c1.*p2;
c1 = -sqrt((1.0d0+p+t).*(1.0d0-p)).*c1.*p1;
if( ierror~=0 )
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
mu = fix(mu - 1);
if( mu<=mu2 )
%
%
dpn(i) = p;
ipn(i) = fix(ip);
i = fix(i - 1);
if( i==0 )
break;
end;
end;
p1 = p2;
ip1 = fix(ip2);
p2 = p;
ip2 = fix(ip);
if( mu<=mu1 )
break;
end;
end;
end;
%
%
k = fix(mu2 - mu1 + 1);
for i = 1 : k;
[dpn(i),ipn(i),ierror]=dxred(dpn(i),ipn(i),ierror);
if( ierror~=0 )
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
end; i = fix(k+1);
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
break;
end;
%
%        SPECIAL CASE WHEN X=-1 OR +1, OR NU=0.
%
k = fix(mu2 - mu1 + 1);
for i = 1 : k;
dpn(i) = 0.0d0;
ipn(i) = 0;
end; i = fix(k+1);
isig = 0;
if( mu1<=0 )
isig = 1;
dpn(1) = sqrt(nu+0.5d0);
ipn(1) = 0;
if( rem(nu,2)~=0 )
if( mode~=1 || darg~=1.0d0 )
if( mode~=2 )
dpn(1) = -dpn(1);
end;
end;
end;
end;
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
end;
end;
%
%          ERROR PRINTOUTS AND TERMINATION.
%
xermsg('SLATEC','DXNRMP','NU, MU1, MU2 or MODE not valid',212,1);
ierror = 212;
end;
dpn_shape=zeros(dpn_shape);dpn_shape(:)=dpn(1:numel(dpn_shape));dpn=dpn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
end %subroutine dxnrmp
%DECK DXPMU
```