Code covered by the BSD License

Highlights fromslatec

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

[nu,mu1,mu2,sarg,mode,spn,ipn,isig,ierror]=xnrmp(nu,mu1,mu2,sarg,mode,spn,ipn,isig,ierror);
```function [nu,mu1,mu2,sarg,mode,spn,ipn,isig,ierror]=xnrmp(nu,mu1,mu2,sarg,mode,spn,ipn,isig,ierror);
persistent c1 c2 i igo ip ip1 ip2 j k mu p p1 p2 p3 rk 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;
if isempty(igo), igo=0; end;
%***BEGIN PROLOGUE  XNRMP
%***PURPOSE  Compute normalized Legendre polynomials.
%***LIBRARY   SLATEC
%***CATEGORY  C3A2, C9
%***TYPE      SINGLE PRECISION (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
%        (DXNRMP is double-precision version)
%        XNRMP 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 XSET.
%
%        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 XNRMP are NU, MU1, MU2, SARG, 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.0 .LE. SARG .LE. 1.0 specifies that
%             Normalized Legendre(NU, MU, SARG) is wanted for MU = MU1,
%             MU1 + 1, ..., MU2.
%         4b. MODE = 2 and -3.14159... .LT. SARG .LT. 3.14159... spec-
%             ifies that Normalized Legendre(NU, MU, COS(SARG)) is want-
%             ed for MU = MU1, MU1 + 1, ..., MU2.
%
%        The output of XNRMP consists of the two vectors SPN and IPN
%        and the error estimate ISIG. The computed values are stored as
%        extended-range numbers such that
%             (SPN(1),IPN(1))=NORMALIZED LEGENDRE(NU,MU1,X)
%             (SPN(2),IPN(2))=NORMALIZED LEGENDRE(NU,MU1+1,X)
%                .
%                .
%             (SPN(K),IPN(K))=NORMALIZED LEGENDRE(NU,MU2,X)
%        where K = MU2 - MU1 + 1 and X = SARG or COS(SARG) 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 SARG 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 (SPN(I),IPN(I)) is SPN(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 SPN(I) and subsequent single-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 single-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 should try using double pre-
%        cision if it has a wider exponent range. If doubleprecision
%        fails, the user could rewrite his/her program to use extended-
%        range arithmetic.
%
%        The interpretation of (SPN(I),IPN(I)) can be changed to
%        SPN(I)*(10**IPN(I)) by calling the extended-range subroutine
%        XCON. 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 XCON(SPN(I), IPN(I),IERROR)
%              IF (IERROR.NE.0) RETURN
%              PRINT 10, SPN(I), IPN(I)
%           10 FORMAT(1X, E30.8 , 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 XCON, (SPN(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=112 or 113, invalid input was provided to XNRMP.
%        If IERROR=101,102,103, or 104, invalid input was provided
%        to XSET.
%        If IERROR=105 or 106, an internal consistency error occurred
%        in XSET (probably due to a software malfunction in the
%        library routine I1MACH).
%        If IERROR=107, an overflow or underflow of an extended-range
%        number was detected in XADJ.
%        If IERROR=108, an overflow or underflow of an extended-range
%        number was detected in XC210.
%
%***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
%   881020  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  XNRMP
spn_shape=size(spn);spn=reshape(spn,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(rk), rk=0; end;
% CALL XSET TO INITIALIZE EXTENDED-RANGE ARITHMETIC (SEE XSET
% LISTING FOR DETAILS)
%***FIRST EXECUTABLE STATEMENT  XNRMP
ierror = 0;
[dumvar1,dumvar2,dumvar3,dumvar4,ierror]=xset(0,0,0.0,0,ierror);
if( ierror==0 )
%
%        TEST FOR PROPER INPUT VALUES.
%
if( nu>=0 )
if( mu1>=0 )
if( mu1<=mu2 )
igo=0;
while (1);
if( nu~=0 )
while (1);
if( mode<1 || mode>2 )
xermsg('SLATEC','XNRMP','NU, MU1, MU2 or MODE not valid',112,1);
ierror = 112;
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
if( mode==2 )
if( abs(sarg)<=4.0.*atan(1.0) )
if( sarg==0.0 )
igo=1;
break;
end;
x = cos(sarg);
sx = abs(sin(sarg));
tx = x./sx;
isig = fix(log10(2.0.*nu.*(5.0+abs(sarg.*tx))));
break;
end;
elseif( abs(sarg)<=1.0 ) ;
if( abs(sarg)==1.0 )
igo=1;
break;
end;
x = sarg;
sx = sqrt((1.0+abs(x)).*((0.5-abs(x))+0.5));
tx = x./sx;
isig = fix(log10(2.0.*nu.*(5.0+tx.^2)));
break;
end;
xermsg('SLATEC','XNRMP','SARG out of range',113,1);
ierror = 113;
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
if(igo==1)
break;
end;
%
%        BEGIN CALCULATION
%
mu = fix(mu2);
i = fix(mu2 - mu1 + 1);
%
%        IF MU.GT.NU, NORMALIZED LEGENDRE(NU,MU,X)=0.
%
while( mu>nu );
spn(i) = 0.0;
ipn(i) = 0;
i = fix(i - 1);
mu = fix(mu - 1);
if( i<=0 )
isig = 0;
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_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.0;
ip1 = 0;
%
%        CALCULATE P2 = NORMALIZED LEGENDRE(NU,NU,X)
%
p2 = 1.0;
ip2 = 0;
p3 = 0.5;
rk = 2.0;
for j = 1 : nu;
p3 =((rk+1.0)./rk).*p3;
p2 = p2.*sx;
if( ierror~=0 )
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
rk = rk + 2.0;
end; j = fix(nu+1);
p2 = p2.*sqrt(p3);
if( ierror~=0 )
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
s = 2.0.*tx;
t = 1.0./nu;
if( mu2>=nu )
spn(i) = p2;
ipn(i) = fix(ip2);
i = fix(i - 1);
if( i==0 )
k = fix(mu2 - mu1 + 1);
for i = 1 : k;
[spn(i),ipn(i),ierror]=xred(spn(i),ipn(i),ierror);
if( ierror~=0 )
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
end; i = fix(k+1);
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
end;
%
%        RECURRENCE PROCESS
%
while( true );
p = mu.*t;
c1 = 1.0./sqrt((1.0-p+t).*(1.0+p));
c2 = s.*p.*c1.*p2;
c1 = -sqrt((1.0+p+t).*(1.0-p)).*c1.*p1;
if( ierror~=0 )
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_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 )
%
%
spn(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;
%
%
k = fix(mu2 - mu1 + 1);
for i = 1 : k;
[spn(i),ipn(i),ierror]=xred(spn(i),ipn(i),ierror);
if( ierror~=0 )
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
return;
end;
end; i = fix(k+1);
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_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;
spn(i) = 0.0;
ipn(i) = 0;
end; i = fix(k+1);
isig = 0;
if( mu1<=0 )
isig = 1;
spn(1) = sqrt(nu+0.5);
ipn(1) = 0;
if( rem(nu,2)~=0 )
if( mode~=1 || sarg~=1.0 )
if( mode~=2 )
spn(1) = -spn(1);
end;
end;
end;
end;
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_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','XNRMP','NU, MU1, MU2 or MODE not valid',112,1);
ierror = 112;
end;
spn_shape=zeros(spn_shape);spn_shape(:)=spn(1:numel(spn_shape));spn=spn_shape;
ipn_shape=zeros(ipn_shape);ipn_shape(:)=ipn(1:numel(ipn_shape));ipn=ipn_shape;
end %subroutine xnrmp
%DECK XPMU
```