C @(#)lfit.for	19.1 (ES0-DMD) 02/25/03 13:30:46
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C   SUBROUTINE  LFIT                VER. 1.0                MAY  17, 1986
C   P.GROSBOL                       ESO/GARCHING
C
C.PURPOSE
C   FIT A STRAIGHT LINE TO A SET OF DATA.
C   ALGORITHM IS TAKEN FROM BEVINGTON 'LINFIT' ROUTINE
C
C.INPUT/OUTPUT
C   THE PARAMETERS ARE
C       X(N)    :  ARRAY WITH INDEPENDENT VARIABLE
C       Y(N)    :  ARRAY WITH DEPENDENT VARIABLE
C       W(I)    :  ARRAY WITH WEIGHTS
C       N       :  NO. OF DATA SETS
C       MODE    :  MODE IS 0:NO WEIGHTS, 1: USE WEIGHTS
C       A       :  CONSTANT IN   Y = A + B*X
C       SA      :  SIGMA FOR A
C       B       :  SLOPE IN EQUATION :  Y = A + B*X
C       SB      :  SIGMA FOR B
C       R       :  REGRESSION COEFFICIENT
C
C-------------------------------------------------------------------------
      SUBROUTINE LFIT(X,Y,W,N,MODE,A,SA,B,SB,R)
C
      IMPLICIT  NONE
      INTEGER   N
      REAL      X(N)
      REAL      Y(N)
      REAL      W(N)
      INTEGER   MODE
      REAL      A
      REAL      SA
      REAL      B
      REAL      SB
      REAL      R
C
      INTEGER   I
      REAL      CN
      REAL      S,SX,SY,SXX,SXY,SYY
      REAL      XX,YY
      REAL      DEL
      REAL      WW,VAR
C
C     ACCUMULATE SUMS
C
      S   = 0.0
      SX  = 0.0
      SY  = 0.0
      SXX = 0.0
      SXY = 0.0
      SYY = 0.0
      DO 100, I = 1,N
         XX = X(I)
         YY = Y(I)
         IF (MODE.EQ.0) THEN
            WW = 1.0
         ELSE
            WW = W(I)
         ENDIF
         S   = S   + WW
         SX  = SX  + WW*XX
         SY  = SY  + WW*YY
         SXX = SXX + WW*XX*XX
         SXY = SXY + WW*XX*YY
         SYY = SYY + WW*YY*YY
  100 CONTINUE
C
C     COMPUTE COEFFICIENTS
C
      DEL = S*SXX - SX*SX
      A   = (SXX*SY - SX*SXY) / DEL
      B   = (SXY*S  - SX*SY) / DEL
      IF (MODE.EQ.0) THEN
         CN  = N - 2
         VAR = (SYY+A*A*S+B*B*SXX-2.0*(A*SY+B*SXY+A*B*SXX))/CN
      ELSE
         VAR = 1.0
      ENDIF
      SA  = SQRT( VAR*SXX/DEL )
      SB  = SQRT( VAR*S/DEL )
      R   = (S*SXY - SX*SY) / SQRT(DEL*(S*SYY-SY*SY))
      RETURN
      END
