zacai.f 3.9 KB

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