*DECK DXLEGF SUBROUTINE DXLEGF (DNU1, NUDIFF, MU1, MU2, THETA, ID, PQA, IPQA, 1 IERROR) C***BEGIN PROLOGUE DXLEGF C***PURPOSE Compute normalized Legendre polynomials and associated C Legendre functions. C***LIBRARY SLATEC C***CATEGORY C3A2, C9 C***TYPE DOUBLE PRECISION (XLEGF-S, DXLEGF-D) C***KEYWORDS LEGENDRE FUNCTIONS C***AUTHOR Smith, John M., (NBS and George Mason University) C***DESCRIPTION C C DXLEGF: Extended-range Double-precision Legendre Functions C C A feature of the DXLEGF subroutine for Legendre functions is C the use of extended-range arithmetic, a software extension of C ordinary floating-point arithmetic that greatly increases the C exponent range of the representable numbers. This avoids the C need for scaling the solutions to lie within the exponent range C of the most restrictive manufacturer's hardware. The increased C exponent range is achieved by allocating an integer storage C location together with each floating-point storage location. C C The interpretation of the pair (X,I) where X is floating-point C and I is integer is X*(IR**I) where IR is the internal radix of C the computer arithmetic. C C This subroutine computes one of the following vectors: C C 1. Legendre function of the first kind of negative order, either C a. P(-MU1,NU,X), P(-MU1-1,NU,X), ..., P(-MU2,NU,X) or C b. P(-MU,NU1,X), P(-MU,NU1+1,X), ..., P(-MU,NU2,X) C 2. Legendre function of the second kind, either C a. Q(MU1,NU,X), Q(MU1+1,NU,X), ..., Q(MU2,NU,X) or C b. Q(MU,NU1,X), Q(MU,NU1+1,X), ..., Q(MU,NU2,X) C 3. Legendre function of the first kind of positive order, either C a. P(MU1,NU,X), P(MU1+1,NU,X), ..., P(MU2,NU,X) or C b. P(MU,NU1,X), P(MU,NU1+1,X), ..., P(MU,NU2,X) C 4. Normalized Legendre polynomials, either C a. PN(MU1,NU,X), PN(MU1+1,NU,X), ..., PN(MU2,NU,X) or C b. PN(MU,NU1,X), PN(MU,NU1+1,X), ..., PN(MU,NU2,X) C C where X = COS(THETA). C C The input values to DXLEGF are DNU1, NUDIFF, MU1, MU2, THETA, C and ID. These must satisfy C C DNU1 is DOUBLE PRECISION and greater than or equal to -0.5; C NUDIFF is INTEGER and non-negative; C MU1 is INTEGER and non-negative; C MU2 is INTEGER and greater than or equal to MU1; C THETA is DOUBLE PRECISION and in the half-open interval (0,PI/2]; C ID is INTEGER and equal to 1, 2, 3 or 4; C C and additionally either NUDIFF = 0 or MU2 = MU1. C C If ID=1 and NUDIFF=0, a vector of type 1a above is computed C with NU=DNU1. C C If ID=1 and MU1=MU2, a vector of type 1b above is computed C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1. C C If ID=2 and NUDIFF=0, a vector of type 2a above is computed C with NU=DNU1. C C If ID=2 and MU1=MU2, a vector of type 2b above is computed C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1. C C If ID=3 and NUDIFF=0, a vector of type 3a above is computed C with NU=DNU1. C C If ID=3 and MU1=MU2, a vector of type 3b above is computed C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1. C C If ID=4 and NUDIFF=0, a vector of type 4a above is computed C with NU=DNU1. C C If ID=4 and MU1=MU2, a vector of type 4b above is computed C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1. C C In each case the vector of computed Legendre function values C is returned in the extended-range vector (PQA(I),IPQA(I)). The C length of this vector is either MU2-MU1+1 or NUDIFF+1. C C Where possible, DXLEGF returns IPQA(I) as zero. In this case the C value of the Legendre function is contained entirely in PQA(I), C so it can be used in subsequent computations without further C consideration of extended-range arithmetic. If IPQA(I) is nonzero, C then the value of the Legendre function is not representable in C floating-point because of underflow or overflow. The program that C calls DXLEGF must test IPQA(I) to ensure correct usage. C C IERROR is an error indicator. If no errors are detected, IERROR=0 C when control returns to the calling routine. If an error is detected, C IERROR is returned as nonzero. The calling routine must check the C value of IERROR. C C If IERROR=210 or 211, invalid input was provided to DXLEGF. C If IERROR=201,202,203, or 204, invalid input was provided to DXSET. C If IERROR=205 or 206, an internal consistency error occurred in C DXSET (probably due to a software malfunction in the library routine C I1MACH). C If IERROR=207, an overflow or underflow of an extended-range number C was detected in DXADJ. C If IERROR=208, an overflow or underflow of an extended-range number C was detected in DXC210. C C***SEE ALSO DXSET C***REFERENCES Olver and Smith, Associated Legendre Functions on the C Cut, J Comp Phys, v 51, n 3, Sept 1983, pp 502--518. C 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 DXPMU, DXPMUP, DXPNRM, DXPQNU, DXQMU, DXQNU, DXRED, C DXSET, XERMSG C***REVISION HISTORY (YYMMDD) C 820728 DATE WRITTEN C 890126 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 DXLEGF DOUBLE PRECISION PQA,DNU1,DNU2,SX,THETA,X,PI2 DIMENSION PQA(*),IPQA(*) C C***FIRST EXECUTABLE STATEMENT DXLEGF IERROR=0 CALL DXSET (0, 0, 0.0D0, 0,IERROR) IF (IERROR.NE.0) RETURN PI2=2.D0*ATAN(1.D0) C C ZERO OUTPUT ARRAYS C L=(MU2-MU1)+NUDIFF+1 DO 290 I=1,L PQA(I)=0.D0 290 IPQA(I)=0 C C CHECK FOR VALID INPUT VALUES C IF(NUDIFF.LT.0) GO TO 400 IF(DNU1.LT.-.5D0) GO TO 400 IF(MU2.LT.MU1) GO TO 400 IF(MU1.LT.0) GO TO 400 IF(THETA.LE.0.D0.OR.THETA.GT.PI2) GO TO 420 IF(ID.LT.1.OR.ID.GT.4) GO TO 400 IF((MU1.NE.MU2).AND.(NUDIFF.GT.0)) GO TO 400 C C IF DNU1 IS NOT AN INTEGER, NORMALIZED P(MU,DNU,X) C CANNOT BE CALCULATED. IF DNU1 IS AN INTEGER AND C MU1.GT.DNU2 THEN ALL VALUES OF P(+MU,DNU,X) AND C NORMALIZED P(MU,NU,X) WILL BE ZERO. C DNU2=DNU1+NUDIFF IF((ID.EQ.3).AND.(MOD(DNU1,1.D0).NE.0.D0)) GO TO 295 IF((ID.EQ.4).AND.(MOD(DNU1,1.D0).NE.0.D0)) GO TO 400 IF((ID.EQ.3.OR.ID.EQ.4).AND.MU1.GT.DNU2) RETURN 295 CONTINUE C X=COS(THETA) SX=1.D0/SIN(THETA) IF(ID.EQ.2) GO TO 300 IF(MU2-MU1.LE.0) GO TO 360 C C FIXED NU, VARIABLE MU C CALL DXPMU TO CALCULATE P(-MU1,NU,X),....,P(-MU2,NU,X) C CALL DXPMU(DNU1,DNU2,MU1,MU2,THETA,X,SX,ID,PQA,IPQA,IERROR) IF (IERROR.NE.0) RETURN GO TO 380 C 300 IF(MU2.EQ.MU1) GO TO 320 C C FIXED NU, VARIABLE MU C CALL DXQMU TO CALCULATE Q(MU1,NU,X),....,Q(MU2,NU,X) C CALL DXQMU(DNU1,DNU2,MU1,MU2,THETA,X,SX,ID,PQA,IPQA,IERROR) IF (IERROR.NE.0) RETURN GO TO 390 C C FIXED MU, VARIABLE NU C CALL DXQNU TO CALCULATE Q(MU,DNU1,X),....,Q(MU,DNU2,X) C 320 CALL DXQNU(DNU1,DNU2,MU1,THETA,X,SX,ID,PQA,IPQA,IERROR) IF (IERROR.NE.0) RETURN GO TO 390 C C FIXED MU, VARIABLE NU C CALL DXPQNU TO CALCULATE P(-MU,DNU1,X),....,P(-MU,DNU2,X) C 360 CALL DXPQNU(DNU1,DNU2,MU1,THETA,ID,PQA,IPQA,IERROR) IF (IERROR.NE.0) RETURN C C IF ID = 3, TRANSFORM P(-MU,NU,X) VECTOR INTO C P(MU,NU,X) VECTOR. C 380 IF(ID.EQ.3) CALL DXPMUP(DNU1,DNU2,MU1,MU2,PQA,IPQA,IERROR) IF (IERROR.NE.0) RETURN C C IF ID = 4, TRANSFORM P(-MU,NU,X) VECTOR INTO C NORMALIZED P(MU,NU,X) VECTOR. C IF(ID.EQ.4) CALL DXPNRM(DNU1,DNU2,MU1,MU2,PQA,IPQA,IERROR) IF (IERROR.NE.0) RETURN C C PLACE RESULTS IN REDUCED FORM IF POSSIBLE C AND RETURN TO MAIN PROGRAM. C 390 DO 395 I=1,L CALL DXRED(PQA(I),IPQA(I),IERROR) IF (IERROR.NE.0) RETURN 395 CONTINUE RETURN C C ***** ERROR TERMINATION ***** C 400 CALL XERMSG ('SLATEC', 'DXLEGF', + 'DNU1, NUDIFF, MU1, MU2, or ID not valid', 210, 1) IERROR=210 RETURN 420 CALL XERMSG ('SLATEC', 'DXLEGF', 'THETA out of range', 211, 1) IERROR=211 RETURN END