204 lines
5.8 KiB
Fortran
204 lines
5.8 KiB
Fortran
SUBROUTINE ZACON(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL,
|
|
* TOL, ELIM, ALIM)
|
|
C***BEGIN PROLOGUE ZACON
|
|
C***REFER TO ZBESK,ZBESH
|
|
C
|
|
C ZACON 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
|
|
C
|
|
C***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,AZABS,ZMLT
|
|
C***END PROLOGUE ZACON
|
|
C COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
|
|
C *S1,S2,Y,Z,ZN
|
|
DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI,
|
|
* CKR, CONER, CPN, CSCL, CSCR, CSGNI, CSGNR, CSPNI, CSPNR,
|
|
* CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR,
|
|
* FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R,
|
|
* SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR,
|
|
* YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, AZABS
|
|
INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
|
|
DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3)
|
|
DATA PI / 3.14159265358979324D0 /
|
|
DATA ZEROR,CONER / 0.0D0,1.0D0 /
|
|
NZ = 0
|
|
ZNR = -ZR
|
|
ZNI = -ZI
|
|
NN = N
|
|
CALL ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, FNUL, TOL,
|
|
* ELIM, ALIM)
|
|
IF (NW.LT.0) GO TO 90
|
|
C-----------------------------------------------------------------------
|
|
C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
|
|
C-----------------------------------------------------------------------
|
|
NN = MIN0(2,N)
|
|
CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
|
|
IF (NW.NE.0) GO TO 90
|
|
S1R = CYR(1)
|
|
S1I = CYI(1)
|
|
FMR = DBLE(FLOAT(MR))
|
|
SGN = -DSIGN(PI,FMR)
|
|
CSGNR = ZEROR
|
|
CSGNI = SGN
|
|
IF (KODE.EQ.1) GO TO 10
|
|
YY = -ZNI
|
|
CPN = DCOS(YY)
|
|
SPN = DSIN(YY)
|
|
CALL ZMLT(CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI)
|
|
10 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
|
|
CPN = DCOS(ARG)
|
|
SPN = DSIN(ARG)
|
|
CSPNR = CPN
|
|
CSPNI = SPN
|
|
IF (MOD(INU,2).EQ.0) GO TO 20
|
|
CSPNR = -CSPNR
|
|
CSPNI = -CSPNI
|
|
20 CONTINUE
|
|
IUF = 0
|
|
C1R = S1R
|
|
C1I = S1I
|
|
C2R = YR(1)
|
|
C2I = YI(1)
|
|
ASCLE = 1.0D+3*D1MACH(1)/TOL
|
|
IF (KODE.EQ.1) GO TO 30
|
|
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
|
|
NZ = NZ + NW
|
|
SC1R = C1R
|
|
SC1I = C1I
|
|
30 CONTINUE
|
|
CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
|
|
CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
|
|
YR(1) = STR + PTR
|
|
YI(1) = STI + PTI
|
|
IF (N.EQ.1) RETURN
|
|
CSPNR = -CSPNR
|
|
CSPNI = -CSPNI
|
|
S2R = CYR(2)
|
|
S2I = CYI(2)
|
|
C1R = S2R
|
|
C1I = S2I
|
|
C2R = YR(2)
|
|
C2I = YI(2)
|
|
IF (KODE.EQ.1) GO TO 40
|
|
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
|
|
NZ = NZ + NW
|
|
SC2R = C1R
|
|
SC2I = C1I
|
|
40 CONTINUE
|
|
CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
|
|
CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
|
|
YR(2) = STR + PTR
|
|
YI(2) = STI + PTI
|
|
IF (N.EQ.2) RETURN
|
|
CSPNR = -CSPNR
|
|
CSPNI = -CSPNI
|
|
AZN = AZABS(ZNR,ZNI)
|
|
RAZN = 1.0D0/AZN
|
|
STR = ZNR*RAZN
|
|
STI = -ZNI*RAZN
|
|
RZR = (STR+STR)*RAZN
|
|
RZI = (STI+STI)*RAZN
|
|
FN = FNU + 1.0D0
|
|
CKR = FN*RZR
|
|
CKI = FN*RZI
|
|
C-----------------------------------------------------------------------
|
|
C SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
|
|
C-----------------------------------------------------------------------
|
|
CSCL = 1.0D0/TOL
|
|
CSCR = TOL
|
|
CSSR(1) = CSCL
|
|
CSSR(2) = CONER
|
|
CSSR(3) = CSCR
|
|
CSRR(1) = CSCR
|
|
CSRR(2) = CONER
|
|
CSRR(3) = CSCL
|
|
BRY(1) = ASCLE
|
|
BRY(2) = 1.0D0/ASCLE
|
|
BRY(3) = D1MACH(2)
|
|
AS2 = AZABS(S2R,S2I)
|
|
KFLAG = 2
|
|
IF (AS2.GT.BRY(1)) GO TO 50
|
|
KFLAG = 1
|
|
GO TO 60
|
|
50 CONTINUE
|
|
IF (AS2.LT.BRY(2)) GO TO 60
|
|
KFLAG = 3
|
|
60 CONTINUE
|
|
BSCLE = BRY(KFLAG)
|
|
S1R = S1R*CSSR(KFLAG)
|
|
S1I = S1I*CSSR(KFLAG)
|
|
S2R = S2R*CSSR(KFLAG)
|
|
S2I = S2I*CSSR(KFLAG)
|
|
CSR = CSRR(KFLAG)
|
|
DO 80 I=3,N
|
|
STR = S2R
|
|
STI = S2I
|
|
S2R = CKR*STR - CKI*STI + S1R
|
|
S2I = CKR*STI + CKI*STR + S1I
|
|
S1R = STR
|
|
S1I = STI
|
|
C1R = S2R*CSR
|
|
C1I = S2I*CSR
|
|
STR = C1R
|
|
STI = C1I
|
|
C2R = YR(I)
|
|
C2I = YI(I)
|
|
IF (KODE.EQ.1) GO TO 70
|
|
IF (IUF.LT.0) GO TO 70
|
|
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
|
|
NZ = NZ + NW
|
|
SC1R = SC2R
|
|
SC1I = SC2I
|
|
SC2R = C1R
|
|
SC2I = C1I
|
|
IF (IUF.NE.3) GO TO 70
|
|
IUF = -4
|
|
S1R = SC1R*CSSR(KFLAG)
|
|
S1I = SC1I*CSSR(KFLAG)
|
|
S2R = SC2R*CSSR(KFLAG)
|
|
S2I = SC2I*CSSR(KFLAG)
|
|
STR = SC2R
|
|
STI = SC2I
|
|
70 CONTINUE
|
|
PTR = CSPNR*C1R - CSPNI*C1I
|
|
PTI = CSPNR*C1I + CSPNI*C1R
|
|
YR(I) = PTR + CSGNR*C2R - CSGNI*C2I
|
|
YI(I) = PTI + CSGNR*C2I + CSGNI*C2R
|
|
CKR = CKR + RZR
|
|
CKI = CKI + RZI
|
|
CSPNR = -CSPNR
|
|
CSPNI = -CSPNI
|
|
IF (KFLAG.GE.3) GO TO 80
|
|
PTR = DABS(C1R)
|
|
PTI = DABS(C1I)
|
|
C1M = DMAX1(PTR,PTI)
|
|
IF (C1M.LE.BSCLE) GO TO 80
|
|
KFLAG = KFLAG + 1
|
|
BSCLE = BRY(KFLAG)
|
|
S1R = S1R*CSR
|
|
S1I = S1I*CSR
|
|
S2R = STR
|
|
S2I = STI
|
|
S1R = S1R*CSSR(KFLAG)
|
|
S1I = S1I*CSSR(KFLAG)
|
|
S2R = S2R*CSSR(KFLAG)
|
|
S2I = S2I*CSSR(KFLAG)
|
|
CSR = CSRR(KFLAG)
|
|
80 CONTINUE
|
|
RETURN
|
|
90 CONTINUE
|
|
NZ = -1
|
|
IF(NW.EQ.(-2)) NZ=-2
|
|
RETURN
|
|
END
|