github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zasyi.f (about) 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,ZABS,ZDIV,ZEXP,ZMLT,ZSQRT 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, ZABS 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 = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 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 ZSQRT(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 ZEXP(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) THEN 55 FDN = DNU2*DNU2 56 END IF 57 EZR = ZR*8.0D0 58 EZI = ZI*8.0D0 59 C----------------------------------------------------------------------- 60 C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE 61 C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE 62 C EXPANSION FOR THE IMAGINARY PART. 63 C----------------------------------------------------------------------- 64 AEZ = 8.0D0*AZ 65 S = TOL/AEZ 66 JL = INT(SNGL(RL+RL)) + 2 67 P1R = ZEROR 68 P1I = ZEROI 69 IF (ZI.EQ.0.0D0) GO TO 30 70 C----------------------------------------------------------------------- 71 C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF 72 C SIGNIFICANCE WHEN FNU OR N IS LARGE 73 C----------------------------------------------------------------------- 74 INU = INT(SNGL(FNU)) 75 ARG = (FNU-DBLE(FLOAT(INU)))*PI 76 INU = INU + N - IL 77 AK = -DSIN(ARG) 78 BK = DCOS(ARG) 79 IF (ZI.LT.0.0D0) BK = -BK 80 P1R = AK 81 P1I = BK 82 IF (MOD(INU,2).EQ.0) GO TO 30 83 P1R = -P1R 84 P1I = -P1I 85 30 CONTINUE 86 DO 70 K=1,IL 87 SQK = FDN - 1.0D0 88 ATOL = S*DABS(SQK) 89 SGN = 1.0D0 90 CS1R = CONER 91 CS1I = CONEI 92 CS2R = CONER 93 CS2I = CONEI 94 CKR = CONER 95 CKI = CONEI 96 AK = 0.0D0 97 AA = 1.0D0 98 BB = AEZ 99 DKR = EZR 100 DKI = EZI 101 DO 40 J=1,JL 102 CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI) 103 CKR = STR*SQK 104 CKI = STI*SQK 105 CS2R = CS2R + CKR 106 CS2I = CS2I + CKI 107 SGN = -SGN 108 CS1R = CS1R + CKR*SGN 109 CS1I = CS1I + CKI*SGN 110 DKR = DKR + EZR 111 DKI = DKI + EZI 112 AA = AA*DABS(SQK)/BB 113 BB = BB + AEZ 114 AK = AK + 8.0D0 115 SQK = SQK - AK 116 IF (AA.LE.ATOL) THEN 117 GO TO 50 118 END IF 119 40 CONTINUE 120 GO TO 110 121 50 CONTINUE 122 S2R = CS1R 123 S2I = CS1I 124 IF (ZR+ZR.GE.ELIM) GO TO 60 125 TZR = ZR + ZR 126 TZI = ZI + ZI 127 CALL ZEXP(-TZR, -TZI, STR, STI) 128 CALL ZMLT(STR, STI, P1R, P1I, STR, STI) 129 CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI) 130 S2R = S2R + STR 131 S2I = S2I + STI 132 60 CONTINUE 133 FDN = FDN + 8.0D0*DFNU + 4.0D0 134 P1R = -P1R 135 P1I = -P1I 136 M = N - IL + K 137 YR(M) = S2R*AK1R - S2I*AK1I 138 YI(M) = S2R*AK1I + S2I*AK1R 139 70 CONTINUE 140 IF (N.LE.2) RETURN 141 NN = N 142 K = NN - 2 143 AK = DBLE(FLOAT(K)) 144 STR = ZR*RAZ 145 STI = -ZI*RAZ 146 RZR = (STR+STR)*RAZ 147 RZI = (STI+STI)*RAZ 148 IB = 3 149 DO 80 I=IB,NN 150 YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2) 151 YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2) 152 AK = AK - 1.0D0 153 K = K - 1 154 80 CONTINUE 155 IF (KODED.EQ.0) RETURN 156 CALL ZEXP(CZR, CZI, CKR, CKI) 157 DO 90 I=1,NN 158 STR = YR(I)*CKR - YI(I)*CKI 159 YI(I) = YR(I)*CKI + YI(I)*CKR 160 YR(I) = STR 161 90 CONTINUE 162 RETURN 163 100 CONTINUE 164 NZ = -1 165 RETURN 166 110 CONTINUE 167 NZ=-2 168 RETURN 169 END