Mercurial > hg > octave-nkf

annotate libcruft/amos/zbinu.f @ 12583:bb29b58e650c release-3-4-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
abandon release-3-4-x branch in favor of workflow using stable and default branches and merging stable to default periodically
author John W. Eaton <jwe@octave.org>
date Fri, 08 Apr 2011 09:06:04 -0400
parents 8b0cb8f79fdc
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 ZBINU(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
2 * TOL, ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
3 C***BEGIN PROLOGUE ZBINU
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
4 C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY
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 ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE
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***ROUTINES CALLED XZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
9 C***END PROLOGUE ZBINU
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
10 DOUBLE PRECISION ALIM, AZ, CWI, CWR, CYI, CYR, DFNU, ELIM, FNU,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
11 * FNUL, RL, TOL, ZEROI, ZEROR, ZI, ZR, XZABS
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
12 INTEGER I, INW, KODE, N, NLAST, NN, NUI, NW, NZ
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
13 DIMENSION CYR(N), CYI(N), CWR(2), CWI(2)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
14 DATA ZEROR,ZEROI / 0.0D0, 0.0D0 /
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
15 C
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
16 NZ = 0
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
17 AZ = XZABS(ZR,ZI)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
18 NN = N
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
19 DFNU = FNU + DBLE(FLOAT(N-1))
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
20 IF (AZ.LE.2.0D0) GO TO 10
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
21 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
22 10 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
23 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
24 C POWER SERIES
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
25 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
26 CALL ZSERI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
27 INW = IABS(NW)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
28 NZ = NZ + INW
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
29 NN = NN - INW
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
30 IF (NN.EQ.0) RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
31 IF (NW.GE.0) GO TO 120
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
32 DFNU = FNU + DBLE(FLOAT(NN-1))
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
33 20 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
34 IF (AZ.LT.RL) GO TO 40
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
35 IF (DFNU.LE.1.0D0) GO TO 30
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
36 IF (AZ+AZ.LT.DFNU*DFNU) GO TO 50
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 C ASYMPTOTIC EXPANSION FOR LARGE Z
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
39 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
40 30 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
41 CALL ZASYI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, RL, TOL, ELIM,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
42 * ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
43 IF (NW.LT.0) GO TO 130
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
44 GO TO 120
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
45 40 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
46 IF (DFNU.LE.1.0D0) GO TO 70
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
47 50 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
48 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
49 C OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
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 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, NN, CYR, CYI, NW, TOL, ELIM,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
52 * ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
53 IF (NW.LT.0) GO TO 130
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
54 NZ = NZ + NW
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
55 NN = NN - NW
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
56 IF (NN.EQ.0) RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
57 DFNU = FNU+DBLE(FLOAT(NN-1))
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
58 IF (DFNU.GT.FNUL) GO TO 110
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
59 IF (AZ.GT.FNUL) GO TO 110
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
60 60 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
61 IF (AZ.GT.RL) GO TO 80
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
62 70 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
63 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
64 C MILLER ALGORITHM NORMALIZED BY THE SERIES
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
65 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
66 CALL ZMLRI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
67 IF(NW.LT.0) GO TO 130
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
68 GO TO 120
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
69 80 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 MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
72 C-----------------------------------------------------------------------
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 C OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
75 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
76 CALL ZUOIK(ZR, ZI, FNU, KODE, 2, 2, CWR, CWI, NW, TOL, ELIM,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
77 * ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
78 IF (NW.GE.0) GO TO 100
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
79 NZ = NN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
80 DO 90 I=1,NN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
81 CYR(I) = ZEROR
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
82 CYI(I) = ZEROI
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
83 90 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
84 RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
85 100 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
86 IF (NW.GT.0) GO TO 130
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
87 CALL ZWRSK(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, CWR, CWI, TOL,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
88 * ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
89 IF (NW.LT.0) GO TO 130
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
90 GO TO 120
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
91 110 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
92 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
93 C INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
94 C-----------------------------------------------------------------------
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
95 NUI = INT(SNGL(FNUL-DFNU)) + 1
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
96 NUI = MAX0(NUI,0)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
97 CALL ZBUNI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, NUI, NLAST, FNUL,
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
98 * TOL, ELIM, ALIM)
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
99 IF (NW.LT.0) GO TO 130
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
100 NZ = NZ + NW
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
101 IF (NLAST.EQ.0) GO TO 120
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
102 NN = NLAST
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
103 GO TO 60
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
104 120 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
105 RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
106 130 CONTINUE
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
107 NZ = -1
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
108 IF(NW.EQ.(-2)) NZ=-2
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
109 RETURN
8b0cb8f79fdc [project @ 1998-11-11 23:32:20 by jwe]
jwe
parents:
diff changeset
110 END