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

     1        SUBROUTINE ZBINU(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL,
     2       * TOL, ELIM, ALIM)
     3  C***BEGIN PROLOGUE  ZBINU
     4  C***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY
     5  C
     6  C     ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE
     7  C
     8  C***ROUTINES CALLED  ZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK
     9  C***END PROLOGUE  ZBINU
    10        DOUBLE PRECISION ALIM, AZ, CWI, CWR, CYI, CYR, DFNU, ELIM, FNU,
    11       * FNUL, RL, TOL, ZEROI, ZEROR, ZI, ZR, ZABS
    12        INTEGER I, INW, KODE, N, NLAST, NN, NUI, NW, NZ
    13        DIMENSION CYR(N), CYI(N), CWR(2), CWI(2)
    14        DATA ZEROR,ZEROI / 0.0D0, 0.0D0 /
    15  C
    16        NZ = 0
    17        AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0)))
    18        NN = N
    19        DFNU = FNU + DBLE(FLOAT(N-1))
    20        IF (AZ.LE.2.0D0) GO TO 10
    21        IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20
    22     10 CONTINUE
    23  C-----------------------------------------------------------------------
    24  C     POWER SERIES
    25  C-----------------------------------------------------------------------
    26        CALL ZSERI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
    27        INW = IABS(NW)
    28        NZ = NZ + INW
    29        NN = NN - INW
    30        IF (NN.EQ.0) RETURN
    31        IF (NW.GE.0) GO TO 120
    32        DFNU = FNU + DBLE(FLOAT(NN-1))
    33     20 CONTINUE
    34        IF (AZ.LT.RL) GO TO 40
    35        IF (DFNU.LE.1.0D0) GO TO 30
    36        IF (AZ+AZ.LT.DFNU*DFNU) GO TO 50
    37  C-----------------------------------------------------------------------
    38  C     ASYMPTOTIC EXPANSION FOR LARGE Z
    39  C-----------------------------------------------------------------------
    40     30 CONTINUE
    41        CALL ZASYI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, RL, TOL, ELIM,
    42       * ALIM)
    43        IF (NW.LT.0) GO TO 130
    44        GO TO 120
    45     40 CONTINUE
    46        IF (DFNU.LE.1.0D0) GO TO 70
    47     50 CONTINUE
    48  C-----------------------------------------------------------------------
    49  C     OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
    50  C-----------------------------------------------------------------------
    51        CALL ZUOIK(ZR, ZI, FNU, KODE, 1, NN, CYR, CYI, NW, TOL, ELIM,
    52       * ALIM)
    53        IF (NW.LT.0) GO TO 130
    54        NZ = NZ + NW
    55        NN = NN - NW
    56        IF (NN.EQ.0) RETURN
    57        DFNU = FNU+DBLE(FLOAT(NN-1))
    58        IF (DFNU.GT.FNUL) GO TO 110
    59        IF (AZ.GT.FNUL) GO TO 110
    60     60 CONTINUE
    61        IF (AZ.GT.RL) GO TO 80
    62     70 CONTINUE
    63  C-----------------------------------------------------------------------
    64  C     MILLER ALGORITHM NORMALIZED BY THE SERIES
    65  C-----------------------------------------------------------------------
    66        CALL ZMLRI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL)
    67        IF(NW.LT.0) GO TO 130
    68        GO TO 120
    69     80 CONTINUE
    70  C-----------------------------------------------------------------------
    71  C     MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
    72  C-----------------------------------------------------------------------
    73  C-----------------------------------------------------------------------
    74  C     OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
    75  C-----------------------------------------------------------------------
    76        CALL ZUOIK(ZR, ZI, FNU, KODE, 2, 2, CWR, CWI, NW, TOL, ELIM,
    77       * ALIM)
    78        IF (NW.GE.0) GO TO 100
    79        NZ = NN
    80        DO 90 I=1,NN
    81          CYR(I) = ZEROR
    82          CYI(I) = ZEROI
    83     90 CONTINUE
    84        RETURN
    85    100 CONTINUE
    86        IF (NW.GT.0) GO TO 130
    87        CALL ZWRSK(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, CWR, CWI, TOL,
    88       * ELIM, ALIM)
    89        IF (NW.LT.0) GO TO 130
    90        GO TO 120
    91    110 CONTINUE
    92  C-----------------------------------------------------------------------
    93  C     INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
    94  C-----------------------------------------------------------------------
    95        NUI = INT(SNGL(FNUL-DFNU)) + 1
    96        NUI = MAX0(NUI,0)
    97        CALL ZBUNI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, NUI, NLAST, FNUL,
    98       * TOL, ELIM, ALIM)
    99        IF (NW.LT.0) GO TO 130
   100        NZ = NZ + NW
   101        IF (NLAST.EQ.0) GO TO 120
   102        NN = NLAST
   103        GO TO 60
   104    120 CONTINUE
   105        RETURN
   106    130 CONTINUE
   107        NZ = -1
   108        IF(NW.EQ.(-2)) NZ=-2
   109        RETURN
   110        END