forked from JuliaLang/julia
/
zacai.f
99 lines (99 loc) · 3.62 KB
/
zacai.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
SUBROUTINE ZACAI(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL,
* ELIM, ALIM)
C***BEGIN PROLOGUE ZACAI
C***REFER TO ZAIRY
C
C ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
C
C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
C MP=PI*MR*CMPLX(0.0,1.0)
C
C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
C HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
C ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
C RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
C IS CALLED FROM ZAIRY.
C
C***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,ZABS
C***END PROLOGUE ZACAI
C COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY
DOUBLE PRECISION ALIM, ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR,
* CSPNI, C1R, C1I, C2R, C2I, CYR, CYI, DFNU, ELIM, FMR, FNU, PI,
* RL, SGN, TOL, YY, YR, YI, ZR, ZI, ZNR, ZNI, D1MACH, ZABS
INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ
DIMENSION YR(N), YI(N), CYR(2), CYI(2)
DATA PI / 3.14159265358979324D0 /
NZ = 0
ZNR = -ZR
ZNI = -ZI
AZ = ZABS(COMPLEX(ZR,ZI))
NN = N
DFNU = FNU + DBLE(FLOAT(N-1))
IF (AZ.LE.2.0D0) GO TO 10
IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20
10 CONTINUE
C-----------------------------------------------------------------------
C POWER SERIES FOR THE I FUNCTION
C-----------------------------------------------------------------------
CALL ZSERI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL, ELIM, ALIM)
GO TO 40
20 CONTINUE
IF (AZ.LT.RL) GO TO 30
C-----------------------------------------------------------------------
C ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
C-----------------------------------------------------------------------
CALL ZASYI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM,
* ALIM)
IF (NW.LT.0) GO TO 80
GO TO 40
30 CONTINUE
C-----------------------------------------------------------------------
C MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
C-----------------------------------------------------------------------
CALL ZMLRI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL)
IF(NW.LT.0) GO TO 80
40 CONTINUE
C-----------------------------------------------------------------------
C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
C-----------------------------------------------------------------------
CALL ZBKNU(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM)
IF (NW.NE.0) GO TO 80
FMR = DBLE(FLOAT(MR))
SGN = -DSIGN(PI,FMR)
CSGNR = 0.0D0
CSGNI = SGN
IF (KODE.EQ.1) GO TO 50
YY = -ZNI
CSGNR = -CSGNI*DSIN(YY)
CSGNI = CSGNI*DCOS(YY)
50 CONTINUE
C-----------------------------------------------------------------------
C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
C WHEN FNU IS LARGE
C-----------------------------------------------------------------------
INU = INT(SNGL(FNU))
ARG = (FNU-DBLE(FLOAT(INU)))*SGN
CSPNR = DCOS(ARG)
CSPNI = DSIN(ARG)
IF (MOD(INU,2).EQ.0) GO TO 60
CSPNR = -CSPNR
CSPNI = -CSPNI
60 CONTINUE
C1R = CYR(1)
C1I = CYI(1)
C2R = YR(1)
C2I = YI(1)
IF (KODE.EQ.1) GO TO 70
IUF = 0
ASCLE = 1.0D+3*D1MACH(1)/TOL
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
NZ = NZ + NW
70 CONTINUE
YR(1) = CSPNR*C1R - CSPNI*C1I + CSGNR*C2R - CSGNI*C2I
YI(1) = CSPNR*C1I + CSPNI*C1R + CSGNR*C2I + CSGNI*C2R
RETURN
80 CONTINUE
NZ = -1
IF(NW.EQ.(-2)) NZ=-2
RETURN
END