C @(#)tdregl.for	19.1 (ESO-DMD) 02/25/03 14:11:19
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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.1  ESO-FORTRAN Conversion, AA  15:59 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION   program TDREGL
C
C.PURPOSE
C  execute the command
C  REGRESSION/LINEAR  table  y[,weight] x1,x2,x3... [f1,f2]
C
C.KEYWORDS
C  REGRESSION, TABLES
C
C.ALGORITHM
C  uses table interface routines
C
C.VERSION
C 020807	last modif
C 
C-----------------------------------------------------------
C
      SUBROUTINE TDREGL
      IMPLICIT NONE
C
      LOGICAL SEL,NEXT,NULL(10)
C
      INTEGER         IVAR(10),NVMAX
      INTEGER         NPAR,ISTAT,I,TID,NCOL,NROW,NSC1,NAC,NAR
      INTEGER         INDX,IW,LL,NAV,IND,IDEP,NIND,NN,I1,I2,NDEP,NTOT
      INTEGER         IROW,N1,J,IDX,STATUS
      INTEGER         PNEVAL, IPAR(13)
      INTEGER         INDEX, KNUL, KUN, TYPE
C
      DOUBLE PRECISION VAL(10),REGM(81),WM(10),TW,TOL,COE(9)
      DOUBLE PRECISION DPAR(20)
C
      REAL            F1,F2,SERR(9),RPAR(13)
C
      CHARACTER*80 TABLE
      CHARACTER*8  METHOD, FORM
      CHARACTER*17 COLUMN
      CHARACTER*80 LINE
      CHARACTER*64 SELE
      CHARACTER    CPAR*20
      CHARACTER*16 MSG
      CHARACTER*90 SDUM
C
       INCLUDE 'MID_INCLUDE:TABLES.INC'
       INCLUDE 'MID_INCLUDE:TABLED.INC'
      DATA           NVMAX/8/
      DATA           CPAR/' '/
      DATA           PNEVAL/4/
      DATA           MSG/'ERR:TREGRLINxxxx'/
C
C ... GET INPUT PARAMETERS
C
      CALL TDPGET(PNEVAL,NPAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 80
      IF (NPAR.LT.3) THEN
          ISTAT  = ERRPAR
          GO TO 80
      END IF
      TABLE  = TPARBF(1)
      METHOD = 'LINEAR'
      DO 10 I = 1,81
          REGM(I) = 0.D0
   10 CONTINUE
      DO 20 I = 1,10
          WM(I)  = 0.D0
          NULL(I) = .FALSE.
   20 CONTINUE
      TW     = 0.D0
C
C ... INIT INPUT TABLE
C
      TID = -1
      CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 80
      CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 80
      CALL TDRSEL(TID,SELE,ISTAT)
C
C ... FIND COLUMNS WITH VARIABLES
C
      INDX   = INDEX(TPARBF(2),',')
      IF (INDX.EQ.0) THEN
          COLUMN = TPARBF(2)
          IW     = 0
      ELSE
          COLUMN = TPARBF(2) (INDX+1:)
          CALL TBCSER(TID,COLUMN,IW,ISTAT)
          IF (IW.EQ.-1) THEN
              ISTAT  = ERRPAR
              GO TO 80
          END IF
          CALL TBFGET(TID,IW,FORM,LL,TYPE,ISTAT)
          IF (TYPE.EQ.D_C_FORMAT) THEN
              ISTAT  = ERRFMT
              GO TO 80
          END IF
          COLUMN = TPARBF(2) (1:INDX-1)
      END IF
      CALL TBCSER(TID,COLUMN,IDEP,ISTAT)
      IF (IDEP.EQ.-1) THEN
          ISTAT  = ERRPAR
          GO TO 80
      END IF
      CALL TBFGET(TID,IDEP,FORM,LL,TYPE,ISTAT)
      IF (TYPE.EQ.D_C_FORMAT) THEN
          ISTAT  = ERRFMT
          GO TO 80
      END IF
      NIND   = 0
      LINE   = ' '
      LINE   = TPARBF(3)
      NN     = INDEX(LINE,' ') - 1
      I1     = 1
      SDUM   = LINE(1:NN)//','
   30 CONTINUE
      I2     = INDEX(SDUM,',') - 1
      SDUM(I2+1:I2+1) = ';'
      COLUMN = LINE(I1:I2)
      CALL TBCSER(TID,COLUMN,IND,ISTAT)
      IF (ISTAT.NE.0) GO TO 80
      IF (IND.EQ.-1) THEN
          ISTAT  = ERRPAR
          GO TO 80
      END IF
      CALL TBFGET(TID,IDEP,FORM,LL,TYPE,ISTAT)
      IF (TYPE.EQ.D_C_FORMAT) THEN
          ISTAT  = ERRFMT
          GO TO 80
      END IF
      NIND   = NIND + 1
      IVAR(NIND) = IND
      I1     = I2 + 2
      IF (I1.LT.NN) GO TO 30
      IF (NIND.GT.NVMAX) THEN
          ISTAT  = ERRPAR
          GO TO 80
      END IF
      NDEP   = NIND + 1
      IVAR(NDEP) = IDEP
      NTOT   = NDEP
C
C ... INCLUDE WEIGHTS
C
      IF (IW.EQ.0) THEN
          IW     = NDEP + 1
          VAL(IW) = 1.D0
      ELSE
          NTOT   = NTOT + 1
          IVAR(NTOT) = IW
          IW     = NTOT
      END IF
C
C ... ITERATE ON ROWS
C
      N1     = 0
      DO 60 IROW = 1,NROW
          CALL TBSGET(TID,IROW,SEL,ISTAT)
          IF (SEL) THEN
              CALL TBRRDD(TID,IROW,NTOT,IVAR,VAL,NULL,ISTAT)
              NEXT   = ( .NOT. NULL(1)) .AND. ( .NOT. NULL(2)) .AND.
     +                 ( .NOT. NULL(3)) .AND. ( .NOT. NULL(4)) .AND.
     +                 ( .NOT. NULL(5)) .AND. ( .NOT. NULL(6)) .AND.
     +                 ( .NOT. NULL(7)) .AND. ( .NOT. NULL(8)) .AND.
     +                 ( .NOT. NULL(9))
              IF (NEXT) THEN
                  N1     = N1 + 1
                  DO 50 I = 1,NDEP
                      WM(I)  = WM(I) + VAL(IW)*VAL(I)
                      DO 40 J = 1,NDEP
                          IDX    = (J-1)*NDEP + I
                          REGM(IDX) = REGM(IDX) + VAL(IW)*VAL(I)*VAL(J)
   40                 CONTINUE
   50             CONTINUE
                  TW     = TW + VAL(IW)
              END IF
          END IF
   60 CONTINUE
      IF (TW.LT.DBLE(NDEP+1)) THEN
          CALL STTPUT(' Error in no. of data...',ISTAT)
          GO TO 80
      END IF
C
C ... DO REGRESSION ANALYSIS
C
      CALL STKRDR('INPUTR',1,1,NAV,F1,KUN,KNUL,ISTAT)
      F2     = F1
      TOL    = 1.D-8
      CALL TDMLRG(REGM,NDEP,WM,TW,COE,SERR,F1,F2,TOL)
C
C ... COPY RESULT OF REGRESSION
C
      DO 65 I = 1, 20
        DPAR(I) = 0.D0
65    CONTINUE
      CPAR(9:16) = TABLE
      CPAR(17:20) = 'LINE'
      IPAR(1) = N1
      IPAR(2) = NIND
      IPAR(3) = IDEP
      RPAR(1) = 0.
      RPAR(2) = 0.
      RPAR(3) = SERR(NDEP)
      DPAR(1) = COE(NDEP)
      DO 70 I = 1,NIND
          IPAR(3+I) = IVAR(I)
          RPAR(3+I) = SERR(I)
          DPAR(1+I) = COE(I)
   70 CONTINUE
      CALL STKWRC('OUTPUTC',1,CPAR,1,20,KUN,ISTAT)
      CALL STKWRI('OUTPUTI',IPAR,1,13,KUN,ISTAT)
      CALL STKWRR('OUTPUTR',RPAR,1,13,KUN,ISTAT)
      CALL STKWRD('OUTPUTD',DPAR,1,20,KUN,ISTAT)
C
C ... PRINT RESULT
C
      IF (SELE(1:1).NE.'-') THEN
          LINE   = ' SELECT  '//SELE
          CALL STTPUT(LINE,ISTAT)
      END IF
      CALL TDRDS1(CPAR,IPAR,RPAR,DPAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 80
C
C ... end
C
      CALL DSCUPT(TID,TID,' ',ISTAT)
      CALL TBTCLO(TID,ISTAT)
   80 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,STATUS)
      END IF
      RETURN
 9000 FORMAT (I4)
      END
