      subroutine fchisq(y,sigmay,npts,nfree,yfit,chisq,modef)
      real*8 y(*),sigmay(*),yfit(*),chisq,z
      chisq=0.D0
c
c     Accumulate chi square
c
      do i=1,npts
       if(sigmay(i).ne.0)then
        z=(y(i)-yfit(i))/sigmay(i)
       else
        write(6,*)'***ERROR: data point ',i,' uncertainty = 0!'
        stop
       endif
       if(modef.eq.1)then
        chisq=chisq+abs(z)                   ! Double exponential
       else if(modef.eq.2)then
        chisq=chisq+log(1.d0+z**2/2.d0)*5.d0 ! Lorentzian
       else
        chisq=chisq+z**2                     ! Gaussian noise
       endif
      enddo
      chisq=chisq/dble(nfree)
      return
      end

