HERWIG/src/hwusta.f

   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