12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394 |
- SUBROUTINE ZWRSK(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI,
- * TOL, ELIM, ALIM)
- C***BEGIN PROLOGUE ZWRSK
- C***REFER TO ZBESI,ZBESK
- C
- C ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
- C NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
- C
- C***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,ZABS
- C***END PROLOGUE ZWRSK
- C COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR
- DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI,
- * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT,
- * STI, STR, TOL, YI, YR, ZRI, ZRR, ZABS, D1MACH
- INTEGER I, KODE, N, NW, NZ
- DIMENSION YR(N), YI(N), CWR(2), CWI(2)
- C-----------------------------------------------------------------------
- C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
- C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
- C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
- C-----------------------------------------------------------------------
- NZ = 0
- CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM)
- IF (NW.NE.0) GO TO 50
- CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL)
- C-----------------------------------------------------------------------
- C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
- C R(FNU+J-1,Z)=Y(J), J=1,...,N
- C-----------------------------------------------------------------------
- CINUR = 1.0D0
- CINUI = 0.0D0
- IF (KODE.EQ.1) GO TO 10
- CINUR = DCOS(ZRI)
- CINUI = DSIN(ZRI)
- 10 CONTINUE
- C-----------------------------------------------------------------------
- C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
- C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
- C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
- C THE RESULT IS ON SCALE.
- C-----------------------------------------------------------------------
- ACW = ZABS(COMPLEX(CWR(2),CWI(2)))
- ASCLE = 1.0D+3*D1MACH(1)/TOL
- CSCLR = 1.0D0
- IF (ACW.GT.ASCLE) GO TO 20
- CSCLR = 1.0D0/TOL
- GO TO 30
- 20 CONTINUE
- ASCLE = 1.0D0/ASCLE
- IF (ACW.LT.ASCLE) GO TO 30
- CSCLR = TOL
- 30 CONTINUE
- C1R = CWR(1)*CSCLR
- C1I = CWI(1)*CSCLR
- C2R = CWR(2)*CSCLR
- C2I = CWI(2)*CSCLR
- STR = YR(1)
- STI = YI(1)
- C-----------------------------------------------------------------------
- C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS
- C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
- C-----------------------------------------------------------------------
- PTR = STR*C1R - STI*C1I
- PTI = STR*C1I + STI*C1R
- PTR = PTR + C2R
- PTI = PTI + C2I
- CTR = ZRR*PTR - ZRI*PTI
- CTI = ZRR*PTI + ZRI*PTR
- ACT = ZABS(COMPLEX(CTR,CTI))
- RACT = 1.0D0/ACT
- CTR = CTR*RACT
- CTI = -CTI*RACT
- PTR = CINUR*RACT
- PTI = CINUI*RACT
- CINUR = PTR*CTR - PTI*CTI
- CINUI = PTR*CTI + PTI*CTR
- YR(1) = CINUR*CSCLR
- YI(1) = CINUI*CSCLR
- IF (N.EQ.1) RETURN
- DO 40 I=2,N
- PTR = STR*CINUR - STI*CINUI
- CINUI = STR*CINUI + STI*CINUR
- CINUR = PTR
- STR = YR(I)
- STI = YI(I)
- YR(I) = CINUR*CSCLR
- YI(I) = CINUI*CSCLR
- 40 CONTINUE
- RETURN
- 50 CONTINUE
- NZ = -1
- IF(NW.EQ.(-2)) NZ=-2
- RETURN
- END
|