*DECK XNRMP SUBROUTINE XNRMP (NU, MU1, MU2, SARG, MODE, SPN, IPN, ISIG, 1 IERROR) C***BEGIN PROLOGUE XNRMP C***PURPOSE Compute normalized Legendre polynomials. C***LIBRARY SLATEC C***CATEGORY C3A2, C9 C***TYPE SINGLE PRECISION (XNRMP-S, DXNRMP-D) C***KEYWORDS LEGENDRE FUNCTIONS C***AUTHOR Lozier, Daniel W., (National Bureau of Standards) C Smith, John M., (NBS and George Mason University) C***DESCRIPTION C C SUBROUTINE TO CALCULATE NORMALIZED LEGENDRE POLYNOMIALS C (DXNRMP is double-precision version) C XNRMP calculates normalized Legendre polynomials of varying C order and fixed argument and degree. The order MU and degree C NU are non-negative integers and the argument is real. Because C the algorithm requires the use of numbers outside the normal C machine range, this subroutine employs a special arithmetic C called extended-range arithmetic. See J.M. Smith, F.W.J. Olver, C and D.W. Lozier, Extended-Range Arithmetic and Normalized C Legendre Polynomials, ACM Transactions on Mathematical Soft- C ware, 93-105, March 1981, for a complete description of the C algorithm and special arithmetic. Also see program comments C in XSET. C C The normalized Legendre polynomials are multiples of the C associated Legendre polynomials of the first kind where the C normalizing coefficients are chosen so as to make the integral C from -1 to 1 of the square of each function equal to 1. See C E. Jahnke, F. Emde and F. Losch, Tables of Higher Functions, C McGraw-Hill, New York, 1960, p. 121. C C The input values to XNRMP are NU, MU1, MU2, SARG, and MODE. C These must satisfy C 1. NU .GE. 0 specifies the degree of the normalized Legendre C polynomial that is wanted. C 2. MU1 .GE. 0 specifies the lowest-order normalized Legendre C polynomial that is wanted. C 3. MU2 .GE. MU1 specifies the highest-order normalized Leg- C endre polynomial that is wanted. C 4a. MODE = 1 and -1.0 .LE. SARG .LE. 1.0 specifies that C Normalized Legendre(NU, MU, SARG) is wanted for MU = MU1, C MU1 + 1, ..., MU2. C 4b. MODE = 2 and -3.14159... .LT. SARG .LT. 3.14159... spec- C ifies that Normalized Legendre(NU, MU, COS(SARG)) is want- C ed for MU = MU1, MU1 + 1, ..., MU2. C C The output of XNRMP consists of the two vectors SPN and IPN C and the error estimate ISIG. The computed values are stored as C extended-range numbers such that C (SPN(1),IPN(1))=NORMALIZED LEGENDRE(NU,MU1,X) C (SPN(2),IPN(2))=NORMALIZED LEGENDRE(NU,MU1+1,X) C . C . C (SPN(K),IPN(K))=NORMALIZED LEGENDRE(NU,MU2,X) C where K = MU2 - MU1 + 1 and X = SARG or COS(SARG) according C to whether MODE = 1 or 2. Finally, ISIG is an estimate of the C number of decimal digits lost through rounding errors in the C computation. For example if SARG is accurate to 12 significant C decimals, then the computed function values are accurate to C 12 - ISIG significant decimals (except in neighborhoods of C zeros). C C The interpretation of (SPN(I),IPN(I)) is SPN(I)*(IR**IPN(I)) C where IR is the internal radix of the computer arithmetic. When C IPN(I) = 0 the value of the normalized Legendre polynomial is C contained entirely in SPN(I) and subsequent single-precision C computations can be performed without further consideration of C extended-range arithmetic. However, if IPN(I) .NE. 0 the corre- C sponding value of the normalized Legendre polynomial cannot be C represented in single-precision because of overflow or under- C flow. THE USER MUST TEST IPN(I) IN HIS/HER PROGRAM. In the case C that IPN(I) is nonzero, the user should try using double pre- C cision if it has a wider exponent range. If double precision C fails, the user could rewrite his/her program to use extended- C range arithmetic. C C The interpretation of (SPN(I),IPN(I)) can be changed to C SPN(I)*(10**IPN(I)) by calling the extended-range subroutine C XCON. This should be done before printing the computed values. C As an example of usage, the Fortran coding C J = K C DO 20 I = 1, K C CALL XCON(SPN(I), IPN(I),IERROR) C IF (IERROR.NE.0) RETURN C PRINT 10, SPN(I), IPN(I) C 10 FORMAT(1X, E30.8 , I15) C IF ((IPN(I) .EQ. 0) .OR. (J .LT. K)) GO TO 20 C J = I - 1 C 20 CONTINUE C will print all computed values and determine the largest J C such that IPN(1) = IPN(2) = ... = IPN(J) = 0. Because of the C change of representation caused by calling XCON, (SPN(I), C IPN(I)) for I = J+1, J+2, ... cannot be used in subsequent C extended-range computations. C C IERROR is an error indicator. If no errors are detected, C IERROR=0 when control returns to the calling routine. If C an error is detected, IERROR is returned as nonzero. The C calling routine must check the value of IERROR. C C If IERROR=112 or 113, invalid input was provided to XNRMP. C If IERROR=101,102,103, or 104, invalid input was provided C to XSET. C If IERROR=105 or 106, an internal consistency error occurred C in XSET (probably due to a software malfunction in the C library routine I1MACH). C If IERROR=107, an overflow or underflow of an extended-range C number was detected in XADJ. C If IERROR=108, an overflow or underflow of an extended-range C number was detected in XC210. C C***SEE ALSO XSET C***REFERENCES Smith, Olver and Lozier, Extended-Range Arithmetic and C Normalized Legendre Polynomials, ACM Trans on Math C Softw, v 7, n 1, March 1981, pp 93--105. C***ROUTINES CALLED XADD, XADJ, XERMSG, XRED, XSET C***REVISION HISTORY (YYMMDD) C 820712 DATE WRITTEN C 881020 Revised to meet SLATEC CML recommendations. (DWL and JMS) C 901019 Revisions to prologue. (DWL and WRB) C 901106 Changed all specific intrinsics to generic. (WRB) C Corrected order of sections in prologue and added TYPE C section. (WRB) C CALLs to XERROR changed to CALLs to XERMSG. (WRB) C 920127 Revised PURPOSE section of prologue. (DWL) C***END PROLOGUE XNRMP INTEGER NU, MU1, MU2, MODE, IPN, ISIG REAL SARG, SPN DIMENSION SPN(*), IPN(*) REAL C1,C2,P,P1,P2,P3,S,SX,T,TX,X,RK C CALL XSET TO INITIALIZE EXTENDED-RANGE ARITHMETIC (SEE XSET C LISTING FOR DETAILS) C***FIRST EXECUTABLE STATEMENT XNRMP IERROR=0 CALL XSET (0, 0, 0.0, 0,IERROR) IF (IERROR.NE.0) RETURN C C TEST FOR PROPER INPUT VALUES. C IF (NU.LT.0) GO TO 110 IF (MU1.LT.0) GO TO 110 IF (MU1.GT.MU2) GO TO 110 IF (NU.EQ.0) GO TO 90 IF (MODE.LT.1 .OR. MODE.GT.2) GO TO 110 GO TO (10, 20), MODE 10 IF (ABS(SARG).GT.1.0) GO TO 120 IF (ABS(SARG).EQ.1.0) GO TO 90 X = SARG SX = SQRT((1.0+ABS(X))*((0.5-ABS(X))+0.5)) TX = X/SX ISIG = LOG10(2.0*NU*(5.0+TX**2)) GO TO 30 20 IF (ABS(SARG).GT.4.0*ATAN(1.0)) GO TO 120 IF (SARG.EQ.0.0) GO TO 90 X = COS(SARG) SX = ABS(SIN(SARG)) TX = X/SX ISIG = LOG10(2.0*NU*(5.0+ABS(SARG*TX))) C C BEGIN CALCULATION C 30 MU = MU2 I = MU2 - MU1 + 1 C C IF MU.GT.NU, NORMALIZED LEGENDRE(NU,MU,X)=0. C 40 IF (MU.LE.NU) GO TO 50 SPN(I) = 0.0 IPN(I) = 0 I = I - 1 MU = MU - 1 IF (I .GT. 0) GO TO 40 ISIG = 0 GO TO 160 50 MU = NU C C P1 = 0. = NORMALIZED LEGENDRE(NU,NU+1,X) C P1 = 0.0 IP1 = 0 C C CALCULATE P2 = NORMALIZED LEGENDRE(NU,NU,X) C P2 = 1.0 IP2 = 0 P3 = 0.5 RK = 2.0 DO 60 J=1,NU P3 = ((RK+1.0)/RK)*P3 P2 = P2*SX CALL XADJ(P2, IP2,IERROR) IF (IERROR.NE.0) RETURN RK = RK + 2.0 60 CONTINUE P2 = P2*SQRT(P3) CALL XADJ(P2, IP2,IERROR) IF (IERROR.NE.0) RETURN S = 2.0*TX T = 1.0/NU IF (MU2.LT.NU) GO TO 70 SPN(I) = P2 IPN(I) = IP2 I = I - 1 IF (I .EQ. 0) GO TO 140 C C RECURRENCE PROCESS C 70 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 CALL XADD(C2, IP2, C1, IP1, P, IP,IERROR) IF (IERROR.NE.0) RETURN MU = MU - 1 IF (MU.GT.MU2) GO TO 80 C C STORE IN ARRAY SPN FOR RETURN TO CALLING ROUTINE. C SPN(I) = P IPN(I) = IP I = I - 1 IF (I .EQ. 0) GO TO 140 80 P1 = P2 IP1 = IP2 P2 = P IP2 = IP IF (MU.LE.MU1) GO TO 140 GO TO 70 C C SPECIAL CASE WHEN X=-1 OR +1, OR NU=0. C 90 K = MU2 - MU1 + 1 DO 100 I=1,K SPN(I) = 0.0 IPN(I) = 0 100 CONTINUE ISIG = 0 IF (MU1.GT.0) GO TO 160 ISIG = 1 SPN(1) = SQRT(NU+0.5) IPN(1) = 0 IF (MOD(NU,2).EQ.0) GO TO 160 IF (MODE.EQ.1 .AND. SARG.EQ.1.0) GO TO 160 IF (MODE.EQ.2) GO TO 160 SPN(1) = -SPN(1) GO TO 160 C C ERROR PRINTOUTS AND TERMINATION. C 110 CALL XERMSG ('SLATEC', 'XNRMP', 'NU, MU1, MU2 or MODE not valid', + 112, 1) IERROR=112 RETURN 120 CALL XERMSG ('SLATEC', 'XNRMP', 'SARG out of range', 113, 1) IERROR=113 RETURN C C RETURN TO CALLING PROGRAM C 140 K = MU2 - MU1 + 1 DO 150 I=1,K CALL XRED(SPN(I),IPN(I),IERROR) IF (IERROR.NE.0) RETURN 150 CONTINUE 160 RETURN END