| [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.
%
%***SEE ALSO XSET
%***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.
%***ROUTINES CALLED XADD, XADJ, XERMSG, XRED, XSET
%***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;
[p2,ip2,ierror]=xadj(p2,ip2,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;
rk = rk + 2.0;
end; j = fix(nu+1);
p2 = p2.*sqrt(p3);
[p2,ip2,ierror]=xadj(p2,ip2,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;
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;
[c2,ip2,c1,ip1,p,ip,ierror]=xadd(c2,ip2,c1,ip1,p,ip,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;
mu = fix(mu - 1);
if( mu<=mu2 )
%
% STORE IN ARRAY SPN FOR RETURN TO CALLING ROUTINE.
%
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;
%
% RETURN TO CALLING PROGRAM
%
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
|
|
