SUBROUTINE runchn(chlen, chnpar, chburn, filenm, chvals, qrand, 
     &                  iear, ifluxar, ibflxar, iigroup, iichanb, 
     &                  iichane, irespon, iesavar, ieffar, nfpar, 
     &                  parasf, parval)

      IMPLICIT NONE

      INTEGER chlen, chnpar, nfpar

      INTEGER chburn
      INTEGER iear, ifluxar, ibflxar, iigroup, iichanb, iichane
      INTEGER irespon, iesavar, ieffar

      DOUBLE PRECISION chvals(chnpar)
      DOUBLE PRECISION parasf(nfpar)
      DOUBLE PRECISION parval(nfpar, 2)

      LOGICAL qrand

      CHARACTER filenm*(*)

c Routine to run an MCMC chain starting at the current parameter values
c Arguments :
c     chlen       i      i: The length of the chain to run
c     chnpar      i      i: The number of parameters (+ 1 for the statistic)
c     chburn      i      i: The burn-in length for the chain
c     filenm      c      i: Output filename
c     chvals      d      w: Workspace array for chain values
c     qrand       l      i: If true then use randomized starting point
c     iear	  i	 i: pointer to upper energy for ith energy range
c     ifluxar	  i      i: pointer to model in energy space
c     ibflxar	  i      i: pointer to background model in energy space
c     iigroup	  i	 i: pointer to number of response groups for ith energy
c     iichanb	  i	 i: pointer to start bin for ith response group
c     iichane	  i	 i: pointer to end bin for ith response group
c     irespon	  i	 i: pointer to response matrix elements
c     iesavar	  i	 i: pointer to save array for model components
c     ieffar	  i	 i: pointer to effective areas
c     nfpar       i      i: number of fit parameters
c     parasf      d      w: work array for saved parameter values
c     parval      d      w: work array for parameter values

      INCLUDE '../../inc/xspec.inc'

      REAL random, oldstat, ftstat

      INTEGER ieftprm, ievalue, ieigenv
      INTEGER ipar, iv, istat, neigen, impar, i, lun
      INTEGER ierr, j

      LOGICAL qincnw, qdone

      CHARACTER wrtstr*2048

      DOUBLE PRECISION fgpval
      REAL xgfstt
      INTEGER fgnfpr, dgdpnt, fgneig, fgnvpr, lenact
      LOGICAL fgqfrz
      EXTERNAL fgpval, xgfstt, fgnfpr, dgdpnt, fgneig, fgqfrz
      EXTERNAL fgnvpr, lenact

c test inputs

      IF ( chlen .LE. 0 )
     &  CALL xwrite('RUNCHN: Invalid chain length', 10)
      IF ( chnpar .LE. 0 )
     &  CALL xwrite('RUNCHN: Invalid number of parameters', 10)
      IF ( chburn .LT. 0 )
     &  CALL xwrite('RUNCHN: Invalid burn-in length', 10)

      IF ( chnpar-1 .NE. fgnvpr() )
     &  CALL xwrite('RUNCHN: Unexpected number of parameters', 10)

      IF ( lenact(filenm) .EQ. 0 )
     &  CALL xwrite('RUNCHN: No filename given', 10)

c save the original parameter values

      nfpar = fgnfpr()
      DO ipar = 1, nfpar
         parasf(ipar) = fgpval(ipar, 'v')
      ENDDO

c get the eigenfunctions and values to use to generate new sets of parameters

      neigen = fgneig()
      ieftprm = -1
      CALL udmget( neigen, 4, ieftprm, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Memory allocation for EFTPRM failed in RUNCHN.', 2)
      ievalue = -1
      CALL udmget( neigen, 7, ievalue, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Memory allocation for EVALUE failed in RUNCHN.', 2)
      ieigenv = -1
      CALL udmget( neigen**2, 7, ieigenv, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Memory allocation for EIGENV failed in RUNCHN.', 2)

      CALL fgeftp(neigen, MEMI(ieftprm), istat)
      CALL fgevec(neigen, MEMD(ieigenv), istat)
      CALL fgeval(neigen, MEMD(ievalue), istat)
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Failed to load eigenvectors in RUNCHN.', 2)

c set the first point in the chain to the current parameter values
c or a randomized starting point

      IF ( qrand ) THEN

c triple the eigenvalues so we get a 3 sigma randomization

         DO i = 1, neigen
            MEMD(ievalue+i-1) = MEMD(ievalue+i-1)*3
         ENDDO

c generate new parameters till we get a valid set

         qdone = .FALSE.
         DO WHILE ( .NOT.qdone )
            CALL simprm(neigen, MEMI(ieftprm), MEMD(ieigenv), 
     &                  MEMD(ievalue), nfpar, parasf, parval(1,1))
            qdone = .TRUE.
            DO ipar = 1, nfpar
               IF ( parval(ipar,1) .LT. fgpval(ipar,'l') .OR.
     &              parval(ipar,1) .GT. fgpval(ipar,'h') ) THEN
                  qdone = .FALSE.
               ENDIF
            ENDDO
         ENDDO

c set the parameters just generated

         DO ipar = 1, nfpar
            CALL fppval(ipar, parval(ipar,1), 'v', istat)
         ENDDO

         CALL evlmod(memr(iear), .FALSE., memr(iesavar), memr(ifluxar), 
     &               memr(ibflxar), memr(ieffar))
         CALL fldmod(memr(ifluxar), memr(ibflxar), memi(iigroup), 
     &               memi(iichanb), memi(iichane), memr(irespon), 
     &               memr(dgdpnt('mod')), memr(ieffar))
         CALL calstat(memr(dgdpnt('mod')), memr(dgdpnt('src')), 
     &                memr(dgdpnt('csv')), memr(dgdpnt('cor')), 
     &                memr(dgdpnt('bkg')), memr(dgdpnt('cbv')), 
     &                memr(dgdpnt('reg')), memr(dgdpnt('brg')), 
     &                memr(dgdpnt('are')), memr(dgdpnt('bre')), 
     &                ftstat)
         chvals(chnpar) = ftstat

c reset the eigenvalues

         DO i = 1, neigen
            MEMD(ievalue+i-1) = MEMD(ievalue+i-1)/3
         ENDDO

      ENDIF

c Set the first elements of the chain. If the randomization was in use this
c will be offset from the fit else it will be the result of the last fit before
c running the chain command.

      iv = 0
      DO ipar = 1, nfpar
         IF ( .NOT. fgqfrz(ipar) ) THEN
            iv = iv + 1
            chvals(iv) = fgpval(ipar, 'v')
         ENDIF
      ENDDO

c If the randomization was not in use we need to set the statistic value

      IF ( .NOT.qrand ) chvals(chnpar) = xgfstt()

c Open the output file and write the first row

      wrtstr = 'Writing chain to '//filenm(:lenact(filenm))
      CALL xwrite(wrtstr, 25)

      CALL getlun(lun)
      CALL openwr(lun, filenm, 'new', ' ', ' ', 0, 1, ierr)
      IF ( ierr .NE. 0 ) THEN
         wrtstr = 'RUNCHN: Failed to open '// 
     &            filenm(:lenact(filenm))
         CALL xwrite(wrtstr, 10)
         WRITE(wrtstr,'(a,i4)') '        iostat = ', ierr
         CALL xwrite(wrtstr, 10)
         CALL frelun(lun)
         RETURN
      ENDIF

      WRITE(wrtstr, *) (chvals(j),j=1,chnpar)
      WRITE(lun, '(a)') wrtstr(:lenact(wrtstr))

c loop over the chain including the burn-in. We accumulate the chain
c through the burn-in then reset the start point so the burn-in data
c is overwritten and the chain ends up with length chlen.

      DO i = 2, chlen+chburn

c generate the proposed new set of parameters - use a multivariate Gaussian
c based on the eigenfunctions of likelihood space at the last calculated
c best fit. This assumes that the chain is being run after a best fit
c procedure with the same data and model.

         DO ipar = 1, nfpar
            parval(ipar,1) = fgpval(ipar, 'v')
         ENDDO

c generate new parameters till we get a valid set

         qdone = .FALSE.
         DO WHILE ( .NOT.qdone )
            CALL simprm(neigen, MEMI(ieftprm), MEMD(ieigenv), 
     &                  MEMD(ievalue), nfpar, parval(1,1), parval(1,2))
            qdone = .TRUE.
            DO ipar = 1, nfpar
               IF ( parval(ipar,2) .LT. fgpval(ipar,'l') .OR.
     &              parval(ipar,2) .GT. fgpval(ipar,'h') ) THEN
                  qdone = .FALSE.
               ENDIF
            ENDDO
         ENDDO

c set the parameters just generated

         DO ipar = 1, nfpar
            CALL fppval(ipar, parval(ipar,2), 'v', istat)
         ENDDO

c calculate the likelihood for these new parameters

         CALL evlmod(memr(iear), .FALSE., memr(iesavar), memr(ifluxar), 
     &               memr(ibflxar), memr(ieffar))
         CALL fldmod(memr(ifluxar), memr(ibflxar), memi(iigroup), 
     &               memi(iichanb), memi(iichane), memr(irespon), 
     &               memr(dgdpnt('mod')), memr(ieffar))
         CALL calstat(memr(dgdpnt('mod')), memr(dgdpnt('src')), 
     &                memr(dgdpnt('csv')), memr(dgdpnt('cor')), 
     &                memr(dgdpnt('bkg')), memr(dgdpnt('cbv')), 
     &                memr(dgdpnt('reg')), memr(dgdpnt('brg')), 
     &                memr(dgdpnt('are')), memr(dgdpnt('bre')), 
     &                ftstat)

c decide whether to include these proposed new parameters

         oldstat = chvals(chnpar)
         qincnw = .FALSE.
         IF ( ftstat .LE. oldstat )  THEN
            qincnw = .TRUE.
         ELSE
            CALL ranlux(random, 1)
            IF ( random .LT. exp(0.5*(oldstat-ftstat)) ) 
     &         qincnw = .TRUE.
         ENDIF

c if we include the new parameters then save in the chain

         IF ( qincnw ) THEN

            iv = 0
            DO ipar = 1, nfpar
               IF ( .NOT. fgqfrz(ipar) ) THEN
                  iv = iv + 1
                  chvals(iv) = fgpval(ipar, 'v')
               ENDIF
            ENDDO
            chvals(chnpar) = ftstat

c if we don't include the new parameters then
c reset the parameters to the old values

         ELSE

            DO ipar = 1, nfpar
               CALL fppval(ipar, parval(ipar,1), 'v', istat)
            ENDDO

         ENDIF

c if this is the end of the burn-in then reset to the start of the file

         IF ( i .EQ. chburn+1 ) REWIND(lun)

c Write the current chain values to the file

         WRITE(wrtstr, *) (chvals(j),j=1,chnpar)
         WRITE(lun, '(a)') wrtstr(:lenact(wrtstr))

c end loop over the chain

      ENDDO

c Close the output file

      CLOSE(lun)
      CALL frelun(lun)

c free up memory from the temporary parameters

      CALL udmfre( ieftprm, 4, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Memory deallocation of EFTPRM failed in RUNCHN.', 2)
      CALL udmfre( ievalue, 7, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Memory deallocation of EVALUE failed in RUNCHN.', 2)
      CALL udmfre( ieigenv, 7, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     & ' WARNING: Memory deallocation of EIGENV failed in RUNCHN.', 2)

c reset to the original parameters and recalculate model

      DO i = 1, nfpar
         CALL fppval(i, parasf(i), 'v', istat)
      ENDDO

      CALL evlmod(memr(iear), .FALSE., memr(iesavar), memr(ifluxar), 
     &            memr(ibflxar), memr(ieffar))
      CALL fldmod(memr(ifluxar), memr(ibflxar), memi(iigroup), 
     &            memi(iichanb), memi(iichane), memr(irespon), 
     &            memr(dgdpnt('mod')), memr(ieffar))
      CALL calstat(memr(dgdpnt('mod')), memr(dgdpnt('src')), 
     &             memr(dgdpnt('csv')), memr(dgdpnt('cor')), 
     &             memr(dgdpnt('bkg')), memr(dgdpnt('cbv')), 
     &             memr(dgdpnt('reg')), memr(dgdpnt('brg')), 
     &             memr(dgdpnt('are')), memr(dgdpnt('bre')), 
     &             ftstat)

      RETURN
      END