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