zacai.f 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899
  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(COMPLEX(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