166 lines
4.9 KiB
Fortran
166 lines
4.9 KiB
Fortran
SUBROUTINE ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
|
|
* ALIM)
|
|
C***BEGIN PROLOGUE ZASYI
|
|
C***REFER TO ZBESI,ZBESK
|
|
C
|
|
C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
|
|
C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
|
|
C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
|
|
C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
|
|
C
|
|
C***ROUTINES CALLED D1MACH,AZABS,ZDIV,AZEXP,ZMLT,AZSQRT
|
|
C***END PROLOGUE ZASYI
|
|
C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
|
|
DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL,
|
|
* AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
|
|
* CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I,
|
|
* P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I,
|
|
* S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, AZABS
|
|
INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
|
|
DIMENSION YR(N), YI(N)
|
|
DATA PI, RTPI /3.14159265358979324D0 , 0.159154943091895336D0 /
|
|
DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
|
|
C
|
|
NZ = 0
|
|
AZ = AZABS(ZR,ZI)
|
|
ARM = 1.0D+3*D1MACH(1)
|
|
RTR1 = DSQRT(ARM)
|
|
IL = MIN0(2,N)
|
|
DFNU = FNU + DBLE(FLOAT(N-IL))
|
|
C-----------------------------------------------------------------------
|
|
C OVERFLOW TEST
|
|
C-----------------------------------------------------------------------
|
|
RAZ = 1.0D0/AZ
|
|
STR = ZR*RAZ
|
|
STI = -ZI*RAZ
|
|
AK1R = RTPI*STR*RAZ
|
|
AK1I = RTPI*STI*RAZ
|
|
CALL AZSQRT(AK1R, AK1I, AK1R, AK1I)
|
|
CZR = ZR
|
|
CZI = ZI
|
|
IF (KODE.NE.2) GO TO 10
|
|
CZR = ZEROR
|
|
CZI = ZI
|
|
10 CONTINUE
|
|
IF (DABS(CZR).GT.ELIM) GO TO 100
|
|
DNU2 = DFNU + DFNU
|
|
KODED = 1
|
|
IF ((DABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20
|
|
KODED = 0
|
|
CALL AZEXP(CZR, CZI, STR, STI)
|
|
CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I)
|
|
20 CONTINUE
|
|
FDN = 0.0D0
|
|
IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
|
|
EZR = ZR*8.0D0
|
|
EZI = ZI*8.0D0
|
|
C-----------------------------------------------------------------------
|
|
C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
|
|
C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
|
|
C EXPANSION FOR THE IMAGINARY PART.
|
|
C-----------------------------------------------------------------------
|
|
AEZ = 8.0D0*AZ
|
|
S = TOL/AEZ
|
|
JL = INT(SNGL(RL+RL)) + 2
|
|
P1R = ZEROR
|
|
P1I = ZEROI
|
|
IF (ZI.EQ.0.0D0) GO TO 30
|
|
C-----------------------------------------------------------------------
|
|
C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
|
|
C SIGNIFICANCE WHEN FNU OR N IS LARGE
|
|
C-----------------------------------------------------------------------
|
|
INU = INT(SNGL(FNU))
|
|
ARG = (FNU-DBLE(FLOAT(INU)))*PI
|
|
INU = INU + N - IL
|
|
AK = -DSIN(ARG)
|
|
BK = DCOS(ARG)
|
|
IF (ZI.LT.0.0D0) BK = -BK
|
|
P1R = AK
|
|
P1I = BK
|
|
IF (MOD(INU,2).EQ.0) GO TO 30
|
|
P1R = -P1R
|
|
P1I = -P1I
|
|
30 CONTINUE
|
|
DO 70 K=1,IL
|
|
SQK = FDN - 1.0D0
|
|
ATOL = S*DABS(SQK)
|
|
SGN = 1.0D0
|
|
CS1R = CONER
|
|
CS1I = CONEI
|
|
CS2R = CONER
|
|
CS2I = CONEI
|
|
CKR = CONER
|
|
CKI = CONEI
|
|
AK = 0.0D0
|
|
AA = 1.0D0
|
|
BB = AEZ
|
|
DKR = EZR
|
|
DKI = EZI
|
|
DO 40 J=1,JL
|
|
CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
|
|
CKR = STR*SQK
|
|
CKI = STI*SQK
|
|
CS2R = CS2R + CKR
|
|
CS2I = CS2I + CKI
|
|
SGN = -SGN
|
|
CS1R = CS1R + CKR*SGN
|
|
CS1I = CS1I + CKI*SGN
|
|
DKR = DKR + EZR
|
|
DKI = DKI + EZI
|
|
AA = AA*DABS(SQK)/BB
|
|
BB = BB + AEZ
|
|
AK = AK + 8.0D0
|
|
SQK = SQK - AK
|
|
IF (AA.LE.ATOL) GO TO 50
|
|
40 CONTINUE
|
|
GO TO 110
|
|
50 CONTINUE
|
|
S2R = CS1R
|
|
S2I = CS1I
|
|
IF (ZR+ZR.GE.ELIM) GO TO 60
|
|
TZR = ZR + ZR
|
|
TZI = ZI + ZI
|
|
CALL AZEXP(-TZR, -TZI, STR, STI)
|
|
CALL ZMLT(STR, STI, P1R, P1I, STR, STI)
|
|
CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI)
|
|
S2R = S2R + STR
|
|
S2I = S2I + STI
|
|
60 CONTINUE
|
|
FDN = FDN + 8.0D0*DFNU + 4.0D0
|
|
P1R = -P1R
|
|
P1I = -P1I
|
|
M = N - IL + K
|
|
YR(M) = S2R*AK1R - S2I*AK1I
|
|
YI(M) = S2R*AK1I + S2I*AK1R
|
|
70 CONTINUE
|
|
IF (N.LE.2) RETURN
|
|
NN = N
|
|
K = NN - 2
|
|
AK = DBLE(FLOAT(K))
|
|
STR = ZR*RAZ
|
|
STI = -ZI*RAZ
|
|
RZR = (STR+STR)*RAZ
|
|
RZI = (STI+STI)*RAZ
|
|
IB = 3
|
|
DO 80 I=IB,NN
|
|
YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
|
|
YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
|
|
AK = AK - 1.0D0
|
|
K = K - 1
|
|
80 CONTINUE
|
|
IF (KODED.EQ.0) RETURN
|
|
CALL AZEXP(CZR, CZI, CKR, CKI)
|
|
DO 90 I=1,NN
|
|
STR = YR(I)*CKR - YI(I)*CKI
|
|
YI(I) = YR(I)*CKI + YI(I)*CKR
|
|
YR(I) = STR
|
|
90 CONTINUE
|
|
RETURN
|
|
100 CONTINUE
|
|
NZ = -1
|
|
RETURN
|
|
110 CONTINUE
|
|
NZ=-2
|
|
RETURN
|
|
END
|