program control
      implicit none
      include 'mkcl.inc'


        real paramval(nparam),psparamval(npsparam)
        real lnL

C In mkcl.inc nparam=2, npsparam=2        
C paramval(nparam) is an array containing values for b and d 
C psparamva(nparam is an array containing values for tilt and amplitude.l
C lnL is the log of the likehood of a given set of parameters


        integer count, accept, ns, l, advnst
 
C count is temp variable that is 0 if value is not accepted, 1 if it is.
C accept keeps a running tally of number of accepted points.
C ns is total number of runs to do.
C l is dummy variable in loop from 1 to ns
C advnst is the value to advance the nest for our random number generator.

        
        double precision dlog, u, ub, ud, test, divide, nsd
        double precision epsb, epsd, b0, d0, prev, cur, ratio
        double precision epstilt, epsamp, amp, tilt, ut, ua


        REAL dumb, RANF

C RANF is random number generator (between 0 and 1) from ranlib.f package at
C http://www.netlib.org/random/
C  dumb is a dummy value used to advance the seed in the random number generator    

        data ns/300000/
        open(2, file='results')




C     Eps values

        epsb=0.0005
        epsd=0.005
        epstilt=.01
        epsamp=0.25

C     Initial param values   - psparamval(1) = amp psparamval(2) = tilt

C        psparamval(1) = 95.3
C        psparamval(2) = 1.1

        amp=95.3
        tilt=1.1
        psparamval(1)=amp
        psparamval(2)=tilt


C     Initial psparam values

C        paramval(1)=0.0159
C        paramval(2)=0.16

        b0=0.0159
        d0=0.16
        paramval(1)=b0
        paramval(2)=d0

        ns = 200000

C Advance the random number seed

        dumb = RANF(332.32)
        dumb = advnst(9)


C Call Lloyd Knox's cmb fast program which takes parameter values as inputs and outputs the log 
C of the likelyhood of having those parameter values. 

        call dgetlnlike(paramval,psparamval,lnL)
        prev = lnL
        accept=0
        do 10 l=1,ns

C Create new b param value and check that it is within the range


           u=RANF()
 20        ub=RANF()
        
          paramval(1)=b0+epsb*(2*ub-1)


          if ((paramval(1).lt.0.015d0).OR.(paramval(1).gt.0.025d0)) then
C             print*, "B Out of Bounds"
             GO TO 20
          else 
          continue
          endif

C Create new d param value and check range


 30       ud=RANF()
          paramval(2)=d0+epsd*(2*ud-1)

          if ((paramval(2).lt.0.1d0).OR.(paramval(2).gt.0.2d0)) then
C             print*, "d out of range"
          GO TO 30
          else 
          continue
          endif

C Create new amp and check range

 40       ua=RANF()
          psparamval(1)=amp+epsamp*(2*ua-1)

          if ((psparamval(1).lt.0.0d0).OR.(psparamval(1).gt.200d0)) then
C         print*, "amp out of range"
          GO TO 40
          else
          continue
          endif

C Create new tilt and check range

 50       ut=RANF()
          psparamval(2)=tilt+epstilt*(2*ut-1)
 
          if ((psparamval(2).lt.0.0d0).OR.(psparamval(2).gt.2.0d0)) then
C         print*, "tilt out of range"
          GO TO 50
          else
          continue
          endif


C Get new log(likelyhood) with new parameters (Lloyd Knox's cmb fast program)

           call dgetlnlike(paramval,psparamval,lnL)
           cur = lnL

           test=dlog(u)+prev-cur

C if test is positive, then accept value; make those parameters are base parameters and
C add one to the count of acceptances.

           if (test.lt.0.d0) then

              prev=cur
                b0=paramval(1)
                d0=paramval(2)
                amp=psparamval(1)
                tilt=psparamval(2)
                 count=1
           else
                count=0

           endif

C Format stuff to write the line number and the parameters to our file.

           write(2, 900) l, b0, d0, amp, tilt
 900       format(I7,1X,F8.6,1X,F8.6,1X,F10.6,1X,F8.6)

C Add one to the acceptance count, if count=1
           accept=accept+count


10       continue

C print the ratio (percent) of acceptances and the totaly acceptance number

         divide = accept
         nsd = ns
         ratio=divide/nsd
         print*,ratio
C         print*,accept/ns
         print*,accept

	END