]>
Commit | Line | Data |
---|---|---|
3820ca8e | 1 | |
2 | CDECK ID>, HWSTAB. | |
3 | ||
4 | *CMZ :- -26/04/91 11.11.56 by Bryan Webber | |
5 | ||
6 | *-- Author : Adapted by Bryan Webber | |
7 | ||
8 | C----------------------------------------------------------------------- | |
9 | ||
10 | FUNCTION HWSTAB(F,AFUN,NN,X,MM) | |
11 | ||
12 | C----------------------------------------------------------------------- | |
13 | ||
14 | C MODIFIED CERN INTERPOLATION ROUTINE DIVDIF | |
15 | ||
16 | C LIKE HWUTAB BUT USES FUNCTION AFUN IN PLACE OF ARRAY A | |
17 | ||
18 | C----------------------------------------------------------------------- | |
19 | ||
20 | IMPLICIT NONE | |
21 | ||
22 | INTEGER NN,MM,MMAX,N,M,MPLUS,IX,IY,MID,NPTS,IP,I,J,L,ISUB | |
23 | ||
24 | DOUBLE PRECISION HWSTAB,AFUN,SUM,X,F(NN),T(20),D(20) | |
25 | ||
26 | LOGICAL EXTRA | |
27 | ||
28 | EXTERNAL AFUN | |
29 | ||
30 | DATA MMAX/10/ | |
31 | ||
32 | N=NN | |
33 | ||
34 | M=MIN(MM,MMAX,N-1) | |
35 | ||
36 | MPLUS=M+1 | |
37 | ||
38 | IX=0 | |
39 | ||
40 | IY=N+1 | |
41 | ||
42 | IF (AFUN(1).GT.AFUN(N)) GOTO 94 | |
43 | ||
44 | 91 MID=(IX+IY)/2 | |
45 | ||
46 | IF (X.GE.AFUN(MID)) GOTO 92 | |
47 | ||
48 | IY=MID | |
49 | ||
50 | GOTO 93 | |
51 | ||
52 | 92 IX=MID | |
53 | ||
54 | 93 IF (IY-IX.GT.1) GOTO 91 | |
55 | ||
56 | GOTO 97 | |
57 | ||
58 | 94 MID=(IX+IY)/2 | |
59 | ||
60 | IF (X.LE.AFUN(MID)) GOTO 95 | |
61 | ||
62 | IY=MID | |
63 | ||
64 | GOTO 96 | |
65 | ||
66 | 95 IX=MID | |
67 | ||
68 | 96 IF (IY-IX.GT.1) GOTO 94 | |
69 | ||
70 | 97 NPTS=M+2-MOD(M,2) | |
71 | ||
72 | IP=0 | |
73 | ||
74 | L=0 | |
75 | ||
76 | GOTO 99 | |
77 | ||
78 | 98 L=-L | |
79 | ||
80 | IF (L.GE.0) L=L+1 | |
81 | ||
82 | 99 ISUB=IX+L | |
83 | ||
84 | IF ((1.LE.ISUB).AND.(ISUB.LE.N)) GOTO 100 | |
85 | ||
86 | NPTS=MPLUS | |
87 | ||
88 | GOTO 101 | |
89 | ||
90 | 100 IP=IP+1 | |
91 | ||
92 | T(IP)=AFUN(ISUB) | |
93 | ||
94 | D(IP)=F(ISUB) | |
95 | ||
96 | 101 IF (IP.LT.NPTS) GOTO 98 | |
97 | ||
98 | EXTRA=NPTS.NE.MPLUS | |
99 | ||
100 | DO 14 L=1,M | |
101 | ||
102 | IF (.NOT.EXTRA) GOTO 12 | |
103 | ||
104 | ISUB=MPLUS-L | |
105 | ||
106 | D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB)) | |
107 | ||
108 | 12 I=MPLUS | |
109 | ||
110 | DO 13 J=L,M | |
111 | ||
112 | ISUB=I-L | |
113 | ||
114 | D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB)) | |
115 | ||
116 | I=I-1 | |
117 | ||
118 | 13 CONTINUE | |
119 | ||
120 | 14 CONTINUE | |
121 | ||
122 | SUM=D(MPLUS) | |
123 | ||
124 | IF (EXTRA) SUM=0.5*(SUM+D(M+2)) | |
125 | ||
126 | J=M | |
127 | ||
128 | DO 15 L=1,M | |
129 | ||
130 | SUM=D(J)+(X-T(J))*SUM | |
131 | ||
132 | J=J-1 | |
133 | ||
134 | 15 CONTINUE | |
135 | ||
136 | HWSTAB=SUM | |
137 | ||
138 | END |