github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zwrsk.f (about)

     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,ZABS
    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, ZABS, 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 = ZABS(CMPLX(CWR(2),CWI(2),kind=KIND(1.0D0)))
    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 = ZABS(CMPLX(CTR,CTI,kind=KIND(1.0D0)))
    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