[u/mrichter/AliRoot.git] / PHOS / shaker / funpre.f

      SUBROUTINE FUNPRE (FUNC,XFCUM,X2LOW,X2HIGH)
C         F. JAMES,    MAY, 1976
C         MODIFIED OCT, 1980 TO ADD PRINTOUT OF INTEGRAL
C         MODIFIED DEC., 1980 TO DELETE LEADING AND TRAILING
C            RANGES OF X WHERE FUNCTION IS ZERO.
C         MODIFIED JUNE,1982 TO FIX POSSIBLE INFINITE LOOP.
C
C         PREPARES THE USER FUNCTION "FUNC" FOR FUNRAN
C         BY FINDING THE PERCENTILES
C         (IN EFFECT, INVERTING THE CUMULATIVE DISTRIBUTION)
      EXTERNAL FUNC
      COMMON/FUNINT/TFTOT
      DIMENSION XFCUM(100)
      DATA NBINS / 99/
      DATA NZ / 10/
      DATA MAXZ / 20/
      DATA NITMAX / 6 /
      DATA IFUNC/0/
      IFUNC = IFUNC + 1
C         FIND MACHINE ACCURACY
      COMP1 = 1.0
      DO 200 I= 1, 100
      COMP1 = COMP1*0.5
      COMP2 = 1.0 - COMP1
      IF(COMP2 .EQ. 1.0) GOTO 210
  200 CONTINUE
      COMP1 = 1.0E-10
  210 PRECIS = COMP1
C         FIND RANGE WHERE FUNCTION IS NON-ZERO.
      CALL FUNZER(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
      XRANGE = XHIGH-XLOW
      IF(XRANGE .LE. 0.) GOTO 900
      RTEPS = MAX(0.0001,PRECIS*10.)
      TFTOT = GAUSS(FUNC,XLOW,XHIGH,RTEPS)
C         PRINT OUT VALUE OF NORMALIZATION INTEGRAL
      WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
      RTEPS = 0.001
      IF(TFTOT .LE. 0.) GOTO 900
      TPCTIL = TFTOT/NBINS
      TZ = TPCTIL/NZ
      TZMAX = TZ * 2.
       XFCUM(1) = XLOW
      XFCUM(NBINS+1) = XHIGH
      X = XLOW
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
      NBINM1 = NBINS - 1
C         LOOP OVER BINS (HUNDREDTH PERCENTILES)
      DO 600 IBIN = 1, NBINM1
      TCUM = 0.
      X1 = X
      F1 = F
      DXMAX = (XHIGH -X) / NZ
      FMIN = TZ/DXMAX
      FMINZ = FMIN
C         LOOP OVER TRAPEZOIDS WITHIN A SUPPOSED PERCENTILE
      DO 500 IZ= 1, MAXZ
      XINCR = TZ/MAX(F1,FMIN,FMINZ)
  350 X = X1 + XINCR
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
      TINCR = (X-X1) * 0.5 * (F+F1)
      IF(TINCR .LT. TZMAX) GOTO 370
      XINCR = XINCR * 0.5
      GOTO 350
  370 CONTINUE
      TCUM = TCUM + TINCR
      IF(TCUM .GE. TPCTIL*0.99) GOTO 520
      FMINZ = TZ*F/ (TPCTIL-TCUM)
      F1 = F
      X1 = X
  500 CONTINUE
      WRITE(6,2000)
 2000 FORMAT('0      FAILURE TO FIND TRAPEZOID   HELP')
C         END OF TRAPEZOID LOOP
C         ADJUST TRAPEZOID INTEGRAL BY GAUSS WITH NEWTON CORRECTION
  520 CONTINUE
      X1 = XFCUM(IBIN)
      XBEST = X
      DTBEST = TPCTIL
      TPART = TPCTIL
C         ALLOW FOR MAXIMUM NITMAX MORE ITERATIONS ON GAUSS
      DO 550 IHOME= 1, NITMAX
      XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
  535 X = XBEST + XINCR
      X2 = X
      TPART2 = GAUSS(FUNC,X1,X2,RTEPS)
      DTPAR2 = TPART2-TPCTIL
      DTABS = ABS(DTPAR2)
      IF(ABS(XINCR) .LT. PRECIS) GOTO 545
      IF(DTABS .LT. DTBEST) GOTO 545
      XINCR = XINCR * 0.5
      GOTO 535
  545 DTBEST = DTABS
      XBEST = X
      IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
      TPART = TPART2
      F = FUNC(X)
  550 CONTINUE
      IHOME = NITMAX
C
  580 CONTINUE
      XFCUM(IBIN+1) = X
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
  600 CONTINUE
C         END OF LOOP OVER BINS
      X1 = XFCUM(NBINS)
      X2 = XHIGH
      TPART = GAUSS(FUNC,X1,X2,RTEPS)
      ABERR = ABS(TPART-TPCTIL)/TFTOT
      WRITE(6,1001) IFUNC,ABERR
      IF(ABERR .GT. RTEPS)  WRITE(6,1002)
      RETURN
  900 WRITE(6,1000) X,F,XLOW,XHIGH
      RETURN
 1000 FORMAT('0FUNPRE FINDS NEGATIVE FUNCTION VALUE OR RANGE OF X'/
     + ,' X=',E15.6,', F=',E15.6,20X,'XLOW=',E15.6,' XHIGH=',E15.6/)
 1001 FORMAT(' SUBROUTINE FUNPRE HAS PREPARED USER FUNCTION NUMBER',I4,
     + ' FOR FUNRAN'/' MAXIMUM RELATIVE ERROR IN CUMULATIVE',
     + 'DISTRIBUTION WILL BE',E15.5)
 1002 FORMAT('+',80X,'WARNING,THIS MAY BE TOO BIG'//)
 1003 FORMAT('0SUBROUTINE FUNPRE FINDS THE INTEGRAL OF USER FUNCTION',
     +  I2,' FROM ',E12.5,' TO ',E12.5,' IS ',E14.6)
      END
Commit	Line	Data
fe4da5cc	1	SUBROUTINE FUNPRE (FUNC,XFCUM,X2LOW,X2HIGH)
	2	C F. JAMES, MAY, 1976
	3	C MODIFIED OCT, 1980 TO ADD PRINTOUT OF INTEGRAL
	4	C MODIFIED DEC., 1980 TO DELETE LEADING AND TRAILING
	5	C RANGES OF X WHERE FUNCTION IS ZERO.
	6	C MODIFIED JUNE,1982 TO FIX POSSIBLE INFINITE LOOP.
	7	C
	8	C PREPARES THE USER FUNCTION "FUNC" FOR FUNRAN
	9	C BY FINDING THE PERCENTILES
	10	C (IN EFFECT, INVERTING THE CUMULATIVE DISTRIBUTION)
	11	EXTERNAL FUNC
	12	COMMON/FUNINT/TFTOT
	13	DIMENSION XFCUM(100)
	14	DATA NBINS / 99/
	15	DATA NZ / 10/
	16	DATA MAXZ / 20/
	17	DATA NITMAX / 6 /
	18	DATA IFUNC/0/
	19	IFUNC = IFUNC + 1
	20	C FIND MACHINE ACCURACY
	21	COMP1 = 1.0
	22	DO 200 I= 1, 100
	23	COMP1 = COMP1*0.5
	24	COMP2 = 1.0 - COMP1
	25	IF(COMP2 .EQ. 1.0) GOTO 210
	26	200 CONTINUE
	27	COMP1 = 1.0E-10
	28	210 PRECIS = COMP1
	29	C FIND RANGE WHERE FUNCTION IS NON-ZERO.
	30	CALL FUNZER(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
	31	XRANGE = XHIGH-XLOW
	32	IF(XRANGE .LE. 0.) GOTO 900
	33	RTEPS = MAX(0.0001,PRECIS*10.)
	34	TFTOT = GAUSS(FUNC,XLOW,XHIGH,RTEPS)
	35	C PRINT OUT VALUE OF NORMALIZATION INTEGRAL
	36	WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
	37	RTEPS = 0.001
	38	IF(TFTOT .LE. 0.) GOTO 900
	39	TPCTIL = TFTOT/NBINS
	40	TZ = TPCTIL/NZ
	41	TZMAX = TZ * 2.
	42	XFCUM(1) = XLOW
	43	XFCUM(NBINS+1) = XHIGH
	44	X = XLOW
	45	F = FUNC(X)
	46	IF(F .LT. 0.) GOTO 900
	47	NBINM1 = NBINS - 1
	48	C LOOP OVER BINS (HUNDREDTH PERCENTILES)
	49	DO 600 IBIN = 1, NBINM1
	50	TCUM = 0.
	51	X1 = X
	52	F1 = F
	53	DXMAX = (XHIGH -X) / NZ
	54	FMIN = TZ/DXMAX
	55	FMINZ = FMIN
	56	C LOOP OVER TRAPEZOIDS WITHIN A SUPPOSED PERCENTILE
	57	DO 500 IZ= 1, MAXZ
	58	XINCR = TZ/MAX(F1,FMIN,FMINZ)
	59	350 X = X1 + XINCR
	60	F = FUNC(X)
	61	IF(F .LT. 0.) GOTO 900
	62	TINCR = (X-X1) * 0.5 * (F+F1)
	63	IF(TINCR .LT. TZMAX) GOTO 370
	64	XINCR = XINCR * 0.5
65	GOTO 350
66	370 CONTINUE
67	TCUM = TCUM + TINCR
68	IF(TCUM .GE. TPCTIL*0.99) GOTO 520
69	FMINZ = TZ*F/ (TPCTIL-TCUM)
70	F1 = F
71	X1 = X
72	500 CONTINUE
73	WRITE(6,2000)
74	2000 FORMAT('0 FAILURE TO FIND TRAPEZOID HELP')
75	C END OF TRAPEZOID LOOP
76	C ADJUST TRAPEZOID INTEGRAL BY GAUSS WITH NEWTON CORRECTION
77	520 CONTINUE
78	X1 = XFCUM(IBIN)
79	XBEST = X
80	DTBEST = TPCTIL
81	TPART = TPCTIL
82	C ALLOW FOR MAXIMUM NITMAX MORE ITERATIONS ON GAUSS
83	DO 550 IHOME= 1, NITMAX
84	XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
85	535 X = XBEST + XINCR
86	X2 = X
87	TPART2 = GAUSS(FUNC,X1,X2,RTEPS)
88	DTPAR2 = TPART2-TPCTIL
89	DTABS = ABS(DTPAR2)
90	IF(ABS(XINCR) .LT. PRECIS) GOTO 545
91	IF(DTABS .LT. DTBEST) GOTO 545
92	XINCR = XINCR * 0.5
93	GOTO 535
94	545 DTBEST = DTABS
95	XBEST = X
96	IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
97	TPART = TPART2
98	F = FUNC(X)
99	550 CONTINUE
100	IHOME = NITMAX
101	C
102	580 CONTINUE
103	XFCUM(IBIN+1) = X
104	F = FUNC(X)
105	IF(F .LT. 0.) GOTO 900
106	600 CONTINUE
107	C END OF LOOP OVER BINS
108	X1 = XFCUM(NBINS)
109	X2 = XHIGH
110	TPART = GAUSS(FUNC,X1,X2,RTEPS)
111	ABERR = ABS(TPART-TPCTIL)/TFTOT
112	WRITE(6,1001) IFUNC,ABERR
113	IF(ABERR .GT. RTEPS) WRITE(6,1002)
114	RETURN
115	900 WRITE(6,1000) X,F,XLOW,XHIGH
116	RETURN
117	1000 FORMAT('0FUNPRE FINDS NEGATIVE FUNCTION VALUE OR RANGE OF X'/
118	+ ,' X=',E15.6,', F=',E15.6,20X,'XLOW=',E15.6,' XHIGH=',E15.6/)
119	1001 FORMAT(' SUBROUTINE FUNPRE HAS PREPARED USER FUNCTION NUMBER',I4,
120	+ ' FOR FUNRAN'/' MAXIMUM RELATIVE ERROR IN CUMULATIVE',
121	+ 'DISTRIBUTION WILL BE',E15.5)
122	1002 FORMAT('+',80X,'WARNING,THIS MAY BE TOO BIG'//)
123	1003 FORMAT('0SUBROUTINE FUNPRE FINDS THE INTEGRAL OF USER FUNCTION',
124	+ I2,' FROM ',E12.5,' TO ',E12.5,' IS ',E14.6)
125	END