2
SUBROUTINE DXNRMP (NU, MU1, MU2, DARG, MODE, DPN, IPN, ISIG,
4
C***BEGIN PROLOGUE DXNRMP
5
C***PURPOSE Compute normalized Legendre polynomials.
8
C***TYPE DOUBLE PRECISION (XNRMP-S, DXNRMP-D)
9
C***KEYWORDS LEGENDRE FUNCTIONS
10
C***AUTHOR Lozier, Daniel W., (National Bureau of Standards)
11
C Smith, John M., (NBS and George Mason University)
14
C SUBROUTINE TO CALCULATE NORMALIZED LEGENDRE POLYNOMIALS
15
C (XNRMP is single-precision version)
16
C DXNRMP calculates normalized Legendre polynomials of varying
17
C order and fixed argument and degree. The order MU and degree
18
C NU are non-negative integers and the argument is real. Because
19
C the algorithm requires the use of numbers outside the normal
20
C machine range, this subroutine employs a special arithmetic
21
C called extended-range arithmetic. See J.M. Smith, F.W.J. Olver,
22
C and D.W. Lozier, Extended-Range Arithmetic and Normalized
23
C Legendre Polynomials, ACM Transactions on Mathematical Soft-
24
C ware, 93-105, March 1981, for a complete description of the
25
C algorithm and special arithmetic. Also see program comments
28
C The normalized Legendre polynomials are multiples of the
29
C associated Legendre polynomials of the first kind where the
30
C normalizing coefficients are chosen so as to make the integral
31
C from -1 to 1 of the square of each function equal to 1. See
32
C E. Jahnke, F. Emde and F. Losch, Tables of Higher Functions,
33
C McGraw-Hill, New York, 1960, p. 121.
35
C The input values to DXNRMP are NU, MU1, MU2, DARG, and MODE.
37
C 1. NU .GE. 0 specifies the degree of the normalized Legendre
38
C polynomial that is wanted.
39
C 2. MU1 .GE. 0 specifies the lowest-order normalized Legendre
40
C polynomial that is wanted.
41
C 3. MU2 .GE. MU1 specifies the highest-order normalized Leg-
42
C endre polynomial that is wanted.
43
C 4a. MODE = 1 and -1.0D0 .LE. DARG .LE. 1.0D0 specifies that
44
C Normalized Legendre(NU, MU, DARG) is wanted for MU = MU1,
46
C 4b. MODE = 2 and -3.14159... .LT. DARG .LT. 3.14159... spec-
47
C ifies that Normalized Legendre(NU, MU, COS(DARG)) is
48
C wanted for MU = MU1, MU1 + 1, ..., MU2.
50
C The output of DXNRMP consists of the two vectors DPN and IPN
51
C and the error estimate ISIG. The computed values are stored as
52
C extended-range numbers such that
53
C (DPN(1),IPN(1))=NORMALIZED LEGENDRE(NU,MU1,DX)
54
C (DPN(2),IPN(2))=NORMALIZED LEGENDRE(NU,MU1+1,DX)
57
C (DPN(K),IPN(K))=NORMALIZED LEGENDRE(NU,MU2,DX)
58
C where K = MU2 - MU1 + 1 and DX = DARG or COS(DARG) according
59
C to whether MODE = 1 or 2. Finally, ISIG is an estimate of the
60
C number of decimal digits lost through rounding errors in the
61
C computation. For example if DARG is accurate to 12 significant
62
C decimals, then the computed function values are accurate to
63
C 12 - ISIG significant decimals (except in neighborhoods of
66
C The interpretation of (DPN(I),IPN(I)) is DPN(I)*(IR**IPN(I))
67
C where IR is the internal radix of the computer arithmetic. When
68
C IPN(I) = 0 the value of the normalized Legendre polynomial is
69
C contained entirely in DPN(I) and subsequent double-precision
70
C computations can be performed without further consideration of
71
C extended-range arithmetic. However, if IPN(I) .NE. 0 the corre-
72
C sponding value of the normalized Legendre polynomial cannot be
73
C represented in double-precision because of overflow or under-
74
C flow. THE USER MUST TEST IPN(I) IN HIS/HER PROGRAM. In the case
75
C that IPN(I) is nonzero, the user could rewrite his/her program
76
C to use extended range arithmetic.
80
C The interpretation of (DPN(I),IPN(I)) can be changed to
81
C DPN(I)*(10**IPN(I)) by calling the extended-range subroutine
82
C DXCON. This should be done before printing the computed values.
83
C As an example of usage, the Fortran coding
86
C CALL DXCON(DPN(I), IPN(I),IERROR)
87
C IF (IERROR.NE.0) RETURN
88
C PRINT 10, DPN(I), IPN(I)
89
C 10 FORMAT(1X, D30.18 , I15)
90
C IF ((IPN(I) .EQ. 0) .OR. (J .LT. K)) GO TO 20
93
C will print all computed values and determine the largest J
94
C such that IPN(1) = IPN(2) = ... = IPN(J) = 0. Because of the
95
C change of representation caused by calling DXCON, (DPN(I),
96
C IPN(I)) for I = J+1, J+2, ... cannot be used in subsequent
97
C extended-range computations.
99
C IERROR is an error indicator. If no errors are detected,
100
C IERROR=0 when control returns to the calling routine. If
101
C an error is detected, IERROR is returned as nonzero. The
102
C calling routine must check the value of IERROR.
104
C If IERROR=212 or 213, invalid input was provided to DXNRMP.
105
C If IERROR=201,202,203, or 204, invalid input was provided
107
C If IERROR=205 or 206, an internal consistency error occurred
108
C in DXSET (probably due to a software malfunction in the
109
C library routine I1MACH).
110
C If IERROR=207, an overflow or underflow of an extended-range
111
C number was detected in DXADJ.
112
C If IERROR=208, an overflow or underflow of an extended-range
113
C number was detected in DXC210.
116
C***REFERENCES Smith, Olver and Lozier, Extended-Range Arithmetic and
117
C Normalized Legendre Polynomials, ACM Trans on Math
118
C Softw, v 7, n 1, March 1981, pp 93--105.
119
C***ROUTINES CALLED DXADD, DXADJ, DXRED, DXSET, XERMSG
120
C***REVISION HISTORY (YYMMDD)
121
C 820712 DATE WRITTEN
122
C 890126 Revised to meet SLATEC CML recommendations. (DWL and JMS)
123
C 901019 Revisions to prologue. (DWL and WRB)
124
C 901106 Changed all specific intrinsics to generic. (WRB)
125
C Corrected order of sections in prologue and added TYPE
127
C CALLs to XERROR changed to CALLs to XERMSG. (WRB)
128
C 920127 Revised PURPOSE section of prologue. (DWL)
129
C***END PROLOGUE DXNRMP
130
INTEGER NU, MU1, MU2, MODE, IPN, ISIG
131
DOUBLE PRECISION DARG, DPN
132
DIMENSION DPN(*), IPN(*)
133
DOUBLE PRECISION C1,C2,P,P1,P2,P3,S,SX,T,TX,X,DK
134
C CALL DXSET TO INITIALIZE EXTENDED-RANGE ARITHMETIC (SEE DXSET
135
C LISTING FOR DETAILS)
136
C***FIRST EXECUTABLE STATEMENT DXNRMP
138
CALL DXSET (0, 0, 0.0D0, 0,IERROR)
139
IF (IERROR.NE.0) RETURN
141
C TEST FOR PROPER INPUT VALUES.
143
IF (NU.LT.0) GO TO 110
144
IF (MU1.LT.0) GO TO 110
145
IF (MU1.GT.MU2) GO TO 110
146
IF (NU.EQ.0) GO TO 90
147
IF (MODE.LT.1 .OR. MODE.GT.2) GO TO 110
149
10 IF (ABS(DARG).GT.1.0D0) GO TO 120
150
IF (ABS(DARG).EQ.1.0D0) GO TO 90
152
SX = SQRT((1.0D0+ABS(X))*((0.5D0-ABS(X))+0.5D0))
154
ISIG = LOG10(2.0D0*NU*(5.0D0+TX**2))
156
20 IF (ABS(DARG).GT.4.0D0*ATAN(1.0D0)) GO TO 120
157
IF (DARG.EQ.0.0D0) GO TO 90
161
ISIG = LOG10(2.0D0*NU*(5.0D0+ABS(DARG*TX)))
168
C IF MU.GT.NU, NORMALIZED LEGENDRE(NU,MU,X)=0.
170
40 IF (MU.LE.NU) GO TO 50
175
IF (I .GT. 0) GO TO 40
180
C P1 = 0. = NORMALIZED LEGENDRE(NU,NU+1,X)
185
C CALCULATE P2 = NORMALIZED LEGENDRE(NU,NU,X)
192
P3 = ((DK+1.0D0)/DK)*P3
194
CALL DXADJ(P2, IP2,IERROR)
195
IF (IERROR.NE.0) RETURN
199
CALL DXADJ(P2, IP2,IERROR)
200
IF (IERROR.NE.0) RETURN
203
IF (MU2.LT.NU) GO TO 70
207
IF (I .EQ. 0) GO TO 140
212
C1 = 1.0D0/SQRT((1.0D0-P+T)*(1.0D0+P))
214
C1 = -SQRT((1.0D0+P+T)*(1.0D0-P))*C1*P1
215
CALL DXADD(C2, IP2, C1, IP1, P, IP,IERROR)
216
IF (IERROR.NE.0) RETURN
218
IF (MU.GT.MU2) GO TO 80
220
C STORE IN ARRAY DPN FOR RETURN TO CALLING ROUTINE.
225
IF (I .EQ. 0) GO TO 140
230
IF (MU.LE.MU1) GO TO 140
233
C SPECIAL CASE WHEN X=-1 OR +1, OR NU=0.
241
IF (MU1.GT.0) GO TO 160
243
DPN(1) = SQRT(NU+0.5D0)
245
IF (MOD(NU,2).EQ.0) GO TO 160
246
IF (MODE.EQ.1 .AND. DARG.EQ.1.0D0) GO TO 160
247
IF (MODE.EQ.2) GO TO 160
251
C ERROR PRINTOUTS AND TERMINATION.
253
110 CALL XERMSG ('SLATEC', 'DXNRMP', 'NU, MU1, MU2 or MODE not valid',
257
120 CALL XERMSG ('SLATEC', 'DXNRMP', 'DARG out of range', 213, 1)
261
C RETURN TO CALLING PROGRAM
263
140 K = MU2 - MU1 + 1
265
CALL DXRED(DPN(I),IPN(I),IERROR)
266
IF (IERROR.NE.0) RETURN
added
removed
deps/openlibm/slatec/dxnrmp.f
