Mercurial > octave-nkf

annotate liboctave/cruft/amos/zacai.f @ 20595:c1a6c31ac29a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
eliminate more simple uses of error_state * ov-classdef.cc: Eliminate simple uses of error_state.
author John W. Eaton <jwe@octave.org>
date Tue, 06 Oct 2015 00:20:02 -0400
parents 648dabbb4c6b
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3217
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
1 SUBROUTINE ZACAI(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
2 * ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
3 C***BEGIN PROLOGUE ZACAI
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
4 C***REFER TO ZAIRY
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
5 C
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
6 C ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
7 C
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
8 C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
9 C MP=PI*MR*CMPLX(0.0,1.0)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
10 C
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
11 C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
12 C HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
13 C ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
14 C RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
15 C IS CALLED FROM ZAIRY.
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
16 C
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
17 C***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,XZABS
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
18 C***END PROLOGUE ZACAI
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
19 C COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
20 DOUBLE PRECISION ALIM, ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
21 * CSPNI, C1R, C1I, C2R, C2I, CYR, CYI, DFNU, ELIM, FMR, FNU, PI,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
22 * RL, SGN, TOL, YY, YR, YI, ZR, ZI, ZNR, ZNI, D1MACH, XZABS
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
23 INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
24 DIMENSION YR(N), YI(N), CYR(2), CYI(2)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
25 DATA PI / 3.14159265358979324D0 /
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
26 NZ = 0
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
27 ZNR = -ZR
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
28 ZNI = -ZI
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
29 AZ = XZABS(ZR,ZI)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
30 NN = N
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
31 DFNU = FNU + DBLE(FLOAT(N-1))
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
32 IF (AZ.LE.2.0D0) GO TO 10
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
33 IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
34 10 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
35 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
36 C POWER SERIES FOR THE I FUNCTION
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
37 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
38 CALL ZSERI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL, ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
39 GO TO 40
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
40 20 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
41 IF (AZ.LT.RL) GO TO 30
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
42 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
43 C ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
44 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
45 CALL ZASYI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
46 * ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
47 IF (NW.LT.0) GO TO 80
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
48 GO TO 40
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
49 30 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
50 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
51 C MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
52 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
53 CALL ZMLRI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
54 IF(NW.LT.0) GO TO 80
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
55 40 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
56 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
57 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
58 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
59 CALL ZBKNU(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
60 IF (NW.NE.0) GO TO 80
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
61 FMR = DBLE(FLOAT(MR))
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
62 SGN = -DSIGN(PI,FMR)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
63 CSGNR = 0.0D0
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
64 CSGNI = SGN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
65 IF (KODE.EQ.1) GO TO 50
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
66 YY = -ZNI
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
67 CSGNR = -CSGNI*DSIN(YY)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
68 CSGNI = CSGNI*DCOS(YY)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
69 50 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
70 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
71 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
72 C WHEN FNU IS LARGE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
73 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
74 INU = INT(SNGL(FNU))
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
75 ARG = (FNU-DBLE(FLOAT(INU)))*SGN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
76 CSPNR = DCOS(ARG)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
77 CSPNI = DSIN(ARG)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
78 IF (MOD(INU,2).EQ.0) GO TO 60
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
79 CSPNR = -CSPNR
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
80 CSPNI = -CSPNI
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
81 60 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
82 C1R = CYR(1)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
83 C1I = CYI(1)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
84 C2R = YR(1)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
85 C2I = YI(1)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
86 IF (KODE.EQ.1) GO TO 70
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
87 IUF = 0
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
88 ASCLE = 1.0D+3*D1MACH(1)/TOL
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
89 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
90 NZ = NZ + NW
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
91 70 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
92 YR(1) = CSPNR*C1R - CSPNI*C1I + CSGNR*C2R - CSGNI*C2I
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
93 YI(1) = CSPNR*C1I + CSPNI*C1R + CSGNR*C2I + CSGNI*C2R
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
94 RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
95 80 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
96 NZ = -1
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
97 IF(NW.EQ.(-2)) NZ=-2
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
98 RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
99 END