[u/mrichter/AliRoot.git] / HERWIG / src / hwstab.f


CDECK  ID>, HWSTAB.

*CMZ :-        -26/04/91  11.11.56  by  Bryan Webber

*-- Author :    Adapted by Bryan Webber

C-----------------------------------------------------------------------

      FUNCTION HWSTAB(F,AFUN,NN,X,MM)

C-----------------------------------------------------------------------

C     MODIFIED CERN INTERPOLATION ROUTINE DIVDIF

C     LIKE HWUTAB BUT USES FUNCTION AFUN IN PLACE OF ARRAY A

C-----------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER NN,MM,MMAX,N,M,MPLUS,IX,IY,MID,NPTS,IP,I,J,L,ISUB

      DOUBLE PRECISION HWSTAB,AFUN,SUM,X,F(NN),T(20),D(20)

      LOGICAL EXTRA

      EXTERNAL AFUN

      DATA MMAX/10/

      N=NN

      M=MIN(MM,MMAX,N-1)

      MPLUS=M+1

      IX=0

      IY=N+1

      IF (AFUN(1).GT.AFUN(N)) GOTO 94

   91 MID=(IX+IY)/2

      IF (X.GE.AFUN(MID)) GOTO 92

      IY=MID

      GOTO 93

   92 IX=MID

   93 IF (IY-IX.GT.1) GOTO 91

      GOTO 97

   94 MID=(IX+IY)/2

      IF (X.LE.AFUN(MID)) GOTO 95

      IY=MID

      GOTO 96

   95 IX=MID

   96 IF (IY-IX.GT.1) GOTO 94

   97 NPTS=M+2-MOD(M,2)

      IP=0

      L=0

      GOTO 99

   98 L=-L

      IF (L.GE.0) L=L+1

   99 ISUB=IX+L

      IF ((1.LE.ISUB).AND.(ISUB.LE.N)) GOTO 100

      NPTS=MPLUS

      GOTO 101

  100 IP=IP+1

      T(IP)=AFUN(ISUB)

      D(IP)=F(ISUB)

  101 IF (IP.LT.NPTS) GOTO 98

      EXTRA=NPTS.NE.MPLUS

      DO 14 L=1,M

      IF (.NOT.EXTRA) GOTO 12

      ISUB=MPLUS-L

      D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB))

   12 I=MPLUS

      DO 13 J=L,M

      ISUB=I-L

      D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB))

      I=I-1

   13 CONTINUE

   14 CONTINUE

      SUM=D(MPLUS)

      IF (EXTRA) SUM=0.5*(SUM+D(M+2))

      J=M

      DO 15 L=1,M

      SUM=D(J)+(X-T(J))*SUM

      J=J-1

   15 CONTINUE

      HWSTAB=SUM

      END
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