*
* $Id: dfnpre.F,v 1.1.1.1 1996/03/08 16:58:53 mclareni Exp $
*
* $Log: dfnpre.F,v $
* Revision 1.1.1.1  1996/03/08 16:58:53  mclareni
* Eurodec
*
*
#include "eurodec/pilot.h"
      SUBROUTINE DFNPRE(FUNC,XFCUM,XLOW,XHIGH)
C.----------------------------------------------------------------------
C.
C.    PREPARES THE USER FUNCTION "FUNC" FOR FUNRAN BY FINDING THE
C.    PERCENTILES (IN EFFECT, INVERTING THE CUMULATIVE DISTRIBUTION).
C.    BY F. JAMES, MAY, 1976 AND H. NEWMAN, JUNE 1978.
C.    LAST UPDATE: 11/11/88
C.
C.----------------------------------------------------------------------
#if defined(CERNLIB_DOUBLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#endif
#include "eurodec/inpout.inc"
      DIMENSION XFCUM(100),DTPAR2(10)
      SAVE IFUNC
      EXTERNAL FUNC
      DATA NBINS,NZ,MAXZ,NITMAX,IFUNC/ 99, 10, 20, 6, 0/
      DATA IPRINT/ 0/
C--
C--   INITIALIZE
      IFUNC=IFUNC+1
      XRANGE=XHIGH-XLOW
C--
C--   CHECK RANGE
      IF (XRANGE.LE.0.0) CALL ERRORD(32,' ',REAL(XRANGE))
C--
C--   CHECK ON NEGATIVE FUNCTION VALUES
      RTEPS=1.E-4
      TFTOT=DGAUSS(FUNC,XLOW,XHIGH,RTEPS)
      RTEPS=1.E-3
      IF (TFTOT.LE.0.0) CALL ERRORD(33,' ',REAL(TFTOT))
      TPCTIL=TFTOT/FLOAT(NBINS)
      TZ=TPCTIL/FLOAT(NZ)
      TZMAX=TZ*2.
      XFCUM(1)=XLOW
      XFCUM(NBINS+1)=XHIGH
      X=XLOW
      F=FUNC(X)
      IF (F.LT.0.0) CALL ERRORD(33,' ',REAL(F))
      NBINM1=NBINS-1
C--
C--   LOOP OVER BINS (HUNDREDTH PERCENTILES)
      DO 90 IBIN=1,NBINM1
        TCUM=0.0
        X1=X
        F1=F
        DXMAX=(XHIGH-X)/FLOAT(NZ)
        FMIN=TZ/DXMAX
        FMINZ=FMIN
C--
C--   LOOP OVER TRAPEZOIDS WITHIN A SUPPOSED PERCENTILE
        DO 30 IZ=1,MAXZ
          XINCR=TZ/MAX(F1,FMIN,FMINZ)
   10     X=X1+XINCR
          F=FUNC(X)
          IF (F.LT.0.0) CALL ERRORD(33,' ',REAL(F))
          TINCR=(X-X1)*0.5*(F+F1)
          IF (TINCR.LT.TZMAX) GOTO 20
          XINCR=XINCR*0.5
          GOTO 10
   20     TCUM=TCUM+TINCR
          IF (TCUM.GE.TPCTIL*0.99) GOTO 40
          FMINZ=TZ*F/(TPCTIL-TCUM)
          F1=F
   30   X1=X
C--
C--   ERROR PRINT
        CALL ERRORD(31,' ',0.)
C--
C--   ADJUST TRAPEZOID INTEGRAL BY GAUSS WITH NEWTON CORRECTION
   40   X1=XFCUM(IBIN)
        XBEST=X
        DTBEST=TPCTIL
        DTCUM=TCUM-TPCTIL
        TPART=TPCTIL
C--
C--   ALLOW FOR MAXIMUM NITMAX MORE ITERATIONS ON GAUSS
        DO 70 IHOME=1,NITMAX
          XINCR=(TPCTIL-TPART)/MAX(F,FMIN)
   50     X=XBEST+XINCR
          X2=X
          TPART2=DGAUSS(FUNC,X1,X2,RTEPS)
          DTPAR2(IHOME)=TPART2-TPCTIL
          DTABS=ABS(DTPAR2(IHOME))
          IF (DTABS.LT.DTBEST) GOTO 60
          XINCR=XINCR*0.5
          GOTO 50
   60     DTBEST=DTABS
          XBEST=X
          IF (DTABS.LT.RTEPS*TPCTIL) GOTO 80
          TPART=TPART2
   70   F=FUNC(X)
        IHOME=NITMAX
   80   XFCUM(IBIN+1)=X
        F=FUNC(X)
        IF (F.LT.0.0) CALL ERRORD(33,' ',REAL(F))
   90 CONTINUE
C--
C--   END OF LOOP OVER BINS
      X1=XFCUM(NBINS)
      X2=XHIGH
      TPART=DGAUSS(FUNC,X1,X2,RTEPS)
      ABERR=ABS(TPART-TPCTIL)/TFTOT
      IF ((IPRINT.NE.0).OR.(ABERR.GT.RTEPS)) THEN
        WRITE(LUN2,9000)
        WRITE(LUN2,9010) IFUNC
        WRITE(LUN2,9020) ABERR
        IF (ABERR.GT.RTEPS) WRITE(LUN2,9030)
      ENDIF
      RETURN
 9000 FORMAT(1H0,80('*'))
 9010 FORMAT(1H0,6X,'SUBROUTINE FUNPRE HAS PREPARED USER FUNCTION',
     &' NUMBER: ',I2,' FOR FUNRAN.')
 9020 FORMAT(1H ,3X,'MAXIMUM RELATIVE ERROR IN CUMULATIVE DISTRIBUTION',
     &'WILL BE: ',F14.10)
 9030 FORMAT(1H ,17X,'****** WARNING: THIS MAY BE TOO BIG !! ******')
      END