[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / v / funlux.F

      SUBROUTINE FUNLXP (FUNC,XFCUM,X2LOW,X2HIGH)
C         F. JAMES,   Sept, 1994
C
C         Prepares the user function FUNC for FUNLUX
C         Inspired by and mostly copied from FUNPRE and FUNRAN
C         except that 
C    1. FUNLUX uses RANLUX underneath, 
C    2. FUNLXP expands the first and last bins to cater for
C              functions with long tails on left and/or right,
C    3. FUNLXP calls FUNPCT to do the actual finding of percentiles.
C    4. both FUNLXP and FUNPCT use RADAPT for Gaussian integration.
C
      EXTERNAL FUNC
      COMMON/FUNINT/TFTOT
      DIMENSION XFCUM(200)
      PARAMETER (RTEPS=0.0002) 
      SAVE IFUNC
      DATA IFUNC/0/
      IFUNC = IFUNC + 1
C         FIND RANGE WHERE FUNCTION IS NON-ZERO.
      CALL FUNLZ(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
      XRANGE = XHIGH-XLOW
      IF(XRANGE .LE. 0.)  THEN
        WRITE (6,'(A,2G15.5)') ' FUNLXP finds function range .LE.0',
     +  XLOW,XHIGH
        GO TO 900
      ENDIF
      CALL RADAPT(FUNC,XLOW,XHIGH,1,RTEPS,0.,TFTOT ,UNCERT)
      WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
 1003 FORMAT(' FUNLXP: integral of USER FUNCTION',
     +  I3,' from ',E12.5,' to ',E12.5,' is ',E14.6)
C
C      WRITE (6,'(A,A)') ' FUNLXP preparing ',
C     + 'first the whole range, then left tail, then right tail.'
      CALL FUNPCT(FUNC,IFUNC,XLOW,XHIGH,XFCUM,1,99,TFTOT,IERR)
      IF (IERR .GT. 0)  GO TO 900
      X2 = XFCUM(3)
      CALL RADAPT(FUNC,XLOW,X2,1,RTEPS,0.,TFTOT1 ,UNCERT)
      CALL FUNPCT(FUNC,IFUNC,XLOW,X2 ,XFCUM,101,49,TFTOT1,IERR)
      IF (IERR .GT. 0)  GO TO 900
      X3 = XFCUM(98)
      CALL RADAPT(FUNC,X3,XHIGH,1,RTEPS,0.,TFTOT2 ,UNCERT)
      CALL FUNPCT(FUNC,IFUNC,X3,XHIGH,XFCUM,151,49,TFTOT2,IERR)
      IF (IERR .GT. 0)  GO TO 900
      WRITE(6,1001) IFUNC,XLOW,XHIGH
 1001 FORMAT(' FUNLXP has prepared USER FUNCTION',I3,
     + ' between',G12.3,' and',G12.3,' for FUNLUX')
      RETURN
  900 CONTINUE
      WRITE (6,'(A)') ' Fatal error in FUNLXP. FUNLUX will not work.'
      END
C
      SUBROUTINE FUNPCT(FUNC,IFUNC,XLOW,XHIGH,XFCUM,NLO,NBINS,TFTOT,
     +                  IERR)
C        Array XFCUM is filled from NLO to NLO+NBINS, which makes
C        the number of values NBINS+1, or the number of bins NBINS
      EXTERNAL FUNC
      DIMENSION XFCUM(*) 
      PARAMETER (RTEPS=0.005, NZ=10, MAXZ=20, NITMAX=6,PRECIS=1.E-6)
CC      DOUBLE PRECISION TPCTIL, TZ, TCUM, XINCR, DTABS,
CC     +  TINCR, TZMAX, XBEST, DTBEST, DTPAR2
C
      IERR = 0
      IF (TFTOT .LE. 0.) GO TO 900
      TPCTIL = TFTOT/NBINS
      TZ = TPCTIL/NZ
      TZMAX = TZ * 2.
      XFCUM(NLO) = XLOW
      XFCUM(NLO+NBINS) = XHIGH
      X = XLOW
      F = FUNC(X)
      IF (F .LT. 0.) GO TO 900
C         Loop over percentile bins
      DO 600 IBIN = NLO, NLO+NBINS-2
      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.) GO TO 900
      TINCR = (X-X1) * 0.5 * (F+F1)
      IF (TINCR .LT. TZMAX) GO TO 370
      XINCR = XINCR * 0.5
      GO TO 350
  370 CONTINUE
      TCUM = TCUM + TINCR
      IF (TCUM .GE. TPCTIL*0.99) GO TO 520
      FMINZ = TZ*F/ (TPCTIL-TCUM)
      F1 = F
      X1 = X
  500 CONTINUE
      WRITE(6,'(A)') ' FUNLUX:  WARNING. FUNPCT fails trapezoid.'
C         END OF TRAPEZOID LOOP
C         Adjust interval using Gaussian integration with 
C             Newton corrections since F is the derivative 
  520 CONTINUE
      X1 = XFCUM(IBIN)
      XBEST = X
      DTBEST = TPCTIL
      TPART = TPCTIL
C         Allow for maximum NITMAX more iterations on RADAPT
      DO 550 IHOME= 1, NITMAX
  535 XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
      X = XBEST + XINCR
      X2 = X
        IF (IHOME .GT. 1 .AND. X2 .EQ. XBEST) THEN
        WRITE (6,'(A,G12.3)')
     +  ' FUNLUX: WARNING from FUNPCT: insufficient precision at X=',X
        GO TO 580
        ENDIF
      REFX = ABS(X)+PRECIS
      CALL RADAPT(FUNC,X1,X2,1,RTEPS,0.,TPART2,UNCERT)
      DTPAR2 = TPART2-TPCTIL
      DTABS = ABS(DTPAR2)
      IF(ABS(XINCR)/REFX .LT. PRECIS) GOTO 545
      IF(DTABS .LT. DTBEST) GOTO 545
      XINCR = XINCR * 0.5
      GOTO 535
  545 DTBEST = DTABS
      XBEST = X
      TPART = TPART2
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
      IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
  550 CONTINUE
      WRITE (6,'(A,I4)') 
     +   ' FUNLUX: WARNING from FUNPCT: cannot converge, bin',IBIN
C
  580 CONTINUE
      XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
      X = XBEST + XINCR
      XFCUM(IBIN+1) = X
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
  600 CONTINUE
C         END OF LOOP OVER BINS
      X1 = XFCUM(NLO+NBINS-1)
      X2 = XHIGH
      CALL RADAPT(FUNC,X1,X2,1,RTEPS,0.,TPART ,UNCERT)
      ABERR = ABS(TPART-TPCTIL)/TFTOT
C      WRITE(6,1001) IFUNC,XLOW,XHIGH
      IF(ABERR .GT. RTEPS)  WRITE(6,1002) ABERR
      RETURN
  900 WRITE(6,1000) X,F
      IERR = 1
      RETURN
 1000 FORMAT(/' FUNLUX fatal error in FUNPCT: function negative:'/
     + ,' at X=',E15.6,', F=',E15.6/)
C 1001 FORMAT(' FUNPCT has prepared USER FUNCTION',I3,
C     + ' between',G12.3,' and',G12.3,' for FUNLUX.')
 1002 FORMAT(' WARNING: Relative error in cumulative distribution',
     + ' may be as big as',F10.7)
      END
      SUBROUTINE FUNLUX(ARRAY,XRAN,LEN)
C         Generation of LEN random numbers in any given distribution,
C         by 4-point interpolation in the inverse cumulative distr.
C         which was previously generated by FUNLXP
      COMMON/FUNINT/XUNI
      DIMENSION ARRAY(200)
      DIMENSION XRAN(LEN)
C  Bin width for main sequence, and its inverse
      PARAMETER (GAP= 1./99.,  GAPINV=99.)
C  Top of left tail, bottom of right tail (each tail replaces 2 bins)
      PARAMETER (TLEFT= 2./99.,BRIGHT=97./99.)
C  Bin width for minor sequences (tails), and its inverse
      PARAMETER (GAPS=TLEFT/49.,  GAPINS=1./GAPS)
C
C   The array ARRAY is assumed to have the following structure: 
C        ARRAY(1-100) contains the 99 bins of the inverse cumulative
C                     distribution of the entire function.
C        ARRAY(101-150) contains the 49-bin blowup of main bins 
C                       1 and 2 (left tail of distribution)
C        ARRAY(151-200) contains the 49-bin blowup of main bins 
C                       98 and 99 (right tail of distribution)
C
      CALL RANLUX(XRAN,LEN)
      
      DO 500 IBUF= 1, LEN
      X = XRAN(IBUF)
      J = INT(  X    *GAPINV) + 1
      IF (J .LT. 3)  THEN
         J1 = INT( X *GAPINS)
             J = J1 + 101
             J = MAX(J,102)
             J = MIN(J,148)
         P = (   X -GAPS*(J1-1)) * GAPINS
         A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
         B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
         XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
      ELSE IF (J .GT. 97)  THEN
         J1 = INT((X-BRIGHT)*GAPINS) 
             J = J1 + 151
             J = MAX(J,152)
             J = MIN(J,198)
         P = (X -BRIGHT -GAPS*(J1-1)) * GAPINS
         A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
         B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
         XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
      ELSE
CC      J = MAX(J,2)
CC      J = MIN(J,98)
         P = (   X -GAP*(J-1)) * GAPINV
         A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
         B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
         XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
      ENDIF
  500 CONTINUE
      XUNI = X
      RETURN
      END
      SUBROUTINE FUNLZ(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
C         FIND RANGE WHERE FUNC IS NON-ZERO.
C         WRITTEN 1980, F. JAMES
C         MODIFIED, NOV. 1985, TO FIX BUG AND GENERALIZE
C         TO FIND SIMPLY-CONNECTED NON-ZERO REGION (XLOW,XHIGH)
C         ANYWHERE WITHIN THE GIVEN REGION (X2LOW,H2HIGH).
C            WHERE 'ANYWHERE' MEANS EITHER AT THE LOWER OR UPPER
C            EDGE OF THE GIVEN REGION, OR, IF IN THE MIDDLE,
C            COVERING AT LEAST 1% OF THE GIVEN REGION.
C         OTHERWISE IT IS NOT GUARANTEED TO FIND THE NON-ZERO REGION.
C         IF FUNCTION EVERYWHERE ZERO, FUNLZ SETS XLOW=XHIGH=0.
      EXTERNAL FUNC
      XLOW = X2LOW
      XHIGH = X2HIGH
C         FIND OUT IF FUNCTION IS ZERO AT ONE END OR BOTH
      XMID = XLOW
      IF (FUNC(XLOW) .GT. 0.) GO TO 120
      XMID = XHIGH
      IF (FUNC(XHIGH) .GT. 0.)  GO TO 50
C         FUNCTION IS ZERO AT BOTH ENDS,
C         LOOK FOR PLACE WHERE IT IS NON-ZERO.
      DO 30 LOGN= 1, 7
      NSLICE = 2**LOGN
      DO 20 I= 1, NSLICE, 2
      XMID = XLOW + I * (XHIGH-XLOW) / NSLICE
      IF (FUNC(XMID) .GT. 0.)  GO TO 50
   20 CONTINUE
   30 CONTINUE
C         FALLING THROUGH LOOP MEANS CANNOT FIND NON-ZERO VALUE
      WRITE(6,554)
      WRITE(6,555) XLOW, XHIGH
      XLOW = 0.
      XHIGH = 0.
      GO TO 220
C
   50 CONTINUE
C         DELETE 'LEADING' ZERO RANGE
      XH = XMID
      XL = XLOW
      DO 70 K= 1, 20
      XNEW = 0.5*(XH+XL)
      IF (FUNC(XNEW) .EQ. 0.) GO TO 68
      XH = XNEW
      GO TO 70
   68 XL = XNEW
   70 CONTINUE
      XLOW = XL
      WRITE(6,555) X2LOW,XLOW
  120 CONTINUE
      IF (FUNC(XHIGH) .GT. 0.) GO TO 220
C         DELETE 'TRAILING' RANGE OF ZEROES
      XL = XMID
      XH = XHIGH
      DO 170 K= 1, 20
      XNEW = 0.5*(XH+XL)
      IF (FUNC(XNEW) .EQ. 0.) GO TO 168
      XL = XNEW
      GO TO 170
  168 XH = XNEW
  170 CONTINUE
      XHIGH = XH
      WRITE(6,555) XHIGH, X2HIGH
C
  220 CONTINUE
      RETURN
  554 FORMAT('0CANNOT FIND NON-ZERO FUNCTION VALUE')
  555 FORMAT(' FUNCTION IS ZERO FROM X=',E12.5,' TO ',E12.5)
      END
Commit	Line	Data
fe4da5cc	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) = AP(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) = AP(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) = AP(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