amos/zacai.f

   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,AZABS
  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, AZABS
  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 = AZABS(ZR,ZI)
  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