zwrsk.f 3.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394
  1. SUBROUTINE ZWRSK(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI,
  2. * TOL, ELIM, ALIM)
  3. C***BEGIN PROLOGUE ZWRSK
  4. C***REFER TO ZBESI,ZBESK
  5. C
  6. C ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
  7. C NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
  8. C
  9. C***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,ZABS
  10. C***END PROLOGUE ZWRSK
  11. C COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR
  12. DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI,
  13. * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT,
  14. * STI, STR, TOL, YI, YR, ZRI, ZRR, ZABS, D1MACH
  15. INTEGER I, KODE, N, NW, NZ
  16. DIMENSION YR(N), YI(N), CWR(2), CWI(2)
  17. C-----------------------------------------------------------------------
  18. C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
  19. C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
  20. C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
  21. C-----------------------------------------------------------------------
  22. NZ = 0
  23. CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM)
  24. IF (NW.NE.0) GO TO 50
  25. CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL)
  26. C-----------------------------------------------------------------------
  27. C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
  28. C R(FNU+J-1,Z)=Y(J), J=1,...,N
  29. C-----------------------------------------------------------------------
  30. CINUR = 1.0D0
  31. CINUI = 0.0D0
  32. IF (KODE.EQ.1) GO TO 10
  33. CINUR = DCOS(ZRI)
  34. CINUI = DSIN(ZRI)
  35. 10 CONTINUE
  36. C-----------------------------------------------------------------------
  37. C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
  38. C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
  39. C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
  40. C THE RESULT IS ON SCALE.
  41. C-----------------------------------------------------------------------
  42. ACW = ZABS(COMPLEX(CWR(2),CWI(2)))
  43. ASCLE = 1.0D+3*D1MACH(1)/TOL
  44. CSCLR = 1.0D0
  45. IF (ACW.GT.ASCLE) GO TO 20
  46. CSCLR = 1.0D0/TOL
  47. GO TO 30
  48. 20 CONTINUE
  49. ASCLE = 1.0D0/ASCLE
  50. IF (ACW.LT.ASCLE) GO TO 30
  51. CSCLR = TOL
  52. 30 CONTINUE
  53. C1R = CWR(1)*CSCLR
  54. C1I = CWI(1)*CSCLR
  55. C2R = CWR(2)*CSCLR
  56. C2I = CWI(2)*CSCLR
  57. STR = YR(1)
  58. STI = YI(1)
  59. C-----------------------------------------------------------------------
  60. C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS
  61. C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
  62. C-----------------------------------------------------------------------
  63. PTR = STR*C1R - STI*C1I
  64. PTI = STR*C1I + STI*C1R
  65. PTR = PTR + C2R
  66. PTI = PTI + C2I
  67. CTR = ZRR*PTR - ZRI*PTI
  68. CTI = ZRR*PTI + ZRI*PTR
  69. ACT = ZABS(COMPLEX(CTR,CTI))
  70. RACT = 1.0D0/ACT
  71. CTR = CTR*RACT
  72. CTI = -CTI*RACT
  73. PTR = CINUR*RACT
  74. PTI = CINUI*RACT
  75. CINUR = PTR*CTR - PTI*CTI
  76. CINUI = PTR*CTI + PTI*CTR
  77. YR(1) = CINUR*CSCLR
  78. YI(1) = CINUI*CSCLR
  79. IF (N.EQ.1) RETURN
  80. DO 40 I=2,N
  81. PTR = STR*CINUR - STI*CINUI
  82. CINUI = STR*CINUI + STI*CINUR
  83. CINUR = PTR
  84. STR = YR(I)
  85. STI = YI(I)
  86. YR(I) = CINUR*CSCLR
  87. YI(I) = CINUI*CSCLR
  88. 40 CONTINUE
  89. RETURN
  90. 50 CONTINUE
  91. NZ = -1
  92. IF(NW.EQ.(-2)) NZ=-2
  93. RETURN
  94. END