HERWIG/src/hwstab.f

   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