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


CDECK  ID>, HWUSTA.

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

*-- Author :    Bryan Webber

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

      SUBROUTINE HWUSTA(NAME)

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

C     MAKES PARTICLE TYPE 'NAME' STABLE

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

      INCLUDE 'HERWIG61.INC'

      INTEGER IPDG,IWIG

      CHARACTER*8 NAME

      CALL HWUIDT(3,IPDG,IWIG,NAME)

      IF (IWIG.EQ.20) CALL HWWARN('HWUSTA',500,*999)

      RSTAB(IWIG)=.TRUE.

      WRITE (6,10) IWIG,NAME

   10 FORMAT(/10X,'PARTICLE TYPE',I4,'=',A8,' SET STABLE')

  999 END

CDECK  ID>, HWUTAB.

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

*-- Author :    Adapted by Bryan Webber

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

      FUNCTION HWUTAB(F,A,NN,X,MM)

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

C     MODIFIED CERN INTERPOLATION ROUTINE DIVDIF

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

      IMPLICIT NONE

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

      DOUBLE PRECISION HWUTAB,SUM,X,F(NN),A(NN),T(20),D(20)

      LOGICAL EXTRA

      DATA MMAX/10/

      N=NN

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

      MPLUS=M+1

      IX=0

      IY=N+1

      IF (A(1).GT.A(N)) GOTO 4

    1 MID=(IX+IY)/2

      IF (X.GE.A(MID)) GOTO 2

      IY=MID

      GOTO 3

    2 IX=MID

    3 IF (IY-IX.GT.1) GOTO 1

      GOTO 7

    4 MID=(IX+IY)/2

      IF (X.LE.A(MID)) GOTO 5

      IY=MID

      GOTO 6

    5 IX=MID

    6 IF (IY-IX.GT.1) GOTO 4

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

      IP=0

      L=0

      GOTO 9

    8 L=-L

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

    9 ISUB=IX+L

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

      NPTS=MPLUS

      GOTO 11

   10 IP=IP+1

      T(IP)=A(ISUB)

      D(IP)=F(ISUB)

   11 IF (IP.LT.NPTS) GOTO 8

      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

      HWUTAB=SUM

      END
Commit	Line	Data
3820ca8e	1
	2	CDECK ID>, HWUSTA.
	3
	4	*CMZ :- -26/04/91 10.18.58 by Bryan Webber
	5
	6	*-- Author : Bryan Webber
	7
	8	C-----------------------------------------------------------------------
	9
	10	SUBROUTINE HWUSTA(NAME)
	11
	12	C-----------------------------------------------------------------------
	13
	14	C MAKES PARTICLE TYPE 'NAME' STABLE
	15
	16	C-----------------------------------------------------------------------
	17
	18	INCLUDE 'HERWIG61.INC'
	19
	20	INTEGER IPDG,IWIG
	21
	22	CHARACTER*8 NAME
	23
	24	CALL HWUIDT(3,IPDG,IWIG,NAME)
	25
	26	IF (IWIG.EQ.20) CALL HWWARN('HWUSTA',500,*999)
	27
	28	RSTAB(IWIG)=.TRUE.
	29
	30	WRITE (6,10) IWIG,NAME
	31
	32	10 FORMAT(/10X,'PARTICLE TYPE',I4,'=',A8,' SET STABLE')
	33
	34	999 END
	35
	36	CDECK ID>, HWUTAB.
	37
	38	*CMZ :- -26/04/91 11.11.56 by Bryan Webber
	39
	40	*-- Author : Adapted by Bryan Webber
	41
	42	C-----------------------------------------------------------------------
	43
	44	FUNCTION HWUTAB(F,A,NN,X,MM)
	45
	46	C-----------------------------------------------------------------------
	47
	48	C MODIFIED CERN INTERPOLATION ROUTINE DIVDIF
	49
	50	C-----------------------------------------------------------------------
	51
	52	IMPLICIT NONE
	53
	54	INTEGER NN,MM,MMAX,N,M,MPLUS,IX,IY,MID,NPTS,IP,I,J,L,ISUB
	55
	56	DOUBLE PRECISION HWUTAB,SUM,X,F(NN),A(NN),T(20),D(20)
	57
	58	LOGICAL EXTRA
	59
	60	DATA MMAX/10/
	61
	62	N=NN
	63
	64	M=MIN(MM,MMAX,N-1)
65
66	MPLUS=M+1
67
68	IX=0
69
70	IY=N+1
71
72	IF (A(1).GT.A(N)) GOTO 4
73
74	1 MID=(IX+IY)/2
75
76	IF (X.GE.A(MID)) GOTO 2
77
78	IY=MID
79
80	GOTO 3
81
82	2 IX=MID
83
84	3 IF (IY-IX.GT.1) GOTO 1
85
86	GOTO 7
87
88	4 MID=(IX+IY)/2
89
90	IF (X.LE.A(MID)) GOTO 5
91
92	IY=MID
93
94	GOTO 6
95
96	5 IX=MID
97
98	6 IF (IY-IX.GT.1) GOTO 4
99
100	7 NPTS=M+2-MOD(M,2)
101
102	IP=0
103
104	L=0
105
106	GOTO 9
107
108	8 L=-L
109
110	IF (L.GE.0) L=L+1
111
112	9 ISUB=IX+L
113
114	IF ((1.LE.ISUB).AND.(ISUB.LE.N)) GOTO 10
115
116	NPTS=MPLUS
117
118	GOTO 11
119
120	10 IP=IP+1
121
122	T(IP)=A(ISUB)
123
124	D(IP)=F(ISUB)
125
126	11 IF (IP.LT.NPTS) GOTO 8
127
128	EXTRA=NPTS.NE.MPLUS
129
130	DO 14 L=1,M
131
132	IF (.NOT.EXTRA) GOTO 12
133
134	ISUB=MPLUS-L
135
136	D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB))
137
138	12 I=MPLUS
139
140	DO 13 J=L,M
141
142	ISUB=I-L
143
144	D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB))
145
146	I=I-1
147
148	13 CONTINUE
149
150	14 CONTINUE
151
152	SUM=D(MPLUS)
153
154	IF (EXTRA) SUM=0.5*(SUM+D(M+2))
155
156	J=M
157
158	DO 15 L=1,M
159
160	SUM=D(J)+(X-T(J))*SUM
161
162	J=J-1
163
164	15 CONTINUE
165
166	HWUTAB=SUM
167
168	END