c-----------------------------------------------------------------------
c LS2loop.F: last modified: edley -> 2/april/96
c-----------------------------------------------------------------------
c
c  WARNING: All the variables must be initialized to 0.  Make sure
c           that you compile the code with an option that forces that.
c           (e.g., -K in the HP-AUX compiler.)
c
c This is LS2.F by  E.Ley, and M.F.J.Steel.  It is used in 
c "On the Estimation of Demand Systems Through Consumption Efficiency"
c by Eduardo Ley & Mark F.J. Steel
c
c http://econwpa.wustl.edu/eprints/em/papers/9503/
c
c-----------------------------------------------------------------------
c
c  - INPUT FILES:
c          		fdata............xxxxxx.dat
c          		fpara............consume.par
c  - OUTPUT FILES:
c			fout.............consumeX.out
c			ftab.............consumeX.tab
c			fmma.............consumeX.mma
c
c  The execution is driven by fpara.  This file must contain:
c 
c   1. A character which will be appended to 'consume' in the .out and 
c      .err files
c   2. A seed for the random number generator
c   3. Number of repetitions
c     ... a bunch of other stuff ---see read(fpara,.. below
c
c  The strategy is to have many par files: A.par, B.par...
c  (each of them with letter A, B, etc, in the 1st position, 1st line)
c  and rename 'X.par' to 'consume.par' when, say, model X 
c  wants to be run.  At runtime, the parameters will always be read
c  from consume.par.  consumeX.out and consumeX.err will then contain the 
c  output.
c  
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      program consume
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
      real seg, secnds, tottime
c      
      character*1 letra
      character*7 name(4)
c
      double precision logIf0,logIf,logl,logprior,logprior0,mu      
c
      integer fdata,fpara,fout,ftab,fmma
      common/files/fdata,fpara,fout,ftab,fmma
      common/seed/idum
      common/data/xm(41), p(41,3), s(41,3)
c     
      parameter (nbins=50, nobs=41)
      parameter (wscale=1.0d0)
c
c  data
c
      dimension x(nobs,3), y(nobs,3)
c
c  beta
c      
      double precision  dummy(nbins), cdf(4,nbins), pdf(nbins)
      double precision  dmin(4),dmax(4),below(4),above(4),ave(4),ssq(4)
      double precision  beta(4),betam(4),betasd(4),alpha(3),c(3)
      double precision  cprior(3),csd(3),cmean(3),csdp(3),cmeanp(3)
      double precision  d(3)
c
c  shares
c
      double precision ss(3), ssp(3),esharem(3),esharesd(3)
c
c  efficiency
c
      dimension cdfrf(nbins),pdfrf(nbins),rs(6)
c
c-----------------------------------------------------------------------
c .par file
c-----------------------------------------------------------------------
c
       open(fpara,file='consume.par')
c
c-----------------------------------------------------------------------
c read parameters:
c-----------------------------------------------------------------------
c
      read(fpara,9) letra
    9 format(a1) 
      read(fpara,*) idum
      read(fpara,*) numrep
c
c sampling stuff
c
      read(fpara,*) (c(i),i=1,3)
c
c  min and max in post plots: beta: alphas and mu
c
      read(fpara,*) dmin(1),dmax(1)
      read(fpara,*) dmin(2),dmax(2)
      read(fpara,*) dmin(3),dmax(3)
      read(fpara,*) dmin(4),dmax(4)
c
      name(1) = 'alpha1'
      name(2) = 'alpha2'
      name(3) = 'alpha3'
      name(4) = 'mu'
c
c  min and max in post plots: efficiency
c
      read(fpara,*) rfmin,rfmax
c
c-----------------------------------------------------------------------
c  output files
c-----------------------------------------------------------------------
c
      open(fout,file='consume'//letra//'.out')
      open(ftab,file='consume'//letra//'.tab')
      open(fmma,file='consume'//letra//'.mma')
c       
c-----------------------------------------------------------------------
c  read  data 
c-----------------------------------------------------------------------
c
      open(fdata,file='consume.dat')
c
      do 40 i =1,nobs
        read(fdata,*) p(i,1), p(i,2), p(i,3), x(i,1), x(i,2), x(i,3)
        do 35 j=1,3
 35     y(i,j) = p(i,j)*x(i,j)
        xm(i) = y(i,1)+y(i,2)+y(i,3)
        do 37 j=1,3
 37     s(i,j) = y(i,j)/xm(i)
 40   continue
c
c  close the .dat file.
c
      close(fdata)
c
c-----------------------------------------------------------------------
c
c set up values for r^* ---prior efficiency
c
      rs(1) = 0.98d0
      rs(2) = 0.95d0
      rs(3) = 0.90d0
      rs(4) = 0.80d0
      rs(5) = 0.70d0
      rs(6) = 0.50d0 
c
c set up values for Dirichlet Hyperparameters
c
      d(1) = 1.0d0
      d(2) = 3.0d0
      d(3) = 30.0d0
c
      index = 0
      tottime = 0.0
      tottime = secnds(tottime)
c
      write(*,*)
      write(*,*)
      write(*,*) 'Monte Carlo Integration in process...'
      write(*,*) numrep, ' drawings per model'
      write(*,*)
c
c LOOPs
c
      do 5000 jloop=1,3
c
c loop over values of the Dirichlet Hyperparameters
c 
      cprior(1) = d(jloop)
      cprior(2) = cprior(1)
      cprior(3) = cprior(1)      
c
c compute means and std
c
      csump = cprior(1)+cprior(2)+cprior(3)
      do 55 i=1,3
      cmeanp(i) = cprior(i)/csump
 55   csdp(i) = sqrt(cprior(i)*(csump-cprior(i))/((csump**2)*(csump+1)))
c
c
c
      do 5000 kloop=1,6
      index = index+1
c      if (kloop.eq.6) then
c       nmontecarlo = 10*numrep
c      else
c       nmontecarlo = numrep
c      endif
c
c write loop status to std out
c
      write(fout,*)
      write(*,*) index,':',jloop,kloop,'(',18-index, ' models to run)'
c
c loop over values of rstar ---prior efficiency
c
       rstar = rs(kloop)
c
c modify importance sampling stuff for r* = 0.5
c
       if (kloop.eq.6) then
         amu = 25.0d0
         bmu = 0.002d0
       else
         amu = 16.0d0
	 bmu = 0.00125d0	
       endif
c
c initialize stuff
c                  
      do 89 i=1,4
      below(i) = 0.0d0
      above(i) = 0.0d0
      ave(i) = 0.0d0
      ssq(i) = 0.0d0
      betam(i) = 0.0d0
 89   betasd(i) = 0.0d0

      do 90 i=1,4
      do 90 j=1,nbins
      cdf(i,j) = 0.0d0
 90   cdfrf(j) = 0.0d0

      do 91 i=1,3
      esharem(i) = 0.0d0
 91   esharesd(i) = 0.0d0

      brf = 0.0d0
      arf = 0.0d0
      averf = 0.0d0
      ssqrf  = 0.0d0    

      rfm  = 0.0d0
      rfsd = 0.0d0

      wmean = 0.0d0
      wss = 0.0d0
      sw = 0.0d0

c
c-----------------------------------------------------------------------
c  write misc input info out to consume.out
c-----------------------------------------------------------------------
c
      call barra('*',75,fout)
      call barra('*',75,fout)
c
      write(fout, *) 'MODEL ',index
      write(fout, *)
      write(fout, *) 'Random seed and Number of Replications:'
      write(fout, *)      
      write(fout, 100) 'Random seed........', idum
      write(fout, 100) 'Number of draws....', nmontecarlo
c
      call barra('_',75,fout)
c
      csum = c(1)+c(2)+c(3)
      do 95 i=1,3
      cmean(i) = c(i)/csum
 95   csd(i) = sqrt(c(i)*(csum-c(i))/((csum**2)*(csum+1)))
c
      write(fout, *)
      write(fout, *) 'Importance Sampling Stuff:'
      write(fout, *)      
      write(fout, 110) 'Gamma: a ...........', amu
      write(fout, 115) 'Gamma: b ...........', bmu
      write(fout, 115) '   a*b: ............', amu*bmu
      write(fout, 115) '   sqrt(a)*b: ......', sqrt(amu)*bmu
      write(fout, 120) 'Dirichlet: .........',(c(i),i=1,3)
      write(fout, 120) '   Mean: ...........',(cmean(i),i=1,3)
      write(fout, 120) '   SD: .............',(csd(i),i=1,3)
c
      call barra('_',75,fout)
c
      write(fout, *)
      write(fout, *) 'Prior Stuff:'
      write(fout, *)      
      write(fout, 110) 'Prior Efficiency ...', rstar
      write(fout, 120) 'Dirichlet: .........',(cprior(i),i=1,3)
      write(fout, 120) '   Mean: ...........',(cmeanp(i),i=1,3)
      write(fout, 120) '   SD: .............',(csdp(i),i=1,3)
c
      call barra('_',75,fout)
  
  100 format(3x,a20,i12)
  110 format(3x,a20,f12.4)
  115 format(3x,a20,e12.4)
  120 format(3x,a20,3f12.4)
c
c-----------------------------------------------------------------------
      write(fout, *)
      write(fout, *) 'Other Stuff:'
      write(fout, *)      
c
c   Compute pieces of the log importance fcn that don't vary across 
c     draws
c
      logIf0 = -amu*log(bmu)-gammln(amu)
     &      + gammln(c(1)+c(2)+c(3)) 
     &      - gammln(c(1)) - gammln(c(2)) - gammln(c(3))
c
      write(fout, 255) exp(logIf0)
  255 format(3x,'Constant in I_f ........................',e20.4)
c      
      logprior0 = log(-log(rstar)) 
     &  + gammln(cprior(1)+cprior(2)+cprior(3))
     &  - gammln(cprior(1))-gammln(cprior(2))-gammln(cprior(3))
c
       write(fout, 265) logprior0
  265 format(3x,'Log Constant in Prior ..................',e20.4)
      write(fout, *)
c
c-----------------------------------------------------------------------
c
      call barra('*',75,fout)
      call barra('*',75,fout)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c ---------------   start the big monte carlo loop ---------------------
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      seg = 0.0
      seg = secnds(seg)
c
c Starting the Big MonteCarlo Loop'
c
c
      count = 0.0d0
      do 2000 iloop=1, nmontecarlo
c
c  1. Take a draw: (mu,alpha(3))
c
      call draw(mu,amu,bmu,alpha,c,logIf,idum)
c
c  2. Evaluate log posterior at the draw, store it in f
c
      call loglike(mu,alpha,logl)
c
c  3. Calculate weight functions for monte carlo
c
      w =  exp(logl - logIf + logprior(rstar,cprior,mu,alpha))
      if (w.gt.0) then 
        count = count + 1.0d0
      else
        goto 2000
      endif
c            
c
c-----------------------------------------------------------------------
c
c calculate means and standard deviations of the parameters
c & accumulate weights in array cdf for posterior plots
c
c
       beta(1) = alpha(1)
       beta(2) = alpha(2)
       beta(3) = alpha(3)      
       beta(4) = mu
c
c  Loop over elements in beta
c
      do 1300 i=1,4
c
c moments --- they're also computed by hist; 
c we leave these redundant statements for checking purposes
c NOTE: beta is (alpha,mu) ---ie, beta(4) = mu   
c
       betam(i)=betam(i) + w*beta(i) 
       betasd(i)=betasd(i) + w*beta(i)**2
c
c cdf's  of beta    
c      
      do 1100 j=1, nbins
 1100 dummy(j) = cdf(i,j) 
      call acc(beta(i),w, dmin(i), dmax(i), nbins, dummy,
     & below(i), above(i), ave(i), ssq(i))     
      do 1200 j=1, nbins
 1200 cdf(i,j) = dummy(j) 
 1300 continue
c
c-----------------------------------------------------------------------
c
c  Shares: E[share_i]
c
       den = 1.0d0+3.0d0*mu
       do 1450 i=1,3
 1450  esharem(i) = esharem(i) + w*(alpha(i) + mu)/den
       den = (alpha(1)+alpha(2)+alpha(3)+3.0d0*mu)/mu
       den = den*(den+1.0d0)
       do 1500 i=1,3
 1500  esharesd(i) = esharesd(i) 
     &  + w*((alpha(i)/mu+1.0d0)**2+(alpha(i)/mu+1.0d0))/den
c
c-----------------------------------------------------------------------
c
c  Efficiency
c
c  1. For each draw of alpha and mu, draw shares from:
c
c     Dirichlet(alpha_1/mu +1,alpha_2/mu +1,alpha_3/mu +1)
c
      do 1600 i=1,3
 1600 ssp(i) = alpha(i)/mu + 1.0d0
c      
      call dir3dev(ss,ssp,idum)
c
c  2. Construct Efficiency:
c
c     r_f = psi/mu =  (s_1/alpha_1)^alpha_1
c                   x (s_2/alpha_2)^alpha_2
c                   x (s_3/alpha_2)^alpha_3
c
      rf = 1.0d0
      do 1700 i=1,3
1700  rf = rf * (ss(i)/alpha(i))**alpha(i)
c
c  3. Accumulate stuff
c
       rfm   = rfm  + w*rf
       rfsd  = rfsd + w*rf**2
c       
       call acc(rf,w,rfmin,rfmax,nbins,cdfrf,brf,arf,averf,ssqrf)
c
c-----------------------------------------------------------------------
c
c  accumulate to compute mean and se of weights
c
      sw = sw + w
      w = w*wscale
      diloop = dble(iloop)
      wss = wss + ((diloop-1.0d0)/diloop)*((w - wmean)**2)
      wmean = wmean + (w - wmean)/diloop
c
c-----------------------------------------------------------------------
c
c
 2000 continue
c
c      
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c ----------------   end of biggg monte carlo looop  -------------------
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      write(fout,*)
      write(fout,2010) (count/dble(nmontecarlo))*100.0d0
 2010 format(2x,'Weights greater than zero: ',f6.2,' %')
c
      write(fout,*)
      write(fout,*) '  POSTERIOR RESULTS:'
      write(fout,*)
c
c  posterior plots: beta
c
      do 2020 i=1,4
      do 2015 j=1,nbins
 2015 dummy(j) = cdf(i,j)
      call barra('-',75,fout)
 2020 call hist(name(i), index, dummy, pdf, dmin(i), dmax(i), nbins,
     & below(i), above(i), ave(i), ssq(i))
c
c-----------------------------------------------------------------------
c
c  Efficiency stuff
c
c  p(e_f|data)
c
      call barra('-',70,fout)
      call hist('Efficiency',index, cdfrf, pdfrf, rfmin, rfmax, 
     & nbins, brf, arf, averf, ssqrf)
c
c-----------------------------------------------------------------------
c 
c-----------------------------------------------------------------------
c
      do 2030 i=1,4
      betam(i) = betam(i)/sw
 2030 betasd(i)=sqrt((betasd(i)/sw) - betam(i)**2)
c
c   Efficiency
c
      rfm   = rfm/sw
      rfsd  = sqrt((rfsd/sw) - rfm**2)
c
c  Expected Shares
c
      do 2040 i=1,3
      esharem(i)   = esharem(i)/sw
 2040 esharesd(i)  = sqrt((esharesd(i)/sw) - esharem(i)**2)
c
      call barra('-',75,ftab) 
      write(ftab,*) index
      write(ftab,2063) rstar,cprior(1),cprior(2),cprior(3)
      write(ftab,*)
      call barra('-',75,fout) 
      write(fout,*)
      write(fout,*) index,' >> Posterior moments:'
      write(fout,*)
      write(fout,*) '  ------------------------------------' 
      write(fout,*) '         vname       MEAN      S.D.' 
      write(fout,*) '  ------------------------------------' 
      do 2050 i=1,4
      write(ftab,2061) name(i), betam(i), betasd(i)
 2050 write(fout,2061) name(i), betam(i), betasd(i)
      write(fout,*)
      write(ftab,*)
      write(fout,2061) 'Efficiency', rfm, rfsd
      write(ftab,2061) 'Efficiency', rfm, rfsd
      do 2055 i=1,3
      write(ftab,2062) 'share',i, esharem(i), esharesd(i)     
 2055 write(fout,2062) 'share',i, esharem(i), esharesd(i)     
 2061 format(2x,a15,2x,f8.4,2x,f8.4)
 2062 format(2x,a14,i1,2x,f8.4,2x,f8.4)
 2063 format(f8.2,4x,'(',f4.1, ',' ,f4.1, ',' ,f4.1, ')' )
      write(fout,*)
c
c-----------------------------------------------------------------------
c
c  weights
c
      write(fout,*)
      call barra('-',75,fout) 
      write(fout,*)
      write(fout,*) '>> WEIGHTS:'
      write(fout,*)
      write(fout,2075)  'Mean', wmean/wscale
      write(fout,2075)  'S.E.', sqrt(wss/diloop)/wscale
      write(fout,2075)  'V.C.', sqrt(wss/diloop)/wmean
 2075 format(2x,a4,' of the weights......', e20.4)     
c
      call barra('.',75,fout) 
c
c-----------------------------------------------------------------------
	write(fout,*)      
	write(fout,*)
	call barra('*',75,fout)
      	call wrtime(fout,seg)
      	call barra('*',75,fout)
       	
c-----------------------------------------------------------------------
        write(fout,*)
 5000   call barra('$',75,fout) 

c-----------------------------------------------------------------------
	write(fout,*)      
	write(fout,*)
	call barra('&',75,fout)
      	call wrtime(fout,tottime)
      	call barra('&',75,fout)
       	
c-----------------------------------------------------------------------
      	
      	
c-----------------------------------------------------------------------
        write(fout,*) '(This file is consume'//letra//'.out)'
c-----------------------------------------------------------------------
      stop
      end
c
c***********************************************************************
c***********************************************************************
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c          end of main program --- subroutines and fcns follow
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c***********************************************************************
c***********************************************************************
c
c-----------------------------------------------------------------------
c  file handling
c-----------------------------------------------------------------------
c
       block data
        common/files/fdata,fpara,fout,ftab,fmma
        integer fdata,fpara,fout,ftab,fmma
        data fdata,fpara,fout,ftab,fmma/10,11,12,13,14/
       end
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: logprior
c-----------------------------------------------------------------------
c
       double precision function logprior(rstar,cprior,mu,c)
       implicit double precision (a-h,o-z)
       double precision mu,rstar,cprior(3),c(3)
c
       logprior = -2.0d0*log(mu)+log(rstar)/mu
     &  +(cprior(1)-1.0d0)*log(c(1)) 
     &  +(cprior(2)-1.0d0)*log(c(2)) 
     &  +(cprior(3)-1.0d0)*log(c(3))
c
      return
      end
c
c-----------------------------------------------------------------------
c          subroutine: loglike
c-----------------------------------------------------------------------
c
      subroutine loglike(mu,a,f)
c
      implicit double precision (a-h,o-z)
      double precision mu, a(3), f, c(3), mui
c      
      common/data/xm(41),p(41,3),s(41,3)
c
c  setup
c
      mui = 1.0d0/mu
      do 50 i=1,3
  50  c(i) = a(i)*mui
c
c   integrating constant
c            
      x = gammln(c(1)+1.0d0)+gammln(c(2)+1.0d0)+gammln(c(3)+1.0d0)
      x = 41*(gammln(3.0d0+mui)-x)
c
c   sum over observations
c
      z = 0.0d0
      do 70 i=1,41
      y = log(p(i,1))+log(p(i,2))-2.0d0*log(xm(i))
  70  z = z + y + c(1)*log(s(i,1))+c(2)*log(s(i,2))+c(3)*log(s(i,3))
c
c  log posterior:
c
      f = x+z
c
      return
      end
c
c-----------------------------------------------------------------------
c

c-----------------------------------------------------------------------
c          function: draw
c-----------------------------------------------------------------------
c
      subroutine draw(mu,amu,bmu,alpha,c,logIf,idum)
c
      implicit double precision (a-h,o-z)
      double precision mu, amu, bmu, alpha(3), c(3), logIf
c
c  1. draw mu from gamma(amu,bmu)
c 
      mu = gam2dev(amu,bmu,idum)
c
c  2. evaluate varying part of log(f_G(mu|amu, bmu))
c
      logIf = (amu-1)*log(mu) - (mu/bmu)
c
c  3. draw alpha from f_D(alpha|c)
c
      call dir3dev(alpha,c,idum)
c
c  4. evaluate varing part of log(f_D(alpha|c)) and sum it to logIf
c  
      logIf =  logIf + (c(1)-1.0d0)*log(alpha(1))
     &      + (c(2)-1.0d0)*log(alpha(2)) + (c(3)-1.0d0)*log(alpha(3))        
c
      return
      end          
c
c
c-----------------------------------------------------------------------
c          function: dir3dev
c-----------------------------------------------------------------------
c
      subroutine dir3dev(x,a,idum)
c
      implicit double precision (a-h,o-z)
c      
      double precision a(3), x(3), gamdev, sum
c
      sum = 0.0d0
      do 10 i=1,3
         x(i) = gamdev(a(i),idum)
 10   sum = sum + x(i)
 
      do 20 i=1,3
 20   x(i) = x(i)/sum
c
      return
      end
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: gam2dev
c-----------------------------------------------------------------------
c
      double precision function gam2dev(a,b,idum)
c
      implicit double precision (a-h,o-z)
      double precision a,b,gamdev          
c
      gam2dev = b*gamdev(a,idum)
      return
      end
c
c-----------------------------------------------------------------------
c          function: gamdev
c-----------------------------------------------------------------------
c
      double precision function gamdev(a,idum)
      implicit double precision (a-h,o-z)
      if (a.le.1) then
         gamdev = gam1(a,idum)
      else
         gamdev = gam2(a,idum)
      endif
      return
      end
c
c-----------------------------------------------------------------------
c          function: gam1
c-----------------------------------------------------------------------
c
      double precision function gam1(a,idum)
c
c  a<1: Stuart's Thm. Devroye p 420
c
      implicit double precision (a-h,o-z)
      x = a+1.0d0
      gam1 = gam2(x,idum)*ran2(idum)**(1.0d0/a)
      return
      end

c
c-----------------------------------------------------------------------
c          function: gam2
c-----------------------------------------------------------------------
c
      double precision function gam2(a,idum)
c
c a>1: Cheng (1977) algorithm.  See Devroye (1986) page 413.
c Note that there is are typos there: x = a*exp(y) not exp(v)!
c and lambda...
c
      implicit double precision (a-h,o-z)
      double precision a,b,c,u,v,x,y,z,r,test1,test2,ran2,lambda
c
c  setup
c
      b = a - log(4.0d0)
      c = a + sqrt(2.0d0*a -1.0d0)
      lambda = 1.0d0/sqrt(2.0d0*a -1.0d0)
c
c  generator
c
 100  continue
      u = ran2(idum)
      v = ran2(idum)
c
      y = (log(v)-log(1.0d0-v))*lambda
      x = a*exp(y)
      z = u*v*v
      r = b + c*y - x
c
c  Here's a quick acceptance step suggested by Cheng
c
      test1 = 4.5d0*z  - 2.5040774d0
      if (r.ge.test1) goto 500
c
c  Here's the real thing...
c
      test2 = log(z)
      if (r.ge.test2) goto 500 
      goto 100
c
 500  gam2 = x
      return
      end
c
c
c-----------------------------------------------------------------------
c          function: gammln
c-----------------------------------------------------------------------
c
c This function computes ln(gamma(x)) for x > 0.
c Taken from Pascal Numerical Recipes p 177
c
c Modified to improve accuracy when x<1 using:
c
c                     pi*e
c Gamma(1-e) = ---------------------
c               Gamma(1+e)*sin(pi*e)
c
      double precision function gammln(xx)
c
      double precision xx
      double precision cof(6),stp,half,one,fpf,x,y,tmp,ser,g,dpi
c
      data cof/76.18009173d0,-86.50532033d0,24.01409822d0, 
     &         -1.231739516d0,0.120858003d-2,-.536382d-5/
      data stp/2.50662827465d0/
      data dpi/3.14159265358979d0/
      data half,one,fpf/0.5d0,1.0d0,5.5d0/
c                 
      if (xx.lt.1) then
      x = one - xx
      tmp =  x + fpf
      tmp = (x + half)*log(tmp)-tmp
      ser = one
      do 1000 j=1,6
         x=x+one
         ser = ser+cof(j)/x
 1000 continue
      g = tmp + log(stp*ser)
      y = dpi*(one-xx)
      gammln = log(y/sin(y)) - g
c
      else
c
      x = xx - one
      tmp =  x + fpf
      tmp = (x + half)*log(tmp)-tmp
      ser = one
      do 1100 j=1,6
         x=x+one
         ser = ser+cof(j)/x
 1100 continue
      gammln = tmp + log(stp*ser)
      endif
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of function gammln
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c     subroutine barra
c-----------------------------------------------------------------------
c
c
      subroutine barra(car,n,fout)
c
c-----------------------------------------------------------------------
c
      integer i, n, fout
      character*1 car
      character*1 bar(80)
      do 1000 i=1,n
 1000 bar(i) = car
      write(fout,*) (bar(i),i=1,n)
      return
      end
c
c-----------------------------------------------------------------------
c   end of subroutine barra
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: ran2
c-----------------------------------------------------------------------
c
      double precision function ran2(idum)
c
c   this function contains a portable random number generator 
c   uniform [0,1]
c   (numerical recipes p. 273)
c
      implicit double precision (a-h,o-z)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      double precision AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1.d0/IM1,IMM1=IM1-1,
     &IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     &NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1.0d0-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END      
c
c-----------------------------------------------------------------------
c          end of function: ran2
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c     subroutine acc
c-----------------------------------------------------------------------
c
      subroutine acc(datum,weight,datmin,datmax,nbins,cdf,b,a,su,sq)
c      
c-----------------------------------------------------------------------
c
      implicit double precision(a-h,o-z)
c      
      double precision datum, weight, datmin, datmax, b, a, su, sq
      integer i, nbins
      double precision cdf(nbins)
c
      double precision x, step
c
      su = su + weight*datum
      sq = sq + weight*datum**2
c
      if (datum.lt.datmin) b = b + weight
      if (datum.ge.datmax) a = a + weight
c
      x = datmin
      step = (datmax-datmin)/nbins
c
      do 1000 i=1,nbins
        x = x + step
        if (datum.le.x) cdf(i)=cdf(i)+weight
 1000 continue
      return
      end 
c
c-----------------------------------------------------------------------
c     end of subroutine acc
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c
      subroutine hist(name,index,cdf,pdf,datmin,datmax,nbins,b,a,su,sq)
c
c-----------------------------------------------------------------------
c
      implicit double precision(a-h,o-z)
c 
      character*(*) name
      integer nbins, fdata,fpara,fout,ftab,fmma, index
      double precision datmin, datmax, b, a, su, sq, prev
      dimension cdf(nbins), pdf(nbins)
c
      common/files/fdata,fpara,fout,ftab,fmma
c      
      double precision max, sum, bound, step, midpoint
      integer i
c
      character*50 height
c 
      height='XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
c      
      step = (datmax-datmin)/nbins
      pdf(1) = cdf(1)
      sum = a + pdf(1)
      do 1000 i=2,nbins
      pdf(i) = cdf(i)-cdf(i-1)
      sum = sum + pdf(i)
 1000 continue
c
      su = su/sum
      sq = sqrt((sq/sum)-su**2)
c
c normalize
c
      max = 0.0d0
      do 2000 i=1,nbins
      pdf(i) = pdf(i)/sum
      if (pdf(i).gt.max) max=pdf(i)
 2000 continue 
      bound = datmin
      write(fout,*)
      write(fout,2060)
      write(fout,2050) name
 2050 format(t15,'%% ',t30, a10, t65,'%%')
       write(fout,2051)
 2051 format(t15, '%% ', t65, '%%')
      write(fout,2052) su,sq
 2052 format(t15,'%%   Mean...',g13.4,4x,'STD...',g10.4,t65, '%%')
      write(fout,2060)
 2060 format(t15,'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
      write(fout,*)
c
      write(fout,2100) '   UBound    PDF   CDF(%)'
      write(fout,2100) '+---------------------------'
 2100 format(t51,a26)
      prev = 0.0d0
      do 3100 i=1,nbins
      bound = bound + step
      if (((cdf(i)/sum).ge.(0.001d0)).and.(prev.lt.(0.999d0))) then
      write(fout,3099) height(1:nint((pdf(i)/max)*50)), bound,
     & pdf(i), (cdf(i)/sum)*100.0d0
 3099 format(a50, '|',f9.3, 1x, f6.4, 1x, f8.4)
      endif
      prev = cdf(i)/sum
 3100 continue
c
c mass to the left and to the right of range bounds
c
      write(fout,*)
      write(fout,4000) datmin, (b/sum)*100.0d0
 4000 format(t15,'Mass to the left of  ',f16.4,' is.....',f8.4,' %')
      write(fout,4010) datmax, (a/sum)*100.0d0
 4010 format(t15,'Mass to the right of ',f16.4,' is.....',f8.4,' %')
c
c print out a table for Mathematica
c
      write(fmma,*) name,index,' ={' 
      midpoint = datmin + step/2.0d0
      do 5100 i=1,nbins-1
      write(fmma,4099) midpoint,pdf(i)
 4099 format('{',f8.4,','f14.10,'},')
 5100 midpoint = midpoint + step
      write(fmma,4100) midpoint,pdf(nbins)
 4100 format('{',f8.4,','f14.10,'}};')
c      
      return
      end
c
c-----------------------------------------------------------------------
c     end of sub hist
c-----------------------------------------------------------------------
c
                real FUNCTION secnds(t)
                secnds = (LONG(362)/60.00)-t
                END
c
c-----------------------------------------------------------------------
c
      subroutine wrtime(filename,t)
c
c-----------------------------------------------------------------------
c
      integer filename, hours, minutes, seconds
      real t
c
c TIME
c
      minutes = 0
      hours = 0
      seconds = int(secnds(t))
      minutes = seconds/60
      hours   = minutes/60
c      
      minutes = minutes - hours*60
      seconds = seconds - minutes*60 - hours*60*60
c 
      write(filename,100) hours, minutes, seconds
c
 100  format('***',t10,'-> Total CPU time: ',i3,' hour(s),',i3,
     & ' minute(s), and ',i3,' seconds.', t73,'***')
c     
      return
      end
c
c-----------------------------------------------------------------------
c

	
    
	

     

---dat file---cut here (incl this line!)

   27.40d0     26.88d0     16.85d0     20.43d0     90.88d0     50.60d0  
   28.77d0     29.18d0     17.77d0     22.85d0     96.60d0     55.48d0  
   29.95d0     27.85d0     18.45d0     25.05d0     94.85d0     58.42d0  
   31.23d0     27.73d0     18.90d0     30.73d0     98.22d0     63.15d0  
   31.85d0     29.73d0     19.40d0     29.85d0    109.15d0     69.03d0  
   31.57d0     29.85d0     20.30d0     29.23d0    114.72d0     75.13d0  
   31.77d0     29.38d0     21.40d0     32.67d0    117.83d0     82.13d0  
   32.10d0     30.02d0     22.18d0     32.10d0    119.67d0     88.05d0  
   33.38d0     30.60d0     23.02d0     38.88d0    124.70d0     94.30d0  
   34.73d0     30.95d0     24.02d0     38.20d0    130.78d0    101.63d0  
   36.50d0     31.90d0     25.02d0     39.65d0    137.10d0    108.55d0  
   36.92d0     32.85d0     26.10d0     37.17d0    141.75d0    115.67d0  
   38.50d0     33.17d0     26.93d0     42.80d0    148.47d0    125.00d0  
   38.77d0     33.63d0     27.80d0     43.42d0    153.20d0    134.00d0  
   38.88d0     34.08d0     28.40d0     41.90d0    157.40d0    141.80d0  
   39.75d0     34.58d0     29.13d0     47.02d0    163.82d0    151.05d0  
   40.15d0     34.90d0     29.85d0     51.80d0    169.35d0    160.63d0  
   40.35d0     35.25d0     30.60d0     56.85d0    179.68d0    172.78d0  
   40.92d0     36.10d0     31.43d0     63.48d0    191.85d0    185.40d0  
   41.55d0     37.35d0     32.60d0     68.53d0    208.45d0    200.30d0  
   42.33d0     38.23d0     33.80d0     70.63d0    216.90d0    216.00d0  
   44.55d0     39.98d0     35.65d0     81.00d0    235.00d0    236.43d0  
   46.10d0     41.98d0     37.70d0     86.22d0    252.18d0    259.43d0  
   47.80d0     44.05d0     40.38d0     85.67d0    270.32d0    284.02d0  
   50.25d0     45.92d0     43.08d0     97.58d0    283.27d0    310.65d0  
   51.05d0     47.90d0     45.58d0    111.22d0    305.10d0    341.27d0  
   50.83d0     53.75d0     48.10d0    124.72d0    339.55d0    372.98d0  
   54.08d0     59.58d0     51.77d0    123.75d0    380.90d0    411.90d0  
   61.15d0     65.13d0     56.38d0    135.35d0    416.20d0    461.23d0  
   65.88d0     67.60d0     60.65d0    161.45d0    451.95d0    515.92d0  
   69.15d0     71.03d0     65.45d0    184.50d0    490.45d0    582.25d0  
   72.78d0     76.00d0     70.30d0    205.57d0    541.80d0    656.10d0  
   78.25d0     83.55d0     75.88d0    218.95d0    613.25d0    734.55d0  
   85.55d0     88.85d0     83.72d0    219.28d0    681.35d0    831.95d0  
   93.65d0     96.78d0     92.30d0    239.88d0    740.58d0    934.70d0  
  100.00d0    100.00d0    100.03d0    252.65d0    771.00d0   1026.97d0  
  101.60d0    102.50d0    106.13d0    289.10d0    816.70d0   1128.75d0  
  102.97d0    106.17d0    111.60d0    335.55d0    867.30d0   1227.63d0  
  102.35d0    108.50d0    117.25d0    368.70d0    913.13d0   1347.52d0  
  101.38d0    110.13d0    122.25d0    402.43d0    939.35d0   1458.05d0  
  100.40d0    114.05d0    127.42d0    413.73d0    982.88d0   1571.22d0


---par file---cut here (incl this line!)
W									letra
-9312773							idum
1000000						    	numrep
30.0d0   90.0d0   80.0d0			c3
0.10d0  0.25d0						alpha1
0.40d0   0.55d0						alpha2
0.30d0	 0.45d0						alpha2
0.010d0  0.080d0  					mumin mumax
0.60d0	1.0d0						efficiency
