MINICERN/mathlib/gen/v/hisran.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:57  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE HISRAN(Y,N,XLO,XWID,XRAN)
  11 C         SUBROUTINE TO GENERATE RANDOM NUMBERS
  12 C         ACCORDING TO AN EMPIRICAL DISTRIBUTION
  13 C         SUPPLIED BY THE USER IN THE FORM OF A HISTOGRAM
  14 C         F. JAMES,    MAY, 1976
  15       DIMENSION Y(*)
  16 #if !defined(CERNLIB_F90)
  17       DATA IERR,NTRY,NXHRAN,NXHPRE/0,3HRAN,3HRAN,3HPRE/
  18 #endif
  19 #if defined(CERNLIB_F90)
  20       INTEGER :: IERR = 0, NTRY = transfer('RAN ', 0),                  &
  21            NXHRAN = transfer('RAN ', 0), NXHPRE = transfer('PRE ', 0)
  22 #endif
  23       IF(Y(N).EQ.1.0) GOTO 200
  24       WRITE(6,1001) Y(N)
  25  1001 FORMAT('0SUBROUTINE HISRAN FINDS Y(N) NOT EQUAL TO 1.0     Y(N)='
  26      +,E15.6/' ASSUMES USER HAS SUPPLIED HISTOGRAM RATHER THAN CUMUL',
  27      +'ATIVE DISTRIBUTION AND HAS FORGOTTEN TO CALL HISPRE'/)
  28       NTRY=NXHRAN
  29       GOTO 50
  30 C         INITIALIZE HISTOGRAM TO FORM CUMULATIVE DISTRIBUTION
  31 #if (defined(CERNLIB_CDC))&&(defined(CERNLIB_F4))
  32       ENTRY HISPRE
  33 #endif
  34 #if !defined(CERNLIB_CDC)||!defined(CERNLIB_F4)
  35       ENTRY HISPRE(Y,N)
  36 #endif
  37       NTRY=NXHPRE
  38    50 CONTINUE
  39       YTOT = 0.
  40       DO 100 I= 1, N
  41       IF(Y(I).LT.0.) GOTO 900
  42       YTOT = YTOT + Y(I)
  43   100 Y(I) = YTOT
  44       IF(YTOT.LE.0.) GOTO 900
  45       YINV = 1.0/YTOT
  46       DO 110 I= 1, N
  47   110 Y(I) = Y(I) * YINV
  48       Y(N) = 1.0
  49       IF(NTRY.EQ.NXHPRE) RETURN
  50 C         NOW GENERATE RANDOM NUMBER BETWEEN 0 AND ONE
  51   200 CONTINUE
  52       YR = RNDM(-1)
  53 C         AND TRANSFORM IT INTO THE CORRESPONDING X-VALUE
  54       L = LOCATF(Y,N,YR)
  55       IF(L.EQ.0) GOTO 240
  56       IF(L.GT.0) GOTO 250
  57 C         USUALLY COME HERE.
  58       L = ABS(L)
  59       XRAN = XLO + XWID * (L +((YR-Y(L))/(Y(L+1)-Y(L))))
  60       RETURN
  61 C         POINT FALLS IN FIRST BIN.  SPECIAL CASE
  62   240 XRAN = XLO + XWID * (YR/Y(1))
  63       RETURN
  64 C         GUARD AGAINST SPECIAL CASE OF FALLING ON EMPTY BIN
  65   250 XRAN = XLO + L * XWID
  66       RETURN
  67   900 CONTINUE
  68       IERR = IERR + 1
  69       IF(IERR.LT.6) WRITE(6,1000)NTRY
  70  1000 FORMAT('0ERROR IN INPUT DATA FOR HIS',A3,' VALUES NOT ALL >0'/)
  71       WRITE(6,1002) (Y(K),K=1,N)
  72  1002 FORMAT(1X,10F13.7)
  73       XRAN = 0.
  74       RETURN
  75       END