amos/zasyi.f

   1       SUBROUTINE ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
   2      * ALIM)
   3 C***BEGIN PROLOGUE  ZASYI
   4 C***REFER TO  ZBESI,ZBESK
   5 C
   6 C     ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
   7 C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
   8 C     REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
   9 C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
  10 C
  11 C***ROUTINES CALLED  D1MACH,AZABS,ZDIV,AZEXP,ZMLT,AZSQRT
  12 C***END PROLOGUE  ZASYI
  13 C     COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
  14       DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL,
  15      * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
  16      * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I,
  17      * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I,
  18      * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, AZABS
  19       INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
  20       DIMENSION YR(N), YI(N)
  21       DATA PI, RTPI  /3.14159265358979324D0 , 0.159154943091895336D0 /
  22       DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
  23 C
  24       NZ = 0
  25       AZ = AZABS(ZR,ZI)
  26       ARM = 1.0D+3*D1MACH(1)
  27       RTR1 = DSQRT(ARM)
  28       IL = MIN0(2,N)
  29       DFNU = FNU + DBLE(FLOAT(N-IL))
  30 C-----------------------------------------------------------------------
  31 C     OVERFLOW TEST
  32 C-----------------------------------------------------------------------
  33       RAZ = 1.0D0/AZ
  34       STR = ZR*RAZ
  35       STI = -ZI*RAZ
  36       AK1R = RTPI*STR*RAZ
  37       AK1I = RTPI*STI*RAZ
  38       CALL AZSQRT(AK1R, AK1I, AK1R, AK1I)
  39       CZR = ZR
  40       CZI = ZI
  41       IF (KODE.NE.2) GO TO 10
  42       CZR = ZEROR
  43       CZI = ZI
  44    10 CONTINUE
  45       IF (DABS(CZR).GT.ELIM) GO TO 100
  46       DNU2 = DFNU + DFNU
  47       KODED = 1
  48       IF ((DABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20
  49       KODED = 0
  50       CALL AZEXP(CZR, CZI, STR, STI)
  51       CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I)
  52    20 CONTINUE
  53       FDN = 0.0D0
  54       IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
  55       EZR = ZR*8.0D0
  56       EZI = ZI*8.0D0
  57 C-----------------------------------------------------------------------
  58 C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
  59 C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
  60 C     EXPANSION FOR THE IMAGINARY PART.
  61 C-----------------------------------------------------------------------
  62       AEZ = 8.0D0*AZ
  63       S = TOL/AEZ
  64       JL = INT(SNGL(RL+RL)) + 2
  65       P1R = ZEROR
  66       P1I = ZEROI
  67       IF (ZI.EQ.0.0D0) GO TO 30
  68 C-----------------------------------------------------------------------
  69 C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
  70 C     SIGNIFICANCE WHEN FNU OR N IS LARGE
  71 C-----------------------------------------------------------------------
  72       INU = INT(SNGL(FNU))
  73       ARG = (FNU-DBLE(FLOAT(INU)))*PI
  74       INU = INU + N - IL
  75       AK = -DSIN(ARG)
  76       BK = DCOS(ARG)
  77       IF (ZI.LT.0.0D0) BK = -BK
  78       P1R = AK
  79       P1I = BK
  80       IF (MOD(INU,2).EQ.0) GO TO 30
  81       P1R = -P1R
  82       P1I = -P1I
  83    30 CONTINUE
  84       DO 70 K=1,IL
  85         SQK = FDN - 1.0D0
  86         ATOL = S*DABS(SQK)
  87         SGN = 1.0D0
  88         CS1R = CONER
  89         CS1I = CONEI
  90         CS2R = CONER
  91         CS2I = CONEI
  92         CKR = CONER
  93         CKI = CONEI
  94         AK = 0.0D0
  95         AA = 1.0D0
  96         BB = AEZ
  97         DKR = EZR
  98         DKI = EZI
  99         DO 40 J=1,JL
 100           CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
 101           CKR = STR*SQK
 102           CKI = STI*SQK
 103           CS2R = CS2R + CKR
 104           CS2I = CS2I + CKI
 105           SGN = -SGN
 106           CS1R = CS1R + CKR*SGN
 107           CS1I = CS1I + CKI*SGN
 108           DKR = DKR + EZR
 109           DKI = DKI + EZI
 110           AA = AA*DABS(SQK)/BB
 111           BB = BB + AEZ
 112           AK = AK + 8.0D0
 113           SQK = SQK - AK
 114           IF (AA.LE.ATOL) GO TO 50
 115    40   CONTINUE
 116         GO TO 110
 117    50   CONTINUE
 118         S2R = CS1R
 119         S2I = CS1I
 120         IF (ZR+ZR.GE.ELIM) GO TO 60
 121         TZR = ZR + ZR
 122         TZI = ZI + ZI
 123         CALL AZEXP(-TZR, -TZI, STR, STI)
 124         CALL ZMLT(STR, STI, P1R, P1I, STR, STI)
 125         CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI)
 126         S2R = S2R + STR
 127         S2I = S2I + STI
 128    60   CONTINUE
 129         FDN = FDN + 8.0D0*DFNU + 4.0D0
 130         P1R = -P1R
 131         P1I = -P1I
 132         M = N - IL + K
 133         YR(M) = S2R*AK1R - S2I*AK1I
 134         YI(M) = S2R*AK1I + S2I*AK1R
 135    70 CONTINUE
 136       IF (N.LE.2) RETURN
 137       NN = N
 138       K = NN - 2
 139       AK = DBLE(FLOAT(K))
 140       STR = ZR*RAZ
 141       STI = -ZI*RAZ
 142       RZR = (STR+STR)*RAZ
 143       RZI = (STI+STI)*RAZ
 144       IB = 3
 145       DO 80 I=IB,NN
 146         YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
 147         YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
 148         AK = AK - 1.0D0
 149         K = K - 1
 150    80 CONTINUE
 151       IF (KODED.EQ.0) RETURN
 152       CALL AZEXP(CZR, CZI, CKR, CKI)
 153       DO 90 I=1,NN
 154         STR = YR(I)*CKR - YI(I)*CKI
 155         YI(I) = YR(I)*CKI + YI(I)*CKR
 156         YR(I) = STR
 157    90 CONTINUE
 158       RETURN
 159   100 CONTINUE
 160       NZ = -1
 161       RETURN
 162   110 CONTINUE
 163       NZ=-2
 164       RETURN
 165       END