SUBROUTINE lfit(x,y,sig,a,maska,covar,chisq,funcs)
	USE nrtype; USE nrutil, ONLY : assert_eq,diagmult,nrerror
	USE nr, ONLY :covsrt,gaussj
	IMPLICIT NONE
	REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig
	REAL(SP), DIMENSION(:), INTENT(INOUT) :: a
	LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska
	REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: covar
	REAL(SP), INTENT(OUT) :: chisq
	INTERFACE
		SUBROUTINE funcs(x,arr)
		USE nrtype
		IMPLICIT NONE
		REAL(SP),INTENT(IN) :: x
		REAL(SP), DIMENSION(:), INTENT(OUT) :: arr
		END SUBROUTINE funcs
	END INTERFACE
	INTEGER(I4B) :: i,j,k,l,ma,mfit,n
	REAL(SP) :: sig2i,wt,ym
	REAL(SP), DIMENSION(size(maska)) :: afunc
	REAL(SP), DIMENSION(size(maska),1) :: beta
	n=assert_eq(size(x),size(y),size(sig),'lfit: n')
	ma=assert_eq(size(maska),size(a),size(covar,1),size(covar,2),'lfit: ma')
	mfit=count(maska)
	if (mfit == 0) call nrerror('lfit: no parameters to be fitted')
	covar(1:mfit,1:mfit)=0.0
	beta(1:mfit,1)=0.0
	do i=1,n
		call funcs(x(i),afunc)
		ym=y(i)
		if (mfit < ma) ym=ym-sum(a(1:ma)*afunc(1:ma), mask=.not. maska)
		sig2i=1.0_sp/sig(i)**2
		j=0
		do l=1,ma
			if (maska(l)) then
				j=j+1
				wt=afunc(l)*sig2i
				k=count(maska(1:l))
				covar(j,1:k)=covar(j,1:k)+wt*pack(afunc(1:l),maska(1:l))
				beta(j,1)=beta(j,1)+ym*wt
			end if
		end do
	end do
	call diagmult(covar(1:mfit,1:mfit),0.5_sp)
	covar(1:mfit,1:mfit)= &
		covar(1:mfit,1:mfit)+transpose(covar(1:mfit,1:mfit))
	call gaussj(covar(1:mfit,1:mfit),beta(1:mfit,1:1))
	a(1:ma)=unpack(beta(1:ma,1),maska,a(1:ma))
	chisq=0.0
	do i=1,n
		call funcs(x(i),afunc)
		chisq=chisq+((y(i)-dot_product(a(1:ma),afunc(1:ma)))/sig(i))**2
	end do
	call covsrt(covar,maska)
	END SUBROUTINE lfit