95 lines
3.3 KiB
Fortran
95 lines
3.3 KiB
Fortran
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,AZABS
|
|
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, AZABS, 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 = AZABS(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 = AZABS(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
|