SUBROUTINE LFIT3(X1,X2,X3,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,
     *                 NCVM,CHISQ,FUNCS,ISTAT)
*
*  Module number: FOS/HRS utility
*
*  Module name: lfit2
*
*  Keyphrase:
*  ----------
*	Linear Least Squares Fit
*
*  Description:
*  ------------
*		Numerical Recipe Routine modified to allow 3 independent
*	variables.
*	Given data points X1(i),X2(i) and Y3(i) with standard deviations
*	SIG(i), minimization of chi squared is used to determine the coef.
*	A where yfit = sum over j of A(j) x AFUNCj(X1,X2,X3)
*
*  FORTRAN name: lfit3
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  	none
*
*  Subroutines Called:
*  -------------------
*  Others:
*	FOS/HRS utility routines COVSRT, GAUSSJ
*
*  History:
*  --------
*  Version      Date        Author          Description
*	1	Aug 87	D. Lindler	Conversion of Numerical Recipe Routine
*-------------------------------------------------------------------------------
C
C INPUT PARAMETERS
C
C	X1 - VECTOR OF X-VALUES (REAL*8)
C	X2 - SECOND VECTOR OF X-VALUES (REAL*8)
C	X3 - THIRD VECTOR OF X-VALUES (REAL*8)
C	Y - VECTOR OF Y-VALUES (REAL*8)
C	SIG - VECTOR OF STANDARD DEV. FOR Y (REAL*8)
C	NDATA - NUMBER OF POINTS IN X,Y, AND SIG (INTEGER)
C	A - VECTOR OF COEFFICIENTS
C	MA - NUMBER OF COEFFICIENTS IN A (INTEGER)
C	LISTA - LIST OF COEFFICIENTS TO FIT (INTEGER)
C	MFIT - NUMBER OF COEFFICIENTS TO FIT
C	NCVM - DIMENSION OF COVARIANCE MATRIX
C	FUNCS - NAME OF FUNCTION SUBROUTINE
C		CALL FUNCS(X1,X2,X3,AFUNC,MA) RETURNS FUNCTION VALUE FOR
C		EACH TERM IN THE EQUATION IN AFUNC
C
C OUTPUT PARAMETERS
C
C	A - VECTOR OF FITTED COEFFICIENTS
C	COVAR - OUTPUT COVARIANCE MATRIX (REAL*8)
C	CHISQ - CHI SQUARED OF THE FIT
C
C------------------------------------------------------------------------
        DOUBLE PRECISION X1,X2,Y,SIG,A,COVAR,BETA,AFUNC,SUM,CHISQ,
     *                  YM,SIG2I,WT,X3
        INTEGER LISTA,I,J,KK,MFIT,ISTAT,K,NDATA,NCVM,IHIT,MA
        DIMENSION X1(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA),
     *          COVAR(NCVM,NCVM),BETA(50),AFUNC(50),X2(NDATA),X3(1)
        EXTERNAL FUNCS
        KK=MFIT+1
        DO 20 J=1,MA
                IHIT=0
                DO 10 K=1,MFIT
                        IF (LISTA(K).EQ.J) IHIT=IHIT+1
10              CONTINUE
                IF (IHIT.EQ.0) THEN
                        LISTA(KK)=J
                        KK=KK+1
                ENDIF
20      CONTINUE
        DO 40 J=1,MFIT
                DO 30 K=1,MFIT
                        COVAR(J,K)=0.
30              CONTINUE
                BETA(J)=0.
40      CONTINUE
        DO 80 I=1,NDATA
                CALL FUNCS(X1(I),X2(I),X3(I),AFUNC,MA)
                YM=Y(I)
                IF(MFIT.LT.MA) THEN
                        DO 50 J=MFIT+1,MA
                                YM=YM-A(LISTA(J))*AFUNC(LISTA(J))
50                      CONTINUE
                ENDIF
                SIG2I=1./SIG(I)**2
                DO 70 J=1,MFIT
                        WT=AFUNC(LISTA(J))*SIG2I
                        DO 60 K=1,J
                                COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K))
60                      CONTINUE
                        BETA(J)=BETA(J)+YM*WT
70              CONTINUE
80      CONTINUE
        IF (MFIT.GT.1) THEN
                DO 100 J=2,MFIT
                        DO 90 K=1,J-1
                                COVAR(K,J)=COVAR(J,K)
90                      CONTINUE
100             CONTINUE
        ENDIF
        CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1,ISTAT)
        IF(ISTAT.GT.0) GO TO 999
        DO 110 J=1,MFIT
                A(LISTA(J))=BETA(J)
110     CONTINUE
        CHISQ=0.
        DO 130 I=1,NDATA
                CALL FUNCS(X1(I),X2(I),X3(I),AFUNC,MA)
                SUM=0.
                DO 120 J=1,MA
                        SUM=SUM+A(J)*AFUNC(J)
120             CONTINUE
                CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
130     CONTINUE
        CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
999     RETURN
        END