amos/zacon.f

   1       SUBROUTINE ZACON(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL,
   2      * TOL, ELIM, ALIM)
   3 C***BEGIN PROLOGUE  ZACON
   4 C***REFER TO  ZBESK,ZBESH
   5 C
   6 C     ZACON 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
  13 C
  14 C***ROUTINES CALLED  ZBINU,ZBKNU,ZS1S2,D1MACH,AZABS,ZMLT
  15 C***END PROLOGUE  ZACON
  16 C     COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
  17 C    *S1,S2,Y,Z,ZN
  18       DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI,
  19      * CKR, CONER, CPN, CSCL, CSCR, CSGNI, CSGNR, CSPNI, CSPNR,
  20      * CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR,
  21      * FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R,
  22      * SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR,
  23      * YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, AZABS
  24       INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
  25       DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3)
  26       DATA PI / 3.14159265358979324D0 /
  27       DATA ZEROR,CONER / 0.0D0,1.0D0 /
  28       NZ = 0
  29       ZNR = -ZR
  30       ZNI = -ZI
  31       NN = N
  32       CALL ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, FNUL, TOL,
  33      * ELIM, ALIM)
  34       IF (NW.LT.0) GO TO 90
  35 C-----------------------------------------------------------------------
  36 C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
  37 C-----------------------------------------------------------------------
  38       NN = MIN0(2,N)
  39       CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
  40       IF (NW.NE.0) GO TO 90
  41       S1R = CYR(1)
  42       S1I = CYI(1)
  43       FMR = DBLE(FLOAT(MR))
  44       SGN = -DSIGN(PI,FMR)
  45       CSGNR = ZEROR
  46       CSGNI = SGN
  47       IF (KODE.EQ.1) GO TO 10
  48       YY = -ZNI
  49       CPN = DCOS(YY)
  50       SPN = DSIN(YY)
  51       CALL ZMLT(CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI)
  52    10 CONTINUE
  53 C-----------------------------------------------------------------------
  54 C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
  55 C     WHEN FNU IS LARGE
  56 C-----------------------------------------------------------------------
  57       INU = INT(SNGL(FNU))
  58       ARG = (FNU-DBLE(FLOAT(INU)))*SGN
  59       CPN = DCOS(ARG)
  60       SPN = DSIN(ARG)
  61       CSPNR = CPN
  62       CSPNI = SPN
  63       IF (MOD(INU,2).EQ.0) GO TO 20
  64       CSPNR = -CSPNR
  65       CSPNI = -CSPNI
  66    20 CONTINUE
  67       IUF = 0
  68       C1R = S1R
  69       C1I = S1I
  70       C2R = YR(1)
  71       C2I = YI(1)
  72       ASCLE = 1.0D+3*D1MACH(1)/TOL
  73       IF (KODE.EQ.1) GO TO 30
  74       CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
  75       NZ = NZ + NW
  76       SC1R = C1R
  77       SC1I = C1I
  78    30 CONTINUE
  79       CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
  80       CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
  81       YR(1) = STR + PTR
  82       YI(1) = STI + PTI
  83       IF (N.EQ.1) RETURN
  84       CSPNR = -CSPNR
  85       CSPNI = -CSPNI
  86       S2R = CYR(2)
  87       S2I = CYI(2)
  88       C1R = S2R
  89       C1I = S2I
  90       C2R = YR(2)
  91       C2I = YI(2)
  92       IF (KODE.EQ.1) GO TO 40
  93       CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
  94       NZ = NZ + NW
  95       SC2R = C1R
  96       SC2I = C1I
  97    40 CONTINUE
  98       CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
  99       CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
 100       YR(2) = STR + PTR
 101       YI(2) = STI + PTI
 102       IF (N.EQ.2) RETURN
 103       CSPNR = -CSPNR
 104       CSPNI = -CSPNI
 105       AZN = AZABS(ZNR,ZNI)
 106       RAZN = 1.0D0/AZN
 107       STR = ZNR*RAZN
 108       STI = -ZNI*RAZN
 109       RZR = (STR+STR)*RAZN
 110       RZI = (STI+STI)*RAZN
 111       FN = FNU + 1.0D0
 112       CKR = FN*RZR
 113       CKI = FN*RZI
 114 C-----------------------------------------------------------------------
 115 C     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
 116 C-----------------------------------------------------------------------
 117       CSCL = 1.0D0/TOL
 118       CSCR = TOL
 119       CSSR(1) = CSCL
 120       CSSR(2) = CONER
 121       CSSR(3) = CSCR
 122       CSRR(1) = CSCR
 123       CSRR(2) = CONER
 124       CSRR(3) = CSCL
 125       BRY(1) = ASCLE
 126       BRY(2) = 1.0D0/ASCLE
 127       BRY(3) = D1MACH(2)
 128       AS2 = AZABS(S2R,S2I)
 129       KFLAG = 2
 130       IF (AS2.GT.BRY(1)) GO TO 50
 131       KFLAG = 1
 132       GO TO 60
 133    50 CONTINUE
 134       IF (AS2.LT.BRY(2)) GO TO 60
 135       KFLAG = 3
 136    60 CONTINUE
 137       BSCLE = BRY(KFLAG)
 138       S1R = S1R*CSSR(KFLAG)
 139       S1I = S1I*CSSR(KFLAG)
 140       S2R = S2R*CSSR(KFLAG)
 141       S2I = S2I*CSSR(KFLAG)
 142       CSR = CSRR(KFLAG)
 143       DO 80 I=3,N
 144         STR = S2R
 145         STI = S2I
 146         S2R = CKR*STR - CKI*STI + S1R
 147         S2I = CKR*STI + CKI*STR + S1I
 148         S1R = STR
 149         S1I = STI
 150         C1R = S2R*CSR
 151         C1I = S2I*CSR
 152         STR = C1R
 153         STI = C1I
 154         C2R = YR(I)
 155         C2I = YI(I)
 156         IF (KODE.EQ.1) GO TO 70
 157         IF (IUF.LT.0) GO TO 70
 158         CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
 159         NZ = NZ + NW
 160         SC1R = SC2R
 161         SC1I = SC2I
 162         SC2R = C1R
 163         SC2I = C1I
 164         IF (IUF.NE.3) GO TO 70
 165         IUF = -4
 166         S1R = SC1R*CSSR(KFLAG)
 167         S1I = SC1I*CSSR(KFLAG)
 168         S2R = SC2R*CSSR(KFLAG)
 169         S2I = SC2I*CSSR(KFLAG)
 170         STR = SC2R
 171         STI = SC2I
 172    70   CONTINUE
 173         PTR = CSPNR*C1R - CSPNI*C1I
 174         PTI = CSPNR*C1I + CSPNI*C1R
 175         YR(I) = PTR + CSGNR*C2R - CSGNI*C2I
 176         YI(I) = PTI + CSGNR*C2I + CSGNI*C2R
 177         CKR = CKR + RZR
 178         CKI = CKI + RZI
 179         CSPNR = -CSPNR
 180         CSPNI = -CSPNI
 181         IF (KFLAG.GE.3) GO TO 80
 182         PTR = DABS(C1R)
 183         PTI = DABS(C1I)
 184         C1M = DMAX1(PTR,PTI)
 185         IF (C1M.LE.BSCLE) GO TO 80
 186         KFLAG = KFLAG + 1
 187         BSCLE = BRY(KFLAG)
 188         S1R = S1R*CSR
 189         S1I = S1I*CSR
 190         S2R = STR
 191         S2I = STI
 192         S1R = S1R*CSSR(KFLAG)
 193         S1I = S1I*CSSR(KFLAG)
 194         S2R = S2R*CSSR(KFLAG)
 195         S2I = S2I*CSSR(KFLAG)
 196         CSR = CSRR(KFLAG)
 197    80 CONTINUE
 198       RETURN
 199    90 CONTINUE
 200       NZ = -1
 201       IF(NW.EQ.(-2)) NZ=-2
 202       RETURN
 203       END