2
SUBROUTINE DXPQNU (NU1, NU2, MU, THETA, ID, PQA, IPQA, IERROR)
3
C***BEGIN PROLOGUE DXPQNU
5
C***PURPOSE To compute the values of Legendre functions for DXLEGF.
6
C This subroutine calculates initial values of P or Q using
7
C power series, then performs forward nu-wise recurrence to
8
C obtain P(-MU,NU,X), Q(0,NU,X), or Q(1,NU,X). The nu-wise
9
C recurrence is stable for P for all mu and for Q for mu=0,1.
12
C***TYPE DOUBLE PRECISION (XPQNU-S, DXPQNU-D)
13
C***KEYWORDS LEGENDRE FUNCTIONS
14
C***AUTHOR Smith, John M., (NBS and George Mason University)
15
C***ROUTINES CALLED DXADD, DXADJ, DXPSI
16
C***COMMON BLOCKS DXBLK1
17
C***REVISION HISTORY (YYMMDD)
19
C 890126 Revised to meet SLATEC CML recommendations. (DWL and JMS)
20
C 901019 Revisions to prologue. (DWL and WRB)
21
C 901106 Changed all specific intrinsics to generic. (WRB)
22
C Corrected order of sections in prologue and added TYPE
24
C 920127 Revised PURPOSE section of prologue. (DWL)
25
C***END PROLOGUE DXPQNU
26
DOUBLE PRECISION A,NU,NU1,NU2,PQ,PQA,DXPSI,R,THETA,W,X,X1,X2,XS,
28
DOUBLE PRECISION DI,DMU,PQ1,PQ2,FACTMU,FLOK
29
DIMENSION PQA(*),IPQA(*)
30
COMMON /DXBLK1/ NBITSF
33
C J0, IPSIK, AND IPSIX ARE INITIALIZED IN THIS SUBROUTINE.
34
C J0 IS THE NUMBER OF TERMS USED IN SERIES EXPANSION
35
C IN SUBROUTINE DXPQNU.
36
C IPSIK, IPSIX ARE VALUES OF K AND X RESPECTIVELY
37
C USED IN THE CALCULATION OF THE DXPSI FUNCTION.
39
C***FIRST EXECUTABLE STATEMENT DXPQNU
45
C FIND NU IN INTERVAL [-.5,.5) IF ID=2 ( CALCULATION OF Q )
47
IF(NU.GE..5D0) NU=NU-1.D0
48
C FIND NU IN INTERVAL (-1.5,-.5] IF ID=1,3, OR 4 ( CALC. OF P )
49
IF(ID.NE.2.AND.NU.GT.-.5D0) NU=NU-1.D0
50
C CALCULATE MU FACTORIAL
58
50 CALL DXADJ(FACTMU,IF,IERROR)
59
IF (IERROR.NE.0) RETURN
60
60 IF(K.EQ.0) FACTMU=1.D0
64
C Y=SIN(THETA/2)**2=(1-X)/2=.5-.5*X
65
C R=TAN(THETA/2)=SQRT((1-X)/(1+X)
71
C USE ASCENDING SERIES TO CALCULATE TWO VALUES OF P OR Q
72
C FOR USE AS STARTING VALUES IN RECURRENCE RELATION.
79
C SERIES FOR P ( ID = 1, 3, OR 4 )
80
C P(-MU,NU,X)=1./FACTORIAL(MU)*SQRT(((1.-X)/(1.+X))**MU)
81
C *SUM(FROM 0 TO J0-1)A(J)*(.5-.5*X)**J
89
A=A*Y*(DI-2.D0-NU)*(DI-1.D0+NU)/((DI-1.D0+DMU)*(DI-1.D0))
90
CALL DXADJ(A,IA,IERROR)
91
IF (IERROR.NE.0) RETURN
92
IF(A.EQ.0.D0) GO TO 66
93
CALL DXADD(PQ,IPQ,A,IA,PQ,IPQ,IERROR)
94
IF (IERROR.NE.0) RETURN
103
77 CALL DXADJ(X1,IPQ,IERROR)
104
IF (IERROR.NE.0) RETURN
107
CALL DXADJ(PQ,IPQ,IERROR)
108
IF (IERROR.NE.0) RETURN
111
C Z=-LN(R)=.5*LN((1+X)/(1-X))
114
W=DXPSI(NU+1.D0,IPSIK,IPSIX)
117
C SERIES SUMMATION FOR Q ( ID = 2 )
118
C Q(0,NU,X)=SUM(FROM 0 TO J0-1)((.5*LN((1+X)/(1-X))
119
C +DXPSI(J+1,IPSIK,IPSIX)-DXPSI(NU+1,IPSIK,IPSIX)))*A(J)*(.5-.5*X)**J
121
C Q(1,NU,X)=-SQRT(1./(1.-X**2))+SQRT((1-X)/(1+X))
122
C *SUM(FROM 0 T0 J0-1)(-NU*(NU+1)/2*LN((1+X)/(1-X))
123
C +(J-NU)*(J+NU+1)/(2*(J+1))+NU*(NU+1)*
124
C (DXPSI(NU+1,IPSIK,IPSIX)-DXPSI(J+1,IPSIK,IPSIX))*A(J)*(.5-.5*X)**J
126
C NOTE, IN THIS LOOP K=J+1
135
A=A*Y*(FLOK-2.D0-NU)*(FLOK-1.D0+NU)/((FLOK-1.D0+DMU)*(FLOK-1.D0))
136
CALL DXADJ(A,IA,IERROR)
137
IF (IERROR.NE.0) RETURN
140
X1=(DXPSI(FLOK,IPSIK,IPSIX)-W+Z)*A
142
CALL DXADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
143
IF (IERROR.NE.0) RETURN
145
83 X1=(NU*(NU+1.D0)*(Z-W+DXPSI(FLOK,IPSIK,IPSIX))+(NU-FLOK+1.D0)
146
1 *(NU+FLOK)/(2.D0*FLOK))*A
148
CALL DXADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
149
IF (IERROR.NE.0) RETURN
153
IF(MU.GE.1) CALL DXADD(PQ,IPQ,-XS,IXS,PQ,IPQ,IERROR)
154
IF (IERROR.NE.0) RETURN
162
IF(NU-1.5D0.LT.NU1) GO TO 120
166
IF(NU.GT.NU2+.5D0) RETURN
169
IF(NU.LT.NU1+.5D0) GO TO 130
173
IF(NU.GT.NU2+.5D0) RETURN
175
C FORWARD NU-WISE RECURRENCE FOR F(MU,NU,X) FOR FIXED MU
177
C (NU+MU+1)*F(MU,NU,X)=(2.*NU+1)*F(MU,NU,X)-(NU-MU)*F(MU,NU-1,X)
178
C WHERE F(MU,NU,X) MAY BE P(-MU,NU,X) OR IF MU IS REPLACED
179
C BY -MU THEN F(MU,NU,X) MAY BE Q(MU,NU,X).
180
C NOTE, IN THIS LOOP, NU=NU+1
182
130 X1=(2.D0*NU-1.D0)/(NU+DMU)*X*PQ1
183
X2=(NU-1.D0-DMU)/(NU+DMU)*PQ2
184
CALL DXADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
185
IF (IERROR.NE.0) RETURN
186
CALL DXADJ(PQ,IPQ,IERROR)
187
IF (IERROR.NE.0) RETURN