C @(#)lfit.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:00 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