1 SUBROUTINE ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, 2 * ALIM) 3C***BEGIN PROLOGUE ZASYI 4C***REFER TO ZBESI,ZBESK 5C 6C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY 7C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE 8C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. 9C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1. 10C 11C***ROUTINES CALLED D1MACH,AZABS,ZDIV,AZEXP,ZMLT,AZSQRT 12C***END PROLOGUE ZASYI 13C 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 / 23C 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)) 30C----------------------------------------------------------------------- 31C OVERFLOW TEST 32C----------------------------------------------------------------------- 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 57C----------------------------------------------------------------------- 58C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE 59C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE 60C EXPANSION FOR THE IMAGINARY PART. 61C----------------------------------------------------------------------- 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 68C----------------------------------------------------------------------- 69C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF 70C SIGNIFICANCE WHEN FNU OR N IS LARGE 71C----------------------------------------------------------------------- 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 166