c-------------------------------------------------------------------- 
c         Find approximation for velocity dispersion
c         using maximum likelihood analysis of dV-dR data
c         Model: 
c              dN = [n(R)*2/sqrt(2pi s(R)**2)*exp(-(V/s(R))**2/2) +
c                   + n_0*R]dV*dR
c              Total Number = n_0*Vmax*Rmax**2/2. +
c                            + int(n(R)*erf(Vmax/(2*s(R)))
c         Functional = Sum_i(ln(N(i;n(R),s(R),n_0))) -
c                     - N*ln(TotalNumber)
      Include 'dvmax.h'
      external FUNC


      Open(1,file='VRdiagram.dat')
      Open(20,file='VRall.dat',position='append')
      Open(3,file='satellites.dat')
      Open(10,file='dvmax.config')

      CALL ReadConfiguration
      CALL ReadData
      CALL DensProfile

         a0     = 2.7e-6  ! initial guess
         dcoef(0) =3.5
         vcoef(0) =160.
         do i=1,Nterms   ! initial guess
             dcoef(i) =0.
             vcoef(i) =0.
         EndDo 
         dcoef(1) =-1.25
         vcoef(1) =-2.

         CALL FindMax
         fInterlop = a0*Vmax*(Rmax**2-Rmin**2)/2.
         Do j=0,2,2
            write (j,'(3(2x,a,f7.2))')' interlopers= ',fInterlop,
     &         'percentage= ',fInterlop/Nsample*100.,
     &         'True= ',Nsample- fInterlop
         EndDo 
      Do j=0,2,2
         write (j,'(a,T12,a,T25,a)')
     &    ' R(kpc)','Surf_dens','dV(km/s)  <= solution'
      EndDo 
      write (*,*)  ' a0=',a0/(2.*pi)*Vmax
      Do i=1,20      ! print new solution
         R = (i/20.)*Rmax
         V = sigma(R)
         D =  dens(R) /(2.*pi*R)
         Dp=  dens(R)*erf(Vmax/sq2/sigma(R))
         Do j=0,2,2
            write (j,'(f8.2,8g12.4)')R,D,V,Dp
         EndDo 
      EndDo 
         R100 = 100.            ! store results for this radius
         V100 = sigma(R100)
         D100 =  dens(R100) /(2.*pi*R100)
      CALL GetTests(Nseed,50)
         i= INT(R100/Rmax*20.)
         rr =(i/20.)*Rmax
         write (*,*)  ' R100=',R100,' must ==',rr
c      write(20,'(3a)')'  Mag Limits    Dmin Dmax(kms)  u-g colors',
c     &    '  Rmin Rmax(kpc) Vmax(kms)  Nsats R(kpc) <dV>    error',
c     &    '  SurfDen     error       L/Lsun      Mstarts/Msun'
      write (20,'(2f7.2,2f8.1,2f6.2,2f6.1,2x,f9.1,i7,3x,3f7.2,4g12.4)')
     &    aMag_min,aMag_max,Veloc_min,Veloc_max,urcolor_min,urcolor_max,
     &     dR_min,dR_max,Vmax,Nsample,R100,V100,SigV(i),D100,SigD(i),
     &     Flux,StMass 
      Stop
      End
c--------------------------------------------------------------------
c                         direct estimate of surface density
      SUBROUTINE DensProfile
c                      
      Include 'dvmax.h'

      dlog =0.07

      Do i=0,100
         Rprof(i) =0.
         Dprof(i) =0.
      EndDo 
      Do i=1,Nsample
         R = Rad(i)
         indx =max(min(INT(Log10(R)/dlog),100),0)
         Dprof(indx) =Dprof(indx) +1.
         Rprof(indx) =Rprof(indx) +R
      EndDo
      Do j=0,2,2
         write (j,'(a,T12,a,T18,a,T27,a,T40,a)')
     &    ' Rin(kpc)','Rmean','Nobjects','Surf_dens','rms'
      EndDo 
      Do i=1,99
         r1 = 10.**(i*dlog)
         r2 = 10.**((i+1)*dlog)
         r  = Rprof(i)/max(Dprof(i),1.d0)
         ss = Dprof(i)/(pi*(r2**2-r1**2))
         ers= sqrt(Dprof(i))/(pi*(r2**2-r1**2))
         if(r2.lt.Rmax.and.Dprof(i).gt.0.d0)Then
         Do j=0,2,2   
            write (j,'(2f8.2,i6,3g12.3)')r1,r,INT(Dprof(i)),ss,ers
         EndDo 
         EndIf 
      EndDo 

      return
      end
c--------------------------------------------------------------------
c            read and parse a line from configuration file
      SUBROUTINE ReadString(x1,x2)
c                      
      Include 'dvmax.h'
      Character Line*50,String*100,L2*100

      read(10,'(a)') String   ! read the string
      write (L2,*)String      ! place it to internal file L2
      read (L2,*,err=3)x1,x2  ! try to get 2 parameters from L2
         ind2 =0
         goto 4
 3    read (L2,*,err=3)x1     ! read only one number, if two failed
         x2   =0.
         ind2 =1
 4       ind1 =0
      do i=3,100              ! find second gap between numbers -> this is
                              !             the beginning of the comment
         if(L2(i-1:i).eq.'  ')ind1 =ind1 +1
         if(L2(i-1:i).ne.'  '.and.ind1.ne.0)Then
            ind1 = 0
            ind2 = ind2 +1
         EndIf 
         if(L2(i-1:i).ne.'  '.and.ind2.eq.2)Then
            ind1 = i
            goto 5
         EndIf 
      enddo
 5    Do i=100,1,-1           ! find position of the last character 
         if(L2(i:i).ne.' ')Then
            ind2 = i
            goto10
         EndIf 
      EndDo 
 10   write(*,'(3x,50("."),T3,a,T50,1P,2g11.4)') L2(ind1:ind2),x1,x2
      write(1,'(3x,50("."),T3,a,T50,1P,2g11.4)') L2(ind1:ind2),x1,x2
      write(2,'(3x,50("."),T3,a,T50,1P,2g11.4)') L2(ind1:ind2),x1,x2

      return
      end
c--------------------------------------------------------------------
c                          read configuration parameters
c                          set criteria for selection
c                          open output file '2'
      SUBROUTINE ReadConfiguration
c                      
      Include 'dvmax.h'
      Character Line*40,FileName*80

      read(10,'(a)') FileName   ! read name of output file
      Open(2,file=FileName)      
      CALL  ReadString( aMag_min,aMag_max)
      CALL  ReadString( Veloc_min,Veloc_max)
      CALL  ReadString( qpar_min,qpar_max)
      CALL  ReadString( urcolor_min,urcolor_max)
      CALL  ReadString( dR_min,dR_max)
      CALL  ReadString( Vmax,dummy)

      Rmax = dR_max
      Rmin = dR_min
      Return
      End
c--------------------------------------------------------------------
c                          read table of SDSS data
      SUBROUTINE ReadData
c                      
      Include 'dvmax.h'
      Character Line*80

      Gsun  = 5.07
      Lines = 5
      Flux  = 0.
      StMass= 0.
      Do i=1,Lines
         read(3,'(a)') Line
         write(*,'(a)') Line
      EndDo 
      Nsample =0
 10   read(3,*,err=100,end=100) n1,n2,aMag,urcolor,stellmass,
     &       aMagSat,x,v,dmag,qparam,Veloc
      If(aMag.lt.aMag_min.and.aMag.gt.aMag_max
     &       .and.Veloc.gt.Veloc_min.and.Veloc.lt.Veloc_max
     &       .and.qparam.gt.qpar_min.and.qparam.lt.qpar_max
     &       .and.urcolor.gt.urcolor_min.and.urcolor.lt.urcolor_max
     &       .and.abs(v).lt.Vmax
     &       .and.x.lt.dR_max.and.x.gt.dR_min)Then
      If(stellmass.gt.13.)
     &       write (*,'(2i8,11g11.3)')  n1,n2,aMag,urcolor,stellmass,
     &    aMagSat,x,v,qparam,Veloc

         Nsample      = Nsample +1
         Rad(Nsample) = x
         Vp(Nsample)  = abs(v)
         Flux         = Flux +10.**(-0.4*(aMag-Gsun))
         StMass       = StMass + 10.**(stellmass)
      write (1,'(i5,2f9.2,2x,7f9.3)')Nsample,aMag,urcolor,Veloc,
     &                            Rad(Nsample),  Vp(Nsample)
      EndIf  
      goto 10

 100  Flux   = Flux/Nsample
      stMass = StMass/Nsample
      Do j=0,2
         write(j,'(1P,3(2x,a,g12.4))')'Luminosity/Lsun=',Flux,
     &          ' Stellar Mass/Msun=',stMass,' M/L=',stMass/Flux
      EndDo 
      Do j=0,2,2
         write (j,'(3x,50("."),T3,a,T50,i6)')'Nsample',Nsample
      EndDo 

      return
      end
c--------------------------------------------------------------------
c        Generates many tests for distribution of dV-dR diagram
c
c
      SUBROUTINE GetTests(Nseed,Niter)
c                      
      Include 'dvmax.h'
      DIMENSION dcoef2(0:Ncoef),vcoef2(0:Ncoef)

      Nseed= 121071

      a1 = a0               ! store previous solution
      Do i=0,Ncoef
         dcoef2(i) =dcoef(i)
         vcoef2(i) =vcoef(i)
      EndDo 
      Do i=0,100
         Rprof(i) =0.
         Vprof(i) =0.
         Dprof(i) =0.
         SigV(i) =0.
         SigD(i) =0.
      EndDo 

      Do iter =0,Niter         ! make realizations and solve each
         a0 = a1               ! re-store previous solution
         Do i=0,Ncoef
            dcoef(i) =dcoef2(i)
            vcoef(i) =vcoef2(i)
         EndDo 

         CALL GetTest(Nseed,iter)
         CALL FindMax
      EndDo 
       Do j=0,2,2
         write (j,'(a,T14,a,T35,a)')
     &   ' R(kpc)','Surf_dens  Error','dV(km/s)    Error <= Monte Carlo'
      EndDo 
      
      Do i=1,20      ! print new solution
         R = (i/20.)*Rmax
         nn = nIter+1
         Rprof(i) =Rprof(i)/nn
         Vprof(i) =Vprof(i)/nn
         Dprof(i) =Dprof(i)/nn
         SigV(i)  =sqrt(SigV(i)/nn-Vprof(i)**2)
         SigD(i)  =sqrt(SigD(i)/nn-Dprof(i)**2)
         Dprof(i) =Dprof(i)/(2.*pi* Rprof(i))
         SigD(i)  =SigD(i)/(2.*pi* Rprof(i))

         write (*,'(5g11.4)')Rprof(i),Dprof(i),SigD(i),Vprof(i),SigV(i)
         write (2,'(5g11.4)')Rprof(i),Dprof(i),SigD(i),Vprof(i),SigV(i)
      EndDo 

      Return
      end
c--------------------------------------------------------------------
c              Generates a test distribution of dV-dR diagram
c              Use Chebyshev polynomials up to order Nterms
c              Set n(R) and sigma_v(R)
c              Make a realization
c              Maxima of  both n(R) and sigma_v(R) are at R=0    
      SUBROUTINE GetTest(Nseed,indx)
c                      
      Include 'dvmax.h'
      EXTERNAL FUNC
c      Do i=0,Ncoef
c         dcoef(i) =0.
c         vcoef(i) =0.
c      EndDo 
c      dcoef(0) =1.
c      dcoef(1) = -0.8*dcoef(0)
c      dcoef(2) =  0.05*dcoef(0)
c      vcoef(0) =160.
c      vcoef(1) =  -0.3*vcoef(0)
c      a0       =2.3d-6      ! defines number of interlopers

      ERR      =1.d-6
      Ninterlop = a0*Vmax*(Rmax**2-Rmin**2)/2.

      CALL DGAUS8 (FUNC, Rmin, Rmax, ERR, ANS, IERR)
      Nsat     = ANS
      Nsample  = Nsat + Ninterlop
      If(indx.le.0)Then
        Do ifile=0,2,2
         write (ifile,'(3(2x,a,i6))')'Test: number of true sats=',Nsat,
     &            'Niterlop=',Ninterlop,' Ntot=',Nsample
         write(ifile,'(a,14g13.4)')' Initial coeff=',
     &             (dcoef(i),i=0,Nterms),
     &             (vcoef(i),i=0,Nterms),a0
        EndDo 
c        Do i=0,20
c          R = (i/20.)*Rmax
c          write (2,'(a,f8.2,3(a,g11.3))')' Rad=',R,' dens=',dens(R),
cc     &          ' velocity=',sigma(R)
c          write (*,'(a,f8.2,3(a,g11.3))')' Rad=',R,' dens=',dens(R),
c     &          ' velocity=',sigma(R)
c        EndDo 
      EndIf 
      dens_max = dens(Rmin)
      Do i=1,1000
         R  =Rmin+(Rmax-Rmin)/1000.*i
         dens_max = max(dens(R),dens_max)
      EndDo 
      If(Nsat.gt.0)Then
         Do i=1,Nsat
            n = 0
 1          R = (Rmax-Rmin)*RANDd(Nseed)+Rmin
            n = n +1
            if(n.gt.1000)Stop ' too many iterations'
c            write (*,'(2i5,3g11.3)')  i,n,R,dens(R),dens_max
            If(dens(R)/dens_max.lt.RANDd(Nseed))Goto1
            Rad(i) = R
            Vp(i)  = sigma(R)*GAUSS(Nseed)
c            write (1,'(2i5,7f9.3)')i,n,Rad(i),Vp(i),dens(R),sigma(R)
         EndDo 
      EndIf                      ! end true satellites
       If(Ninterlop.gt.0)Then
         Do i=1,Ninterlop
            n = 0
 2          R = (Rmax-Rmin)*RANDd(Nseed)+Rmin
            n = n +1
            if(n.gt.1000)Stop ' too many iterations'
c            write (*,'(2i5,3g11.3)')  i,n,R,dens(R),dens_max
            If(R/Rmax.lt.RANDd(Nseed))Goto2
            Rad(i+Nsat) = R
            Vp(i+Nsat)  = Vmax*RANDd(Nseed)
c            write (1,'(2i5,7f9.3)')i,n,Rad(i+Nsat),Vp(i+Nsat)
         EndDo 
      EndIf                      ! end interlopers
      return
      end
c--------------------------------------------------------------------
c                    find maximum of likelihood function
      SUBROUTINE FindMax
c                    pmin and pmax are arrays of current limits 
c                            for search of maximum Likelihood
c                    pmax  are the current steps
c      Mapping of pmin/max/step to parameters of Cheb.approximation:
c         1         2     ...  Nterms+2 ...   2*Nterms+3 <= pmin/max/step indexes
c      dcoef(0)  dcoef(1) ...  Vcoef(0) ...   a0
c
      Include 'dvmax.h'
      DIMENSION ddcoef(0:Ncoef),dvcoef(0:Ncoef)
      Logical boundary,decrease
      DOUBLE PRECISION NumbTot


c      a0     =2.7e-6  ! initial guess
c      dcoef(0) =3.5
c      vcoef(0) =120.
     
c      do i=1,Nterms   ! initial guess
c          dcoef(i) =0.
c          vcoef(i) =0.
c      EndDo 
c      dcoef(1) =-1.25
c      vcoef(1) =-10.

      da0    =2.e-6        ! initial step for interlopers

      CALL Likehood

      If(Debug)write(*,'(a,g16.8,a,12g13.4)')' Like=',aLike,
     &    ' coeff:',(dcoef(i),i=0,Nterms),
     &             (vcoef(i),i=0,Nterms),a0
      do i=0,Nterms   ! initial steps for other coefficients
          ddcoef(i) =0.05
          dvcoef(i) =0.2
      EndDo 
      pmin(iSh2)  = 0.         ! set limits for initial search of interlopers
      pmax(iSh2)  = da0*10.
      pstep(iSh2) = da0
      do i=0,Nterms        ! set limits for initial search of other parameters
         pmin(i+1)    = dcoef(i)-20.* ddcoef(i)
         pmax(i+1)    = dcoef(i)+ 20.* ddcoef(i)
         pstep(i+1)   =       ddcoef(i)
         pmin(i+iSh1) = vcoef(i)-20.* dvcoef(i)
         pmax(i+iSh1) = vcoef(i)+ 20.* dvcoef(i)
         pstep(i+iSh1)=       dvcoef(i)
      EndDo

      dd_change =1.
      Do iter =1,nItermax
c         write (*,*)  '<===== Iteration ',iter,' ======>'
         decrease = .true.
      Do indx=2,Npars        ! loop over all parameters
         a = pmin(indx)  ! find maximum for current steps pstep
         b = pmax(indx)
         dd= pstep(indx)
         CALL FindMaxOne(indx,a,b,dd,alike_min,a_min,boundary)
         If(boundary)decrease =.false.
c         write (*,'(a,g14.5,a,i4,5g12.4)')' L=', alike_min,
c     &        ' parameter=',indx,a_min,pmin(indx),pmax(indx),pstep(indx)

      EndDo
      CALL Likehood
      If(Debug)write(*,'(a,g16.8,a,14g13.4)')' Like=',aLike,
     &    ' coeff:',(dcoef(i),i=0,Nterms),
     &             (vcoef(i),i=0,Nterms),a0
      dd_change =1.     ! deside to change steps or not
      if(decrease)dd_change =(sqrt(2.))
      Do indx=2,Npars        ! change limits for the search
         a = pmin(indx)
         b = pmax(indx)
         width = (b-a)/dd_change
         CALL GetValue(a_min,indx)
         pmin(indx)= a_min-width/2.
         pmax(indx)= a_min+width/2.
         pstep(indx)=pstep(indx)/dd_change
      EndDo
      
      EndDo                    ! end search of maximum

      Npredict =NumbTot()      ! restore total number of objects
      write(*,*) Npredict,Nsample
      Do i=0,Nterms
         dcoef(i) =dcoef(i)*(Nsample/float(Npredict))
      EndDo 
      a0 = a0*(Nsample/float(Npredict))
      If(Debug)write(*,'(a,g16.8,a,14g13.4)')' Like=',aLike,
     &    ' coeff:',(dcoef(i),i=0,Nterms),
     &             (vcoef(i),i=0,Nterms),a0
       Do i=1,20      ! print new solution
         R = (i/20.)*Rmax
         V = sigma(R)
         D =  dens(R)
         If(Debug)write (*,'(a,f8.2,3(a,g11.3))')' Rad=',R,' dens=',
     &                  D,' velocity=',V
         Rprof(i) =Rprof(i)+R
         Vprof(i) =Vprof(i)+V
         Dprof(i) =Dprof(i)+D
         SigV(i)  =SigV(i) +V**2
         SigD(i)  =SigD(i) +D**2
      EndDo 
      return
      end

c--------------------------------------------------------------------
c               Assign value to parameter indx    
      SUBROUTINE AssignValue(a,indx)
c            indx = 1-(Nterms+1) -> dcoef
c                 = Nterms+2     -> vcoef
c                 = 2(Nterms+1)+1-> interlopers
c
      Include 'dvmax.h'

      Select Case(indx)
      CASE (1:Nterms+1)
         dcoef(indx-1) =a
      CASE(iSh1:Npars-1)
         vcoef(indx-iSh1) =a
      CASE(Npars)
         a0 =a
      CASE DEFAULT
         write (*,*)  ' wrong selection of parameter'
         stop
      End Select
      Return
      End
c--------------------------------------------------------------------
c               Retreive value of parameter indx    
      SUBROUTINE GetValue(a,indx)
c
      Include 'dvmax.h'

      Select Case(indx)
      CASE (1:Nterms+1)
         a =dcoef(indx-1)
      CASE(iSh1:Npars-1)
         a= vcoef(indx-iSh1)
      CASE(Npars)
         a= a0
      CASE DEFAULT
         write (*,*)  ' wrong selection of parameter'
         stop
      End Select
      Return
      End
c--------------------------------------------------------------------
c            find maximum of likelihood function: 
c            change only one argument given by indx
c            dd = step
      SUBROUTINE FindMaxOne(indx,a,b,dd,alike_min,a_min,boundary)
c
      Include 'dvmax.h'
      Logical boundary
      
c      write (*,*)  ' inside MaxOne:',indx,a,b
      iter_max = 2  ! number of iterations
                    ! each iteration decreases dd by dd_change
      dd_change = sqrt(2.)
      iter  = 1
 10   s     = -1.d20
      u     = a
      width = (b-a)
      nn    = INT(width/dd+0.5) ! adjust dd
      dd    = width/nn          ! so that b is the last point
 20   CALL AssignValue(u,indx)  ! find min of aLike for step dd
      CALL Likehood
         if(aLike.gt.s)Then
            s    = aLike
            umin = u
         EndIf 
         u =u +dd
         If(u.lt.b)goto 20
         iter = iter +1          ! change dd if not on boundary
      If(umin.gt.a.and.umin.lt.b)Then  ! inside the boundaries
         a = umin-width/dd_change/2.
         b = umin+width/dd_change/2.
         dd= dd/dd_change
         boundary = .false.
      else 
         a = umin-width/2.  ! On boundary: do not change step
         b = umin+width/2.  ! re-center the search box
         boundary = .true.
      EndIf 
c         write (*,'(3x,a,g14.6,i4,g13.4,a,3g12.3)')
c     &                    ' FunctMin=',s,iter,umin,
c     &                    ' limits: ab,dd=',a,b,dd   
      If(iter.lt.iter_max)goto 10

      alike_min = s
      a_min     = umin
      a         = umin-width/2.  ! re-center the search box
      b         = umin+width/2.
      CALL AssignValue(umin,indx)



      return
      end

c--------------------------------------------------------------------
c                    Likelihood function
      SUBROUTINE Likehood
c
      Include 'dvmax.h'
      DOUBLE PRECISION NumbTot

c      write (*,*)  ' Find Likelehood function'
c      write (*,*)  ' Nsample =',Nsample

      sum = 0.
      Do i=1,Nsample
         R = Rad(i)
         V = abs(Vp(i))
         s = sigma(R)
         dd= dens(R)
         If(s.le.0..or.dd.le.0..or.a0.lt.0.)Then ! give large penalty for
            w = -1.e3                ! unphysical density or sigma
            s = abs(s)
            dd= abs(dd)
         else
            w =0.
         endif
         sum =sum +LOG10(gq2*dd/s/exp(0.5*(V/s)**2)+
     &                a0*R)+w
c         if(w.gt.1.)
c     &    write (*,'(i5,4g11.3,a,4g12.3)')i,R,V,s,sum,
c     &         ' dens,exp=',dens(R),0.5*(V/s)**2,w
      EndDo 
      gg = NumbTot()
      aLike = (sum -Nsample*LOG10(gg))/Nsample+5.2
c      write (*,'(10x,a,4g13.5,a,8g11.3)')
c     &    ' L=',aLike,sum,-Nsample*LOG(gg),gg,' params:',a0,
c     &    (dcoef(i),i=0,1),(vcoef(i),i=0,1)
      return
      end

c--------------------------------------------------------------------
c                   Normalization: total number of objects
c                   for given set of parameters
      DOUBLE PRECISION FUNCTION NumbTot()
c                      
      Include 'dvmax.h'
      External FUNC
      ERR =1.d-6

      CALL DGAUS8 (FUNC, Rmin, Rmax, ERR, ANS, IERR)
      NumbTot = a0*Vmax*(Rmax**2-Rmin**2)/2. +ANS
      return
      end
c--------------------------------------------------------------------
c                    Number-density of true satellites
      DOUBLE PRECISION FUNCTION dens(R)
c
      Include 'dvmax.h'
      s =0.
      Do i=0,Ncoef
         s =s +dcoef(i)*T((R-Rmin)/(Rmax-Rmin),i)
      EndDo
      dens = s
      return
      end
c--------------------------------------------------------------------
c                   Velocity dispersion as function of radius
      DOUBLE PRECISION FUNCTION sigma(R)
c
      Include 'dvmax.h'
      s =0.
      Do i=0,Ncoef
         s =s +vcoef(i)*T((R-Rmin)/(Rmax-Rmin),i)
      EndDo
      sigma = s
      return
      end
c--------------------------------------------------------------------
c                    Integrant for the number of true satellites
      DOUBLE PRECISION FUNCTION FUNC(R)
c
      Include 'dvmax.h'
         FUNC = dens(R)*erf(Vmax/sq2/sigma(R))
      return
      end
c--------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION T(x,n)
      IMPLICIT DOUBLE PRECISION(a-h,o-z)

         If(n.eq.0)Then
            T = 1.
         Else If(n.eq.1)Then
            T = 2.*x-1.
         Else If(n.eq.2)Then
            T = 8.*x**2-8.*x+1.
         Else If(n.eq.3)Then
            T = 32*x**3-48.*x**2+18.*x-1.
         Else If(n.eq.4)Then
            T = 128.*x**4-256.*x**3+160.*x**2-32.*x+1.
         EndIf 

      return
      end

c-------------------------------------------------------------------- 
      SUBROUTINE DGAUS8 (FUN, A, B, ERR, ANS, IERR)
C***BEGIN PROLOGUE  DGAUS8
C***PURPOSE  Integrate a real function of one variable over a finite
C            interval using an adaptive 8-point Legendre-Gauss
C            algorithm.  Intended primarily for high accuracy
C            integration or integration of smooth functions.
C***LIBRARY   SLATEC
C***CATEGORY  H2A1A1
C***TYPE      DOUBLE PRECISION (GAUS8-S, DGAUS8-D)
C***KEYWORDS  ADAPTIVE QUADRATURE, AUTOMATIC INTEGRATOR,
C             GAUSS QUADRATURE, NUMERICAL INTEGRATION
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract  *** a DOUBLE PRECISION routine ***
C        DGAUS8 integrates real functions of one variable over finite
C        intervals using an adaptive 8-point Legendre-Gauss algorithm.
C        DGAUS8 is intended primarily for high accuracy integration
C        or integration of smooth functions.
C
C        The maximum number of significant digits obtainable in ANS
C        is the smaller of 18 and the number of digits carried in
C        double precision arithmetic.
C
C     Description of Arguments
C
C        Input--* FUN, A, B, ERR are DOUBLE PRECISION *
C        FUN - name of external function to be integrated.  This name
C              must be in an EXTERNAL statement in the calling program.
C              FUN must be a DOUBLE PRECISION function of one DOUBLE
C              PRECISION argument.  The value of the argument to FUN
C              is the variable of integration which ranges from A to B.
C        A   - lower limit of integration
C        B   - upper limit of integration (may be less than A)
C        ERR - is a requested pseudorelative error tolerance.  Normally
C              pick a value of ABS(ERR) so that DTOL .LT. ABS(ERR) .LE.
C              1.0D-3 where DTOL is the larger of 1.0D-18 and the
C              double precision unit roundoff D1MACH(4).  ANS will
C              normally have no more error than ABS(ERR) times the
C              integral of the absolute value of FUN(X).  Usually,
C              smaller values of ERR yield more accuracy and require
C              more function evaluations.
C
C              A negative value for ERR causes an estimate of the
C              absolute error in ANS to be returned in ERR.  Note that
C              ERR must be a variable (not a constant) in this case.
C              Note also that the user must reset the value of ERR
C              before making any more calls that use the variable ERR.
C
C        Output--* ERR,ANS are double precision *
C        ERR - will be an estimate of the absolute error in ANS if the
C              input value of ERR was negative.  (ERR is unchanged if
C              the input value of ERR was non-negative.)  The estimated
C              error is solely for information to the user and should
C              not be used as a correction to the computed integral.
C        ANS - computed value of integral
C        IERR- a status code
C            --Normal codes
C               1 ANS most likely meets requested error tolerance,
C                 or A=B.
C              -1 A and B are too nearly equal to allow normal
C                 integration.  ANS is set to zero.
C            --Abnormal code
C               2 ANS probably does not meet requested error tolerance.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, I1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   810223  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  DGAUS8
      INTEGER IERR, K, KML, KMX, L, LMN, LMX, LR, MXL, NBITS,
     1 NIB, NLMN, NLMX
      INTEGER I1MACH
      DOUBLE PRECISION A,AA,AE,ANIB,ANS,AREA,B,C,CE,EE,EF,
     1 EPS, ERR, EST, GL, GLR, GR, HH, SQ2, TOL, VL, VR, W1, W2, W3,
     2 W4, X1, X2, X3, X4, X, H,
     3 D1MACH5, D1MACH4
      DOUBLE PRECISION D1MACH, G8, FUN
      DIMENSION AA(60), HH(60), LR(60), VL(60), GR(60)
      SAVE X1, X2, X3, X4, W1, W2, W3, W4, SQ2,
     1 NLMN, KMX, KML
      DATA X1, X2, X3, X4/
     1     1.83434642495649805D-01,     5.25532409916328986D-01,
     2     7.96666477413626740D-01,     9.60289856497536232D-01/
      DATA W1, W2, W3, W4/
     1     3.62683783378361983D-01,     3.13706645877887287D-01,
     2     2.22381034453374471D-01,     1.01228536290376259D-01/
      DATA SQ2/1.41421356D0/
      DATA NLMN/1/,KMX/5000/,KML/6/
      G8(X,H)=H*((W1*(FUN(X-X1*H) + FUN(X+X1*H))
     1           +W2*(FUN(X-X2*H) + FUN(X+X2*H)))
     2          +(W3*(FUN(X-X3*H) + FUN(X+X3*H))
     3           +W4*(FUN(X-X4*H) + FUN(X+X4*H))))
C***FIRST EXECUTABLE STATEMENT  DGAUS8
C
C     Initialize
C
      D1MACH4 = 1.d-17 
c      K = I1MACH(14)
c      ANIB = D1MACH(5)*K/0.30102000D0
c      NBITS = ANIB
c      NLMX = MIN(60,(NBITS*5)/8)
      NBITS = 16 
      NLMX = MIN(60,(NBITS*5)/8)
      ANS = 0.0D0
      IERR = 1
      CE = 0.0D0
      IF (A .EQ. B) GO TO 140
      LMX = NLMX
      LMN = NLMN
      IF (B .EQ. 0.0D0) GO TO 10
      IF (SIGN(1.0D0,B)*A .LE. 0.0D0) GO TO 10
      C = ABS(1.0D0-A/B)
      IF (C .GT. 0.1D0) GO TO 10
      IF (C .LE. 0.0D0) GO TO 140
      ANIB = 0.5D0 - LOG(C)/0.69314718D0
      NIB = ANIB
      LMX = MIN(NLMX,NBITS-NIB-7)
      IF (LMX .LT. 1) GO TO 130
      LMN = MIN(LMN,LMX)
   10 TOL = MAX(ABS(ERR),2.0D0**(5-NBITS))/2.0D0
      IF (ERR .EQ. 0.0D0) TOL = SQRT(D1MACH4)
      EPS = TOL
      HH(1) = (B-A)/4.0D0
      AA(1) = A
      LR(1) = 1
      L = 1
      EST = G8(AA(L)+2.0D0*HH(L),2.0D0*HH(L))
      K = 8
      AREA = ABS(EST)
      EF = 0.5D0
      MXL = 0
C
C     Compute refined estimates, estimate the error, etc.
C
   20 GL = G8(AA(L)+HH(L),HH(L))
      GR(L) = G8(AA(L)+3.0D0*HH(L),HH(L))
      K = K + 16
      AREA = AREA + (ABS(GL)+ABS(GR(L))-ABS(EST))
C     IF (L .LT .LMN) GO TO 11
      GLR = GL + GR(L)
      EE = ABS(EST-GLR)*EF
      AE = MAX(EPS*AREA,TOL*ABS(GLR))
      IF (EE-AE) 40, 40, 50
   30 MXL = 1
   40 CE = CE + (EST-GLR)
      IF (LR(L)) 60, 60, 80
C
C     Consider the left half of this level
C
   50 IF (K .GT. KMX) LMX = KML
      IF (L .GE. LMX) GO TO 30
      L = L + 1
      EPS = EPS*0.5D0
      EF = EF/SQ2
      HH(L) = HH(L-1)*0.5D0
      LR(L) = -1
      AA(L) = AA(L-1)
      EST = GL
      GO TO 20
C
C     Proceed to right half at this level
C
   60 VL(L) = GLR
   70 EST = GR(L-1)
      LR(L) = 1
      AA(L) = AA(L) + 4.0D0*HH(L)
      GO TO 20
C
C     Return one level
C
   80 VR = GLR
   90 IF (L .LE. 1) GO TO 120
      L = L - 1
      EPS = EPS*2.0D0
      EF = EF*SQ2
      IF (LR(L)) 100, 100, 110
  100 VL(L) = VL(L+1) + VR
      GO TO 70
  110 VR = VL(L+1) + VR
      GO TO 90
C
C     Exit
C
  120 ANS = VR
      IF ((MXL.EQ.0) .OR. (ABS(CE).LE.2.0D0*TOL*AREA)) GO TO 140
      IERR = 2
c      write (*,*)  'DGAUS8: ',
c     +   'ANS is probably insufficiently accurate.'
      GO TO 140
  130 IERR = -1
      write (*,*)   'DGAUS8: ',
     +   'A and B are too nearly equal to allow normal integration.'
  140 IF (ERR .LT. 0.0D0) ERR = CE
      RETURN
      END
C------------------------------------------------
C				                                       random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4
     +	/453815927,2147483648.,2147483647,536870912,131072,256,
     +	 16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END
C--------------------------------------
C				normal random numbers
      DOUBLE PRECISION FUNCTION GAUSS(M)
C--------------------------------------
      DOUBLE PRECISION X,X2
      X=0.
      DO  I=1,5
         X=X+RANDd(M)
      EndDo
      X2   =1.5491933*(X-2.5)
      GAUSS=X2*(1.-0.01*(3.-X2*X2))
      RETURN
      END