src/numerical/slatec/fortran/zwrsk.f

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