C @(#)gradet.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15:41
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.IDENTIFICATION
C  subroutine GRADET         version 1.1       830908
C  A. Kruszewski             ESO Garching
C.KEYWORDS
C  gradient regression
C.PURPOSE
C  calculates regression coefficient between vectors "X" and "Y"
C  only these values of arrays "X" and "Y" are used for which
C  corresponding element of array "IW" is equal 1
C.INPUT/OUTPUT
C  X           real*4 array      independent variable
C  Y           real*4 array      dependent variable
C  IW          integer*4 array   array indicating which values
C                                of arrays "X" and "Y" are to be used
C  N           integer*4         number of elements in "X" and "Y"
C                                minus one
C  output arguments
C  AVGR        real*4           regression coefficient
C  ZRLV        real*4           zero level
C  DVGR        real*4 array     deviations of elements of vector "X"
C                               from the regression line
C  SIGMA       real*4           sigma of single point
C-----------------------------------------------------------------------
      SUBROUTINE GRADET(N, X, Y, IW, AVGR, ZRLV, DVGR, SIGMA)
      IMPLICIT  NONE
C
      INTEGER   N
      REAL      X(0:N)
      REAL      Y(0:N)
      INTEGER   IW(0:N)
      REAL      AVGR
      REAL      ZRLV
      REAL      DVGR(0:N)
      REAL      SIGMA
C
      INTEGER   NN 
      INTEGER   L
      REAL      A11, A12, A22
      REAL      V1, V2, V12
      REAL      B11, B12, B13
      REAL      B22, B222, B23
      REAL      B33
      REAL      C11, C12, C22
C
C ***
      A11 = 0.0
      A12 = 0.0
      A22 = 0.0
      V1 = 0.0
      V2 = 0.0
      V12 = 0.0
      NN = 0
C
C ******      Prepare normal equations.
C
      DO 10 L = 0 , N
            IF ( IW(L) .EQ. 1 ) THEN
                  A11 = A11 + 1
                  A12 = A12 + X(L)
                  A22 = A22 + X(L) * X(L)
                  V1 = V1 + Y(L)
                  V2 = V2 + X(L) * Y(L)
                  V12 = V12 + Y(L) * Y(L)
                  NN = NN + 1
            ENDIF
   10      CONTINUE
C
C ******      Invert matrix of coefficients in normal equation.
C
      B11 = SQRT( A11 )
      B12 = A12 / B11
      B13 = V1 / B11
      B222 = A22 - B12 * B12
      IF ( B222 .GT. 0.0 ) THEN
            B22 = SQRT( B222 )
            B23 = ( V2 - B12 * B13 ) / B22
            B33 = V12 - B13 * B13 - B23 * B23
            C11 = 1.0 / B11
            C22 = 1.0 / B22
            C12 = -(C11*B12) / B22
C
C ******            Calculate unknowns ZRLV and AVGR.
C
            ZRLV = C11*B13 + C12*B23
            AVGR = C22*B23
C
C ******            Calculates mean error of
C ******            single data point SIGMA. 
C
            IF ( NN .GT. 2 .AND. B33 .GT. 0.0 ) THEN
                  B33 = SQRT( B33 / (NN-2) )
                  SIGMA = B33
            ELSE
                  SIGMA = 0.0
            ENDIF
C
C ******            Calculate deviations of data points
C ******            from the fitted regression line.
C
            IF ( SIGMA .GT. 0.0 ) THEN
                  DO 20 L = 0 , N
                        DVGR(L) = Y(L) - ZRLV - X(L)*AVGR
   20                  CONTINUE
               ELSE
                  DO 30 L = 0 , N
                        DVGR(L) = 0.0
   30                  CONTINUE
               ENDIF
      ELSE
            AVGR = 0.0
            SIGMA = 0.0
      ENDIF
C
      RETURN
      END