amos/zwrsk.f

   1       SUBROUTINE ZWRSK(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI,
   2      * TOL, ELIM, ALIM)
   3 C***BEGIN PROLOGUE  ZWRSK
   4 C***REFER TO  ZBESI,ZBESK
   5 C
   6 C     ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
   7 C     NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
   8 C
   9 C***ROUTINES CALLED  D1MACH,ZBKNU,ZRATI,AZABS
  10 C***END PROLOGUE  ZWRSK
  11 C     COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR
  12       DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI,
  13      * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT,
  14      * STI, STR, TOL, YI, YR, ZRI, ZRR, AZABS, D1MACH
  15       INTEGER I, KODE, N, NW, NZ
  16       DIMENSION YR(N), YI(N), CWR(2), CWI(2)
  17 C-----------------------------------------------------------------------
  18 C     I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
  19 C     Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
  20 C     WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
  21 C-----------------------------------------------------------------------
  22       NZ = 0
  23       CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM)
  24       IF (NW.NE.0) GO TO 50
  25       CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL)
  26 C-----------------------------------------------------------------------
  27 C     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
  28 C     R(FNU+J-1,Z)=Y(J),  J=1,...,N
  29 C-----------------------------------------------------------------------
  30       CINUR = 1.0D0
  31       CINUI = 0.0D0
  32       IF (KODE.EQ.1) GO TO 10
  33       CINUR = DCOS(ZRI)
  34       CINUI = DSIN(ZRI)
  35    10 CONTINUE
  36 C-----------------------------------------------------------------------
  37 C     ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
  38 C     THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
  39 C     SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
  40 C     THE RESULT IS ON SCALE.
  41 C-----------------------------------------------------------------------
  42       ACW = AZABS(CWR(2),CWI(2))
  43       ASCLE = 1.0D+3*D1MACH(1)/TOL
  44       CSCLR = 1.0D0
  45       IF (ACW.GT.ASCLE) GO TO 20
  46       CSCLR = 1.0D0/TOL
  47       GO TO 30
  48    20 CONTINUE
  49       ASCLE = 1.0D0/ASCLE
  50       IF (ACW.LT.ASCLE) GO TO 30
  51       CSCLR = TOL
  52    30 CONTINUE
  53       C1R = CWR(1)*CSCLR
  54       C1I = CWI(1)*CSCLR
  55       C2R = CWR(2)*CSCLR
  56       C2I = CWI(2)*CSCLR
  57       STR = YR(1)
  58       STI = YI(1)
  59 C-----------------------------------------------------------------------
  60 C     CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS
  61 C     UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
  62 C-----------------------------------------------------------------------
  63       PTR = STR*C1R - STI*C1I
  64       PTI = STR*C1I + STI*C1R
  65       PTR = PTR + C2R
  66       PTI = PTI + C2I
  67       CTR = ZRR*PTR - ZRI*PTI
  68       CTI = ZRR*PTI + ZRI*PTR
  69       ACT = AZABS(CTR,CTI)
  70       RACT = 1.0D0/ACT
  71       CTR = CTR*RACT
  72       CTI = -CTI*RACT
  73       PTR = CINUR*RACT
  74       PTI = CINUI*RACT
  75       CINUR = PTR*CTR - PTI*CTI
  76       CINUI = PTR*CTI + PTI*CTR
  77       YR(1) = CINUR*CSCLR
  78       YI(1) = CINUI*CSCLR
  79       IF (N.EQ.1) RETURN
  80       DO 40 I=2,N
  81         PTR = STR*CINUR - STI*CINUI
  82         CINUI = STR*CINUI + STI*CINUR
  83         CINUR = PTR
  84         STR = YR(I)
  85         STI = YI(I)
  86         YR(I) = CINUR*CSCLR
  87         YI(I) = CINUI*CSCLR
  88    40 CONTINUE
  89       RETURN
  90    50 CONTINUE
  91       NZ = -1
  92       IF(NW.EQ.(-2)) NZ=-2
  93       RETURN
  94       END