zbinu.f 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110
  1. SUBROUTINE ZBINU(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL,
  2. * TOL, ELIM, ALIM)
  3. C***BEGIN PROLOGUE ZBINU
  4. C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY
  5. C
  6. C ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE
  7. C
  8. C***ROUTINES CALLED ZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK
  9. C***END PROLOGUE ZBINU
  10. DOUBLE PRECISION ALIM, AZ, CWI, CWR, CYI, CYR, DFNU, ELIM, FNU,
  11. * FNUL, RL, TOL, ZEROI, ZEROR, ZI, ZR, ZABS
  12. INTEGER I, INW, KODE, N, NLAST, NN, NUI, NW, NZ
  13. DIMENSION CYR(N), CYI(N), CWR(2), CWI(2)
  14. DATA ZEROR,ZEROI / 0.0D0, 0.0D0 /
  15. C
  16. NZ = 0
  17. AZ = ZABS(COMPLEX(ZR,ZI))
  18. NN = N
  19. DFNU = FNU + DBLE(FLOAT(N-1))
  20. IF (AZ.LE.2.0D0) GO TO 10
  21. IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20
  22. 10 CONTINUE
  23. C-----------------------------------------------------------------------
  24. C POWER SERIES
  25. C-----------------------------------------------------------------------
  26. CALL ZSERI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
  27. INW = IABS(NW)
  28. NZ = NZ + INW
  29. NN = NN - INW
  30. IF (NN.EQ.0) RETURN
  31. IF (NW.GE.0) GO TO 120
  32. DFNU = FNU + DBLE(FLOAT(NN-1))
  33. 20 CONTINUE
  34. IF (AZ.LT.RL) GO TO 40
  35. IF (DFNU.LE.1.0D0) GO TO 30
  36. IF (AZ+AZ.LT.DFNU*DFNU) GO TO 50
  37. C-----------------------------------------------------------------------
  38. C ASYMPTOTIC EXPANSION FOR LARGE Z
  39. C-----------------------------------------------------------------------
  40. 30 CONTINUE
  41. CALL ZASYI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, RL, TOL, ELIM,
  42. * ALIM)
  43. IF (NW.LT.0) GO TO 130
  44. GO TO 120
  45. 40 CONTINUE
  46. IF (DFNU.LE.1.0D0) GO TO 70
  47. 50 CONTINUE
  48. C-----------------------------------------------------------------------
  49. C OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
  50. C-----------------------------------------------------------------------
  51. CALL ZUOIK(ZR, ZI, FNU, KODE, 1, NN, CYR, CYI, NW, TOL, ELIM,
  52. * ALIM)
  53. IF (NW.LT.0) GO TO 130
  54. NZ = NZ + NW
  55. NN = NN - NW
  56. IF (NN.EQ.0) RETURN
  57. DFNU = FNU+DBLE(FLOAT(NN-1))
  58. IF (DFNU.GT.FNUL) GO TO 110
  59. IF (AZ.GT.FNUL) GO TO 110
  60. 60 CONTINUE
  61. IF (AZ.GT.RL) GO TO 80
  62. 70 CONTINUE
  63. C-----------------------------------------------------------------------
  64. C MILLER ALGORITHM NORMALIZED BY THE SERIES
  65. C-----------------------------------------------------------------------
  66. CALL ZMLRI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL)
  67. IF(NW.LT.0) GO TO 130
  68. GO TO 120
  69. 80 CONTINUE
  70. C-----------------------------------------------------------------------
  71. C MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
  72. C-----------------------------------------------------------------------
  73. C-----------------------------------------------------------------------
  74. C OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
  75. C-----------------------------------------------------------------------
  76. CALL ZUOIK(ZR, ZI, FNU, KODE, 2, 2, CWR, CWI, NW, TOL, ELIM,
  77. * ALIM)
  78. IF (NW.GE.0) GO TO 100
  79. NZ = NN
  80. DO 90 I=1,NN
  81. CYR(I) = ZEROR
  82. CYI(I) = ZEROI
  83. 90 CONTINUE
  84. RETURN
  85. 100 CONTINUE
  86. IF (NW.GT.0) GO TO 130
  87. CALL ZWRSK(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, CWR, CWI, TOL,
  88. * ELIM, ALIM)
  89. IF (NW.LT.0) GO TO 130
  90. GO TO 120
  91. 110 CONTINUE
  92. C-----------------------------------------------------------------------
  93. C INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
  94. C-----------------------------------------------------------------------
  95. NUI = INT(SNGL(FNUL-DFNU)) + 1
  96. NUI = MAX0(NUI,0)
  97. CALL ZBUNI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, NUI, NLAST, FNUL,
  98. * TOL, ELIM, ALIM)
  99. IF (NW.LT.0) GO TO 130
  100. NZ = NZ + NW
  101. IF (NLAST.EQ.0) GO TO 120
  102. NN = NLAST
  103. GO TO 60
  104. 120 CONTINUE
  105. RETURN
  106. 130 CONTINUE
  107. NZ = -1
  108. IF(NW.EQ.(-2)) NZ=-2
  109. RETURN
  110. END