This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / v / funlux.F
1       SUBROUTINE FUNLXP (FUNC,XFCUM,X2LOW,X2HIGH)
2 C         F. JAMES,   Sept, 1994
3 C
4 C         Prepares the user function FUNC for FUNLUX
5 C         Inspired by and mostly copied from FUNPRE and FUNRAN
6 C         except that 
7 C    1. FUNLUX uses RANLUX underneath, 
8 C    2. FUNLXP expands the first and last bins to cater for
9 C              functions with long tails on left and/or right,
10 C    3. FUNLXP calls FUNPCT to do the actual finding of percentiles.
11 C    4. both FUNLXP and FUNPCT use RADAPT for Gaussian integration.
12 C
13       EXTERNAL FUNC
14       COMMON/FUNINT/TFTOT
15       DIMENSION XFCUM(200)
16       PARAMETER (RTEPS=0.0002) 
17       SAVE IFUNC
18       DATA IFUNC/0/
19       IFUNC = IFUNC + 1
20 C         FIND RANGE WHERE FUNCTION IS NON-ZERO.
21       CALL FUNLZ(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
22       XRANGE = XHIGH-XLOW
23       IF(XRANGE .LE. 0.)  THEN
24         WRITE (6,'(A,2G15.5)') ' FUNLXP finds function range .LE.0',
25      +  XLOW,XHIGH
26         GO TO 900
27       ENDIF
28       CALL RADAPT(FUNC,XLOW,XHIGH,1,RTEPS,0.,TFTOT ,UNCERT)
29       WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
30  1003 FORMAT(' FUNLXP: integral of USER FUNCTION',
31      +  I3,' from ',E12.5,' to ',E12.5,' is ',E14.6)
32 C
33 C      WRITE (6,'(A,A)') ' FUNLXP preparing ',
34 C     + 'first the whole range, then left tail, then right tail.'
35       CALL FUNPCT(FUNC,IFUNC,XLOW,XHIGH,XFCUM,1,99,TFTOT,IERR)
36       IF (IERR .GT. 0)  GO TO 900
37       X2 = XFCUM(3)
38       CALL RADAPT(FUNC,XLOW,X2,1,RTEPS,0.,TFTOT1 ,UNCERT)
39       CALL FUNPCT(FUNC,IFUNC,XLOW,X2 ,XFCUM,101,49,TFTOT1,IERR)
40       IF (IERR .GT. 0)  GO TO 900
41       X3 = XFCUM(98)
42       CALL RADAPT(FUNC,X3,XHIGH,1,RTEPS,0.,TFTOT2 ,UNCERT)
43       CALL FUNPCT(FUNC,IFUNC,X3,XHIGH,XFCUM,151,49,TFTOT2,IERR)
44       IF (IERR .GT. 0)  GO TO 900
45       WRITE(6,1001) IFUNC,XLOW,XHIGH
46  1001 FORMAT(' FUNLXP has prepared USER FUNCTION',I3,
47      + ' between',G12.3,' and',G12.3,' for FUNLUX')
48       RETURN
49   900 CONTINUE
50       WRITE (6,'(A)') ' Fatal error in FUNLXP. FUNLUX will not work.'
51       END
52 C
53       SUBROUTINE FUNPCT(FUNC,IFUNC,XLOW,XHIGH,XFCUM,NLO,NBINS,TFTOT,
54      +                  IERR)
55 C        Array XFCUM is filled from NLO to NLO+NBINS, which makes
56 C        the number of values NBINS+1, or the number of bins NBINS
57       EXTERNAL FUNC
58       DIMENSION XFCUM(*) 
59       PARAMETER (RTEPS=0.005, NZ=10, MAXZ=20, NITMAX=6,PRECIS=1.E-6)
60 CC      DOUBLE PRECISION TPCTIL, TZ, TCUM, XINCR, DTABS,
61 CC     +  TINCR, TZMAX, XBEST, DTBEST, DTPAR2
62 C
63       IERR = 0
64       IF (TFTOT .LE. 0.) GO TO 900
65       TPCTIL = TFTOT/NBINS
66       TZ = TPCTIL/NZ
67       TZMAX = TZ * 2.
68       XFCUM(NLO) = XLOW
69       XFCUM(NLO+NBINS) = XHIGH
70       X = XLOW
71       F = FUNC(X)
72       IF (F .LT. 0.) GO TO 900
73 C         Loop over percentile bins
74       DO 600 IBIN = NLO, NLO+NBINS-2
75       TCUM = 0.
76       X1 = X
77       F1 = F
78       DXMAX = (XHIGH -X) / NZ
79       FMIN = TZ/DXMAX
80       FMINZ = FMIN
81 C         Loop over trapezoids within a supposed percentile
82       DO 500 IZ= 1, MAXZ
83       XINCR = TZ/MAX(F1,FMIN,FMINZ)
84   350 X = X1 + XINCR
85       F = FUNC(X)
86       IF (F .LT. 0.) GO TO 900
87       TINCR = (X-X1) * 0.5 * (F+F1)
88       IF (TINCR .LT. TZMAX) GO TO 370
89       XINCR = XINCR * 0.5
90       GO TO 350
91   370 CONTINUE
92       TCUM = TCUM + TINCR
93       IF (TCUM .GE. TPCTIL*0.99) GO TO 520
94       FMINZ = TZ*F/ (TPCTIL-TCUM)
95       F1 = F
96       X1 = X
97   500 CONTINUE
98       WRITE(6,'(A)') ' FUNLUX:  WARNING. FUNPCT fails trapezoid.'
99 C         END OF TRAPEZOID LOOP
100 C         Adjust interval using Gaussian integration with 
101 C             Newton corrections since F is the derivative 
102   520 CONTINUE
103       X1 = XFCUM(IBIN)
104       XBEST = X
105       DTBEST = TPCTIL
106       TPART = TPCTIL
107 C         Allow for maximum NITMAX more iterations on RADAPT
108       DO 550 IHOME= 1, NITMAX
109   535 XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
110       X = XBEST + XINCR
111       X2 = X
112         IF (IHOME .GT. 1 .AND. X2 .EQ. XBEST) THEN
113         WRITE (6,'(A,G12.3)')
114      +  ' FUNLUX: WARNING from FUNPCT: insufficient precision at X=',X
115         GO TO 580
116         ENDIF
117       REFX = ABS(X)+PRECIS
118       CALL RADAPT(FUNC,X1,X2,1,RTEPS,0.,TPART2,UNCERT)
119       DTPAR2 = TPART2-TPCTIL
120       DTABS = ABS(DTPAR2)
121       IF(ABS(XINCR)/REFX .LT. PRECIS) GOTO 545
122       IF(DTABS .LT. DTBEST) GOTO 545
123       XINCR = XINCR * 0.5
124       GOTO 535
125   545 DTBEST = DTABS
126       XBEST = X
127       TPART = TPART2
128       F = FUNC(X)
129       IF(F .LT. 0.) GOTO 900
130       IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
131   550 CONTINUE
132       WRITE (6,'(A,I4)') 
133      +   ' FUNLUX: WARNING from FUNPCT: cannot converge, bin',IBIN
134 C
135   580 CONTINUE
136       XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
137       X = XBEST + XINCR
138       XFCUM(IBIN+1) = X
139       F = FUNC(X)
140       IF(F .LT. 0.) GOTO 900
141   600 CONTINUE
142 C         END OF LOOP OVER BINS
143       X1 = XFCUM(NLO+NBINS-1)
144       X2 = XHIGH
145       CALL RADAPT(FUNC,X1,X2,1,RTEPS,0.,TPART ,UNCERT)
146       ABERR = ABS(TPART-TPCTIL)/TFTOT
147 C      WRITE(6,1001) IFUNC,XLOW,XHIGH
148       IF(ABERR .GT. RTEPS)  WRITE(6,1002) ABERR
149       RETURN
150   900 WRITE(6,1000) X,F
151       IERR = 1
152       RETURN
153  1000 FORMAT(/' FUNLUX fatal error in FUNPCT: function negative:'/
154      + ,' at X=',E15.6,', F=',E15.6/)
155 C 1001 FORMAT(' FUNPCT has prepared USER FUNCTION',I3,
156 C     + ' between',G12.3,' and',G12.3,' for FUNLUX.')
157  1002 FORMAT(' WARNING: Relative error in cumulative distribution',
158      + ' may be as big as',F10.7)
159       END
160       SUBROUTINE FUNLUX(ARRAY,XRAN,LEN)
161 C         Generation of LEN random numbers in any given distribution,
162 C         by 4-point interpolation in the inverse cumulative distr.
163 C         which was previously generated by FUNLXP
164       COMMON/FUNINT/XUNI
165       DIMENSION ARRAY(200)
166       DIMENSION XRAN(LEN)
167 C  Bin width for main sequence, and its inverse
168       PARAMETER (GAP= 1./99.,  GAPINV=99.)
169 C  Top of left tail, bottom of right tail (each tail replaces 2 bins)
170       PARAMETER (TLEFT= 2./99.,BRIGHT=97./99.)
171 C  Bin width for minor sequences (tails), and its inverse
172       PARAMETER (GAPS=TLEFT/49.,  GAPINS=1./GAPS)
173 C
174 C   The array ARRAY is assumed to have the following structure: 
175 C        ARRAY(1-100) contains the 99 bins of the inverse cumulative
176 C                     distribution of the entire function.
177 C        ARRAY(101-150) contains the 49-bin blowup of main bins 
178 C                       1 and 2 (left tail of distribution)
179 C        ARRAY(151-200) contains the 49-bin blowup of main bins 
180 C                       98 and 99 (right tail of distribution)
181 C
182       CALL RANLUX(XRAN,LEN)
183       
184       DO 500 IBUF= 1, LEN
185       X = XRAN(IBUF)
186       J = INT(  X    *GAPINV) + 1
187       IF (J .LT. 3)  THEN
188          J1 = INT( X *GAPINS)
189              J = J1 + 101
190              J = MAX(J,102)
191              J = MIN(J,148)
192          P = (   X -GAPS*(J1-1)) * GAPINS
193          A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
194          B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
195          XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
196       ELSE IF (J .GT. 97)  THEN
197          J1 = INT((X-BRIGHT)*GAPINS) 
198              J = J1 + 151
199              J = MAX(J,152)
200              J = MIN(J,198)
201          P = (X -BRIGHT -GAPS*(J1-1)) * GAPINS
202          A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
203          B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
204          XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
205       ELSE
206 CC      J = MAX(J,2)
207 CC      J = MIN(J,98)
208          P = (   X -GAP*(J-1)) * GAPINV
209          A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
210          B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
211          XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
212       ENDIF
213   500 CONTINUE
214       XUNI = X
215       RETURN
216       END
217       SUBROUTINE FUNLZ(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
218 C         FIND RANGE WHERE FUNC IS NON-ZERO.
219 C         WRITTEN 1980, F. JAMES
220 C         MODIFIED, NOV. 1985, TO FIX BUG AND GENERALIZE
221 C         TO FIND SIMPLY-CONNECTED NON-ZERO REGION (XLOW,XHIGH)
222 C         ANYWHERE WITHIN THE GIVEN REGION (X2LOW,H2HIGH).
223 C            WHERE 'ANYWHERE' MEANS EITHER AT THE LOWER OR UPPER
224 C            EDGE OF THE GIVEN REGION, OR, IF IN THE MIDDLE,
225 C            COVERING AT LEAST 1% OF THE GIVEN REGION.
226 C         OTHERWISE IT IS NOT GUARANTEED TO FIND THE NON-ZERO REGION.
227 C         IF FUNCTION EVERYWHERE ZERO, FUNLZ SETS XLOW=XHIGH=0.
228       EXTERNAL FUNC
229       XLOW = X2LOW
230       XHIGH = X2HIGH
231 C         FIND OUT IF FUNCTION IS ZERO AT ONE END OR BOTH
232       XMID = XLOW
233       IF (FUNC(XLOW) .GT. 0.) GO TO 120
234       XMID = XHIGH
235       IF (FUNC(XHIGH) .GT. 0.)  GO TO 50
236 C         FUNCTION IS ZERO AT BOTH ENDS,
237 C         LOOK FOR PLACE WHERE IT IS NON-ZERO.
238       DO 30 LOGN= 1, 7
239       NSLICE = 2**LOGN
240       DO 20 I= 1, NSLICE, 2
241       XMID = XLOW + I * (XHIGH-XLOW) / NSLICE
242       IF (FUNC(XMID) .GT. 0.)  GO TO 50
243    20 CONTINUE
244    30 CONTINUE
245 C         FALLING THROUGH LOOP MEANS CANNOT FIND NON-ZERO VALUE
246       WRITE(6,554)
247       WRITE(6,555) XLOW, XHIGH
248       XLOW = 0.
249       XHIGH = 0.
250       GO TO 220
251 C
252    50 CONTINUE
253 C         DELETE 'LEADING' ZERO RANGE
254       XH = XMID
255       XL = XLOW
256       DO 70 K= 1, 20
257       XNEW = 0.5*(XH+XL)
258       IF (FUNC(XNEW) .EQ. 0.) GO TO 68
259       XH = XNEW
260       GO TO 70
261    68 XL = XNEW
262    70 CONTINUE
263       XLOW = XL
264       WRITE(6,555) X2LOW,XLOW
265   120 CONTINUE
266       IF (FUNC(XHIGH) .GT. 0.) GO TO 220
267 C         DELETE 'TRAILING' RANGE OF ZEROES
268       XL = XMID
269       XH = XHIGH
270       DO 170 K= 1, 20
271       XNEW = 0.5*(XH+XL)
272       IF (FUNC(XNEW) .EQ. 0.) GO TO 168
273       XL = XNEW
274       GO TO 170
275   168 XH = XNEW
276   170 CONTINUE
277       XHIGH = XH
278       WRITE(6,555) XHIGH, X2HIGH
279 C
280   220 CONTINUE
281       RETURN
282   554 FORMAT('0CANNOT FIND NON-ZERO FUNCTION VALUE')
283   555 FORMAT(' FUNCTION IS ZERO FROM X=',E12.5,' TO ',E12.5)
284       END