source: ZHANGProjects/ICOSIM/CPP/trunk/source/crystalprotonfuns.f @ 17

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  Math functions used in crystal_dan_FINAL_VERSION_CHvsCV.f
C      From track_dan_new.f
C
C
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC



      double precision function RAN_GAUSS(cut)
!*********************************************************************
!
! RAN_GAUSS - will generate a normal distribution from a uniform
!   distribution between [0,1].
!   See "Communications of the ACM", V. 15 (1972), p. 873.
!
! cut - double precision - cut for distribution in units of sigma
!                the cut must be greater than 0.5
!
!*********************************************************************
      implicit none
 
      logical flag
      real rndm4
      double precision x, u1, u2, twopi, r,cut
      save
 
c      write(*,*) "In RAN_GAUSS: "
c      write(*,*) "flag = ",flag
c      write(*,*) "U1 = ", U1,"U2 = ", U2
            twopi=8d0*atan(1d0)
    1       if (flag) then
              r = dble(rndm4( ))
              r = max(r, 0.5d0**32)
              r = min(r, 1d0-0.5d0**32)
              u1 = sqrt(-2d0*log( r ))
              u2 = dble(rndm4( ))
              x = u1 * cos(twopi*u2)
            else
              x = u1 * sin(twopi*u2)
            endif
 
c            write(*,*) "flag = ",flag
c            write(*,*) "U1 = ", U1,"U2 = ", U2, "x = ", x

          flag = .not. flag
 
!  cut the distribution if cut > 0.5
          if (cut .gt. 0.5d0 .and. abs(x) .gt. cut) goto 1
 
          RAN_GAUSS = x

c          write(*,*) "ran_gauss = ", RAN_GAUSS

        return
      end

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      function rndm4()
      implicit none
      integer len, in
      real rndm4, a
      save
      parameter ( len =  30000 )
      dimension a(len)
      data in/1/
!
      if ( in.eq.1 ) then
         call ranlux(a,len)
         rndm4=a(1)
         in=2
!        write(6,'('' LEN: '',i5)')LEN
      else
         rndm4=a(in)
         in=in+1
         if(in.eq.len+1)in=1
      endif
      return
      end



      subroutine ranlux(rvec,lenv)
!         Subtract-and-borrow random number generator proposed by
!         Marsaglia and Zaman, implemented by F. James with the name
!         RCARRY in 1991, and later improved by Martin Luescher
!         in 1993 to produce "Luxury Pseudorandom Numbers".
!     Fortran 77 coded by F. James, 1993
!
!   LUXURY LEVELS.
!   ------ ------      The available luxury levels are:
!
!  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
!           and Zaman, very long period, but fails many tests.
!  level 1  (p=48): considerable improvement in quality over level 0,
!           now passes the gap test, but still fails spectral test.
!  level 2  (p=97): passes all known tests, but theoretically still
!           defective.
!  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
!           correlations have very small chance of being observed.
!  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
!
!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!!!  Calling sequences for RANLUX:                                  ++
!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
!!!                   32-bit random floating point numbers between  ++
!!!                   zero (not included) and one (also not incl.). ++
!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
!!!               which is integer between zero and MAXLEV, or if   ++
!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
!!!               should be set to zero unless restarting at a break++
!!!               point given by output of RLUXAT (see RLUXAT).     ++
!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
!!!               which can be used to restart the RANLUX generator ++
!!!               at the current point by calling RLUXGO.  K1 and K2++
!!!               specify how many numbers were generated since the ++
!!!               initialization with LUX and INT.  The restarting  ++
!!!               skips over  K1+K2*E9   numbers, so it can be long.++
!!!   A more efficient but less convenient way of restarting is by: ++
!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
!!!                 32-bit integer seeds, to be used for restarting ++
!!!      ISVEC must be dimensioned 25 in the calling program        ++
!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit none
      integer lenv,isdext,iseeds,maxlev,ndskip,itwo24,next,j24,i24,     &
     &inseed,mkount,kount,in24,nskip,lxdflt,jsdflt,jseed,lp,i,k,icons,  &
     &inner,izip,izip2,ivec,isk,igiga,isd,k2,k1,inout,lout,ins,lux,ilx, &
     &iouter
      real rvec,seeds,twop12,twom12,twom24,carry,uni
      dimension rvec(lenv)
      dimension seeds(24), iseeds(24), isdext(25)
      parameter (maxlev=4, lxdflt=3)
      dimension ndskip(0:maxlev)
      dimension next(24)
      parameter (twop12=4096., igiga=1000000000,jsdflt=314159265)
      parameter (itwo24=2**24, icons=2147483563)
      save notyet, i24, j24, carry, seeds, twom24, twom12, luxlev
      save nskip, ndskip, in24, next, kount, mkount, inseed
      integer luxlev
      logical notyet
      data notyet, luxlev, in24, kount, mkount /.true., lxdflt, 0,0,0/
      data i24,j24,carry/24,10,0./
!                               default
!  Luxury Level   0     1     2   *3*    4
      data ndskip/0,   24,   73,  199,  365 /
!Corresponds to p=24    48    97   223   389
!     time factor 1     2     3     6    10   on slow workstation
!                 1    1.5    2     3     5   on fast mainframe
!
!  NOTYET is .TRUE. if no initialization has been performed yet.
!              Default Initialization by Multiplicative Congruential
      if (notyet) then
         notyet = .false.
         jseed = jsdflt
         inseed = jseed
         write(*,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',jseed
         luxlev = lxdflt
         nskip = ndskip(luxlev)
         lp = nskip + 24
         in24 = 0
         kount = 0
         mkount = 0
!         WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
!     &        LUXLEV,'      p =',LP
            twom24 = 1.
         do 25 i= 1, 24
            twom24 = twom24 * 0.5
         k = jseed/53668
         jseed = 40014*(jseed-k*53668) -k*12211
         if (jseed .lt. 0)  jseed = jseed+icons
         iseeds(i) = mod(jseed,itwo24)
   25    continue
         twom12 = twom24 * 4096.
         do 50 i= 1,24
         seeds(i) = real(iseeds(i))*twom24
         next(i) = i-1
   50    continue
         next(1) = 24
         i24 = 24
         j24 = 10
         carry = 0.
         if (seeds(24) .eq. 0.) carry = twom24
      endif
!
!          The Generator proper: "Subtract-with-borrow",
!          as proposed by Marsaglia and Zaman,
!          Florida State University, March, 1989
!
      do 100 ivec= 1, lenv
      uni = seeds(j24) - seeds(i24) - carry
      if (uni .lt. 0.)  then
         uni = uni + 1.
         carry = twom24
      else
         carry = 0.
      endif
      seeds(i24) = uni
      i24 = next(i24)
      j24 = next(j24)
      rvec(ivec) = uni
!  small numbers (with less than 12 "significant" bits) are "padded".
      if (uni .lt. twom12)  then
         rvec(ivec) = rvec(ivec) + twom24*seeds(j24)
!        and zero is forbidden in case someone takes a logarithm
         if (rvec(ivec) .eq. 0.)  rvec(ivec) = twom24*twom24
      endif
!        Skipping to luxury.  As proposed by Martin Luscher.
      in24 = in24 + 1
      if (in24 .eq. 24)  then
         in24 = 0
         kount = kount + nskip
         do 90 isk= 1, nskip
         uni = seeds(j24) - seeds(i24) - carry
         if (uni .lt. 0.)  then
            uni = uni + 1.
            carry = twom24
         else
            carry = 0.
         endif
         seeds(i24) = uni
         i24 = next(i24)
         j24 = next(j24)
   90    continue
      endif
  100 continue
      kount = kount + lenv
      if (kount .ge. igiga)  then
         mkount = mkount + 1
         kount = kount - igiga
      endif
      return
!
!           Entry to input and float integer seeds from previous run
      entry rluxin(isdext)
         notyet = .false.
         twom24 = 1.
         do 195 i= 1, 24
         next(i) = i-1
  195    twom24 = twom24 * 0.5
         next(1) = 24
         twom12 = twom24 * 4096.
      write(*,*) ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
      write(*,'(5X,5I12)') isdext
      do 200 i= 1, 24
      seeds(i) = real(isdext(i))*twom24
  200 continue
      carry = 0.
      if (isdext(25) .lt. 0)  carry = twom24
      isd = iabs(isdext(25))
      i24 = mod(isd,100)
      isd = isd/100
      j24 = mod(isd,100)
      isd = isd/100
      in24 = mod(isd,100)
      isd = isd/100
      luxlev = isd
        if (luxlev .le. maxlev) then
          nskip = ndskip(luxlev)
          write(*,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',  &
     &luxlev
        else  if (luxlev .ge. 24) then
          nskip = luxlev - 24
          write(*,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',luxlev
        else
          nskip = ndskip(maxlev)
          write(*,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',luxlev
          luxlev = maxlev
        endif
      inseed = -1
      return
!
!                    Entry to ouput seeds as integers
      entry rluxut(isdext)
      do 300 i= 1, 24
         isdext(i) = int(seeds(i)*twop12*twop12)
  300 continue
      isdext(25) = i24 + 100*j24 + 10000*in24 + 1000000*luxlev
      if (carry .gt. 0.)  isdext(25) = -isdext(25)
      return
!
!                    Entry to output the "convenient" restart point
      entry rluxat(lout,inout,k1,k2)
      lout = luxlev
      inout = inseed
      k1 = kount
      k2 = mkount
      return
!
!                    Entry to initialize from one or three integers
      entry rluxgo(lux,ins,k1,k2)
         if (lux .lt. 0) then
            luxlev = lxdflt
         else if (lux .le. maxlev) then
            luxlev = lux
         else if (lux .lt. 24 .or. lux .gt. 2000) then
            luxlev = maxlev
            write(*,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',lux
         else
            luxlev = lux
            do 310 ilx= 0, maxlev
              if (lux .eq. ndskip(ilx)+24)  luxlev = ilx
  310       continue
         endif
      if (luxlev .le. maxlev)  then
         nskip = ndskip(luxlev)
         write(*,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', &
     &luxlev,'     P=', nskip+24
      else
          nskip = luxlev - 24
          write(*,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',luxlev
      endif
      in24 = 0
      if (ins .lt. 0)  write(*,*)                                       &
     &' Illegal initialization by RLUXGO, negative input seed'
      if (ins .gt. 0)  then
        jseed = ins
        write(*,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', &
     &jseed, k1,k2
      else
        jseed = jsdflt
        write(*,*)' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
      endif
      inseed = jseed
      notyet = .false.
      twom24 = 1.
         do 325 i= 1, 24
           twom24 = twom24 * 0.5
         k = jseed/53668
         jseed = 40014*(jseed-k*53668) -k*12211
         if (jseed .lt. 0)  jseed = jseed+icons
         iseeds(i) = mod(jseed,itwo24)
  325    continue
      twom12 = twom24 * 4096.
         do 350 i= 1,24
         seeds(i) = real(iseeds(i))*twom24
         next(i) = i-1
  350    continue
      next(1) = 24
      i24 = 24
      j24 = 10
      carry = 0.
      if (seeds(24) .eq. 0.) carry = twom24
!        If restarting at a break point, skip K1 + IGIGA*K2
!        Note that this is the number of numbers delivered to
!        the user PLUS the number skipped (if luxury .GT. 0).
      kount = k1
      mkount = k2
      if (k1+k2 .ne. 0)  then
        do 500 iouter= 1, k2+1
          inner = igiga
          if (iouter .eq. k2+1)  inner = k1
          do 450 isk= 1, inner
            uni = seeds(j24) - seeds(i24) - carry
            if (uni .lt. 0.)  then
               uni = uni + 1.
               carry = twom24
            else
               carry = 0.
            endif
            seeds(i24) = uni
            i24 = next(i24)
            j24 = next(j24)
  450     continue
  500   continue
!         Get the right value of IN24 by direct calculation
        in24 = mod(kount, nskip+24)
        if (mkount .gt. 0)  then
           izip = mod(igiga, nskip+24)
           izip2 = mkount*izip + in24
           in24 = mod(izip2, nskip+24)
        endif
!       Now IN24 had better be between zero and 23 inclusive
        if (in24 .gt. 23) then
           write(*,'(A/A,3I11,A,I5)')                                   &
     &'  Error in RESTARTING with RLUXGO:','  The values', ins,         &
     &k1, k2, ' cannot occur at luxury level', luxlev
           in24 = 0
        endif
      endif
      return
      end
 
!cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine funlxp (func,xfcum,x2low,x2high)
!         F. JAMES,   Sept, 1994
!
!         Prepares the user function FUNC for FUNLUX
!         Inspired by and mostly copied from FUNPRE and FUNRAN
!         except that
!    1. FUNLUX uses RANLUX underneath,
!    2. FUNLXP expands the first and last bins to cater for
!              functions with long tails on left and/or right,
!    3. FUNLXP calls FUNPCT to do the actual finding of percentiles.
!    4. both FUNLXP and FUNPCT use RADAPT for Gaussian integration.
!
      implicit none
      external func
      integer ifunc,ierr
      real x2high,x2low,xfcum,rteps,xhigh,xlow,xrange,uncert,x2,tftot1, &
     &x3,tftot2,func
      real tftot
      common/funint/tftot
      dimension xfcum(200)
      parameter (rteps=0.0002)
      save ifunc
      data ifunc/0/
      ifunc = ifunc + 1
!         FIND RANGE WHERE FUNCTION IS NON-ZERO.
      call funlz(func,x2low,x2high,xlow,xhigh)
      xrange = xhigh-xlow
      if(xrange .le. 0.)  then
        write(*,'(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)
!
!      WRITE (6,'(A,A)') ' FUNLXP preparing ',
!     + '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(*,*) ' Fatal error in FUNLXP. FUNLUX will not work.'
      end
!


      subroutine funlux(array,xran,len)
!         Generation of LEN random numbers in any given distribution,
!         by 4-point interpolation in the inverse cumulative distr.
!         which was previously generated by FUNLXP
      implicit none
      real tftot
      common/funint/tftot
      integer len,ibuf,j,j1
      real array,xran,gap,gapinv,tleft,bright,gaps,gapins,x,p,a,b
      dimension array(200)
      dimension xran(len)
!  Bin width for main sequence, and its inverse
      parameter (gap= 1./99.,  gapinv=99.)
!  Top of left tail, bottom of right tail (each tail replaces 2 bins)
      parameter (tleft= 2./99.,bright=97./99.)
!  Bin width for minor sequences (tails), and its inverse
      parameter (gaps=tleft/49.,  gapins=1./gaps)
!
!   The array ARRAY is assumed to have the following structure:
!        ARRAY(1-100) contains the 99 bins of the inverse cumulative
!                     distribution of the entire function.
!        ARRAY(101-150) contains the 49-bin blowup of main bins
!                       1 and 2 (left tail of distribution)
!        ARRAY(151-200) contains the 49-bin blowup of main bins
!                       98 and 99 (right tail of distribution)
!
      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.)*(p-2.)*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.)*(p-2.)*0.5
      else
!      J = MAX(J,2)
!      J = MIN(J,98)
         p = (   x -gap*(j-1)) * gapinv
         a = (p+1.) * array(j+2) - (p-2.)*array(j-1)
         b = (p-1.) * array(j) - p * array(j+1)
         xran(ibuf) = a*p*(p-1.)*0.16666667 + b*(p+1.)*(p-2.)*0.5
      endif
  500 continue
      tftot = x
      return
      end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine funlz(func,x2low,x2high,xlow,xhigh)
!         FIND RANGE WHERE FUNC IS NON-ZERO.
!         WRITTEN 1980, F. JAMES
!         MODIFIED, NOV. 1985, TO FIX BUG AND GENERALIZE
!         TO FIND SIMPLY-CONNECTED NON-ZERO REGION (XLOW,XHIGH)
!         ANYWHERE WITHIN THE GIVEN REGION (X2LOW,H2HIGH).
!            WHERE 'ANYWHERE' MEANS EITHER AT THE LOWER OR UPPER
!            EDGE OF THE GIVEN REGION, OR, IF IN THE MIDDLE,
!            COVERING AT LEAST 1% OF THE GIVEN REGION.
!         OTHERWISE IT IS NOT GUARANTEED TO FIND THE NON-ZERO REGION.
!         IF FUNCTION EVERYWHERE ZERO, FUNLZ SETS XLOW=XHIGH=0.
      implicit none
      external func
      integer logn,nslice,i,k
      real xhigh,xlow,x2high,x2low,func,xmid,xh,xl,xnew
      xlow = x2low
      xhigh = x2high
!         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
!         FUNCTION IS ZERO AT BOTH ENDS,
!         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
!         FALLING THROUGH LOOP MEANS CANNOT FIND NON-ZERO VALUE
      write(*,554)
      write(*,555) xlow, xhigh
      xlow = 0.
      xhigh = 0.
      go to 220
!
   50 continue
!         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(*,555) x2low,xlow
  120 continue
      if (func(xhigh) .gt. 0.) go to 220
!         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(*,555) xhigh, x2high
!
  220 continue
      return
  554 format('0CANNOT FIND NON-ZERO FUNCTION VALUE')
  555 format(' FUNCTION IS ZERO FROM X=',e12.5,' TO ',e12.5)
      end


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
! $Id: radapt.F,v 1.1.1.1 1996/04/01 15:02:13 mclareni Exp $
!
! $Log: radapt.F,v $
! Revision 1.1.1.1  1996/04/01 15:02:13  mclareni
! Mathlib gen
!
!
      subroutine radapt(f,a,b,nseg,reltol,abstol,res,err)
 
!     RES = Estimated Integral of F from A to B,
!     ERR = Estimated absolute error on RES.
!     NSEG  specifies how the adaptation is to be done:
!        =0   means use previous binning,
!        =1   means fully automatic, adapt until tolerance attained.
!        =n>1 means first split interval into n equal segments,
!             then adapt as necessary to attain tolerance.
!     The specified tolerances are:
!            relative: RELTOL ;  absolute: ABSTOL.
!        It stop s when one OR the other is satisfied, or number of
!        segments exceeds NDIM.  Either TOLA or TOLR (but not both!)
!        can be set to zero, in which case only the other is used.
 
      implicit none
      external f
      integer nseg,ndim,nter,nsegd,i,iter,ibig
      real err,res,abstol,reltol,b,a,xlo,xhi,tval,ters,te,root,xhib,    &
     &bin,xlob,bige,hf,xnew,r1,f
      double precision tvals,terss
 
      parameter (ndim=100)
      parameter (r1 = 1., hf = r1/2.)
 
      dimension xlo(ndim),xhi(ndim),tval(ndim),ters(ndim)
      save xlo,xhi,tval,ters,nter
      data nter /0/
 
      if(nseg .le. 0)  then
       if(nter .eq. 0) then
        nsegd=1
        go to 2
       endif
       tvals=0d0
       terss=0d0
       do 1 i = 1,nter
       call rgs56p(f,xlo(i),xhi(i),tval(i),te)
       ters(i)=te**2
       tvals=tvals+tval(i)
       terss=terss+ters(i)
    1  continue
       root= sqrt(2.*terss)
       go to 9
      endif
      nsegd=min(nseg,ndim)
    2 xhib=a
      bin=(b-a)/nsegd
      do 3 i = 1,nsegd
      xlo(i)=xhib
      xlob=xlo(i)
      xhi(i)=xhib+bin
      if(i .eq. nsegd) xhi(i)=b
      xhib=xhi(i)
      call rgs56p(f,xlob,xhib,tval(i),te)
      ters(i)=te**2
    3 continue
      nter=nsegd
      do 4 iter = 1,ndim
      tvals=tval(1)
      terss=ters(1)
      do 5 i = 2,nter
      tvals=tvals+tval(i)
      terss=terss+ters(i)
    5 continue
      root= sqrt(2.*terss)
      if(root .le. abstol .or. root .le. reltol*abs(tvals)) go to 9
      if(nter .eq. ndim) go to 9
      bige=ters(1)
      ibig=1
      do 6 i = 2,nter
      if(ters(i) .gt. bige) then
       bige=ters(i)
       ibig=i
      endif
    6 continue
      nter=nter+1
      xhi(nter)=xhi(ibig)
      xnew=hf*(xlo(ibig)+xhi(ibig))
      xhi(ibig)=xnew
      xlo(nter)=xnew
      call rgs56p(f,xlo(ibig),xhi(ibig),tval(ibig),te)
      ters(ibig)=te**2
      call rgs56p(f,xlo(nter),xhi(nter),tval(nter),te)
      ters(nter)=te**2
    4 continue
    9 res=tvals
      err=root
      return
      end


!cccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
!
! $Id: rgs56p.F,v 1.1.1.1 1996/04/01 15:02:14 mclareni Exp $
!
! $Log: rgs56p.F,v $
! Revision 1.1.1.1  1996/04/01 15:02:14  mclareni
! Mathlib gen
!
!
      subroutine rgs56p(f,a,b,res,err)
      implicit none
      integer i
      real err,res,b,a,f,w6,x6,w5,x5,rang,r1,hf
      double precision e5,e6
 
      parameter (r1 = 1., hf = r1/2.)
      dimension x5(5),w5(5),x6(6),w6(6)
 
      data (x5(i),w5(i),i=1,5)                                          &
     &/4.6910077030668004e-02, 1.1846344252809454e-01,                  &
     &2.3076534494715846e-01, 2.3931433524968324e-01,                   &
     &5.0000000000000000e-01, 2.8444444444444444e-01,                   &
     &7.6923465505284154e-01, 2.3931433524968324e-01,                   &
     &9.5308992296933200e-01, 1.1846344252809454e-01/
 
      data (x6(i),w6(i),i=1,6)                                          &
     &/3.3765242898423989e-02, 8.5662246189585178e-02,                  &
     &1.6939530676686775e-01, 1.8038078652406930e-01,                   &
     &3.8069040695840155e-01, 2.3395696728634552e-01,                   &
     &6.1930959304159845e-01, 2.3395696728634552e-01,                   &
     &8.3060469323313225e-01, 1.8038078652406930e-01,                   &
     &9.6623475710157601e-01, 8.5662246189585178e-02/
 
      rang=b-a
      e5=0d0
      e6=0d0
      do 1 i = 1,5
      e5=e5+w5(i)*f(a+rang*x5(i))
      e6=e6+w6(i)*f(a+rang*x6(i))
    1 continue
      e6=e6+w6(6)*f(a+rang*x6(6))
      res=hf*(e6+e5)*rang
      err=abs((e6-e5)*rang)
      return
      end
!GRD


!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine funpct(func,ifunc,xlow,xhigh,xfcum,nlo,nbins,tftot,    &
     &ierr)
!        Array XFCUM is filled from NLO to NLO+NBINS, which makes
!        the number of values NBINS+1, or the number of bins NBINS
      implicit none
      external func
      integer ierr,nbins,nlo,ifunc,nz,ibin,maxz,iz,nitmax,ihome
      real tftot,xhigh,xlow,func,xfcum,rteps,tpctil,tz,tzmax,x,f,tcum,  &
     &x1,f1,dxmax,fmin,fminz,xincr,tincr,xbest,dtbest,tpart,x2,precis,  &
     &refx,uncert,tpart2,dtpar2,dtabs,aberr
      dimension xfcum(*)
      parameter (rteps=0.005, nz=10, maxz=20, nitmax=6,precis=1e-6)
!      DOUBLE PRECISION TPCTIL, TZ, TCUM, XINCR, DTABS,
!     &  TINCR, TZMAX, XBEST, DTBEST, DTPAR2
!
      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
!         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
!         Loop over trapezoids within a supposed percentil
      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(*,*) ' FUNLUX:  WARNING. FUNPCT fails trapezoid.'
!         END OF TRAPEZOID LOOP
!         Adjust interval using Gaussian integration with
!             Newton corrections since F is the derivative
  520 continue
      x1 = xfcum(ibin)
      xbest = x
      dtbest = tpctil
      tpart = tpctil
!         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(*,'(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(*,'(A,I4)')                                                 &
     &' FUNLUX: WARNING from FUNPCT: cannot converge, bin',ibin
!
  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
!         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
!      WRITE(6,1001) IFUNC,XLOW,XHIGH
      if(aberr .gt. rteps)  write(*,1002) aberr
      return
  900 write(*,1000) x,f
      ierr = 1
      return
 1000 format(/' FUNLUX fatal error in FUNPCT: function negative:'/      &
     &,' at X=',e15.6,', F=',e15.6/)
! 1001 FORMAT(' FUNPCT has prepared USER FUNCTION',I3,
!     + ' between',G12.3,' and',G12.3,' for FUNLUX.')
 1002 format(' WARNING: Relative error in cumulative distribution',     &
     &' may be as big as',f10.7)
      end