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

     1        SUBROUTINE ZACAI(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL,
     2       * ELIM, ALIM)
     3  C***BEGIN PROLOGUE  ZACAI
     4  C***REFER TO  ZAIRY
     5  C
     6  C     ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
     7  C
     8  C         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
     9  C                 MP=PI*MR*CMPLX(0.0,1.0)
    10  C
    11  C     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
    12  C     HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
    13  C     ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
    14  C     RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
    15  C     IS CALLED FROM ZAIRY.
    16  C
    17  C***ROUTINES CALLED  ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,ZABS
    18  C***END PROLOGUE  ZACAI
    19  C     COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY
    20        DOUBLE PRECISION ALIM, ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR,
    21       * CSPNI, C1R, C1I, C2R, C2I, CYR, CYI, DFNU, ELIM, FMR, FNU, PI,
    22       * RL, SGN, TOL, YY, YR, YI, ZR, ZI, ZNR, ZNI, D1MACH, ZABS
    23        INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ
    24        DIMENSION YR(N), YI(N), CYR(2), CYI(2)
    25        DATA PI / 3.14159265358979324D0 /
    26        NZ = 0
    27        ZNR = -ZR
    28        ZNI = -ZI
    29        AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0)))
    30        NN = N
    31        DFNU = FNU + DBLE(FLOAT(N-1))
    32        IF (AZ.LE.2.0D0) GO TO 10
    33        IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20
    34     10 CONTINUE
    35  C-----------------------------------------------------------------------
    36  C     POWER SERIES FOR THE I FUNCTION
    37  C-----------------------------------------------------------------------
    38        CALL ZSERI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL, ELIM, ALIM)
    39        GO TO 40
    40     20 CONTINUE
    41        IF (AZ.LT.RL) GO TO 30
    42  C-----------------------------------------------------------------------
    43  C     ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
    44  C-----------------------------------------------------------------------
    45        CALL ZASYI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM,
    46       * ALIM)
    47        IF (NW.LT.0) GO TO 80
    48        GO TO 40
    49     30 CONTINUE
    50  C-----------------------------------------------------------------------
    51  C     MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
    52  C-----------------------------------------------------------------------
    53        CALL ZMLRI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL)
    54        IF(NW.LT.0) GO TO 80
    55     40 CONTINUE
    56  C-----------------------------------------------------------------------
    57  C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
    58  C-----------------------------------------------------------------------
    59        CALL ZBKNU(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM)
    60        IF (NW.NE.0) GO TO 80
    61        FMR = DBLE(FLOAT(MR))
    62        SGN = -DSIGN(PI,FMR)
    63        CSGNR = 0.0D0
    64        CSGNI = SGN
    65        IF (KODE.EQ.1) GO TO 50
    66        YY = -ZNI
    67        CSGNR = -CSGNI*DSIN(YY)
    68        CSGNI = CSGNI*DCOS(YY)
    69     50 CONTINUE
    70  C-----------------------------------------------------------------------
    71  C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
    72  C     WHEN FNU IS LARGE
    73  C-----------------------------------------------------------------------
    74        INU = INT(SNGL(FNU))
    75        ARG = (FNU-DBLE(FLOAT(INU)))*SGN
    76        CSPNR = DCOS(ARG)
    77        CSPNI = DSIN(ARG)
    78        IF (MOD(INU,2).EQ.0) GO TO 60
    79        CSPNR = -CSPNR
    80        CSPNI = -CSPNI
    81     60 CONTINUE
    82        C1R = CYR(1)
    83        C1I = CYI(1)
    84        C2R = YR(1)
    85        C2I = YI(1)
    86        IF (KODE.EQ.1) GO TO 70
    87        IUF = 0
    88        ASCLE = 1.0D+3*D1MACH(1)/TOL
    89        CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
    90        NZ = NZ + NW
    91     70 CONTINUE
    92        YR(1) = CSPNR*C1R - CSPNI*C1I + CSGNR*C2R - CSGNI*C2I
    93        YI(1) = CSPNR*C1I + CSPNI*C1R + CSGNR*C2I + CSGNI*C2R
    94        RETURN
    95     80 CONTINUE
    96        NZ = -1
    97        IF(NW.EQ.(-2)) NZ=-2
    98        RETURN
    99        END