C===========================================================================
C Copyright (C) 1995-2005 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.0  ESO-FORTRAN Conversion, AA  16:02 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  program TDREGP
C
C.PURPOSE
C
C  Execute the command
C  REGRESSION/POL table dep_var ind_var1[,ind_var2] deg1[,deg2]  [outkey]
C  OR
C  REGRESSION/POL table dep_var ind_var1[,ind_var2,...]   ?      [outkey]
C
C.KEYWORDS
C
C  Tables, Multiple linear regression, POLYNOMIAL FIT
C
C.ALGORITHM
C
C  uses table interface routines
C  fitting algorithm using Householder transformation
C
C.VERSION
C 051019	last modif
C 
C-----------------------------------------------------------
C
      SUBROUTINE TDREGP
      IMPLICIT NONE
C
      LOGICAL NULL(9), SEL, IPOL, NEXT
C
      INTEGER PARVAL
      INTEGER IVAR(9),IDEG(2),IPAR(13)
      INTEGER NROW,NCOL, NAC, NAR, I, ISTAT, STATUS
      INTEGER IAV, TID, NSC1, NN, NVAR, I1, I2, K, L, N, LL, N1
      INTEGER IROW, KMIN, NPAR, NM1, K1
      INTEGER INDEX, KUN, KNUL, TYPE, MAXPAR, INDKEY
      INTEGER NOELM,BYTELM,DEGMAX
      PARAMETER (MAXPAR=100)
C
      REAL ERROR,RPAR1(4),RPAR(13)
C
      DOUBLE PRECISION XMIN,XMAX,YMIN,YMAX
      DOUBLE PRECISION VAL(9)
      DOUBLE PRECISION G(MAXPAR,MAXPAR),X(MAXPAR)
C
      CHARACTER*16 MSG
      CHARACTER*80 TABLE,LINE,SELE
      CHARACTER*90 SDUM
      CHARACTER*17 COLUMN(9)
      CHARACTER*8  FORM,TYPKEY
      CHARACTER*20 CPAR, KEYNAM, INPNAM
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C 
      DATA MSG/'ERR:TREGRMUL'/
      DATA PARVAL/5/
      DATA CPAR/' '/
      DATA IPAR/13*0/
      DATA RPAR/13*0./
C
C ... GET INPUT PARAMETERS
C
      DO 10, I=1,9
          NULL(I) = .FALSE.
          IVAR(I) = 0
          VAL(I) = 0
   10 CONTINUE

      CALL TDPGET(PARVAL,NPAR,ISTAT)

      IF (ISTAT.NE.0) GO TO 70
      IF (NPAR.LT.4) THEN
          ISTAT  = ERRPAR
          GO TO 70
      END IF

      TABLE  = TPARBF(1)
      KEYNAM = TPARBF(5)
      INDKEY = INDEX(KEYNAM,' ') - 1

      DO 15 I = 1, MAXPAR
          X(I) = 0.D0
   15 CONTINUE
      INPNAM = KEYNAM(1:INDKEY)//'D'
      CALL STKFND(INPNAM,TYPKEY,NOELM,BYTELM,ISTAT)
      CALL STKWRD(INPNAM,X,1,NOELM,KUN,ISTAT)

      CALL STKRDR('INPUTR',1,2,IAV,RPAR1,KUN,KNUL,ISTAT)
      IDEG(1) = RPAR1(1)
      IDEG(2) = RPAR1(2)

      DEGMAX = 0
      IF (IDEG(2).GT.0)   DEGMAX = IDEG(1)

      IF (DEGMAX.GT.8) THEN
         CALL STETER(9,'REGRES/POLY: degree along X is limited to 8')
      ENDIF

      IF (IDEG(2).GT.8) THEN
         CALL STETER(9,'REGRES/POLY: degree along Y is limited to 8')
      ENDIF
C
C ... READ INPUT TABLE
C
      TID = -1
      CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
      CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
      CALL TDRSEL(TID,SELE,ISTAT)
C
C ... DECODE PARAMETERS
C
      COLUMN(1) = TPARBF(2)
      LINE   = ' '
      LINE   = TPARBF(3)
      NN     = INDEX(LINE,' ') - 1
      NVAR   = 1
      I1     = 1
      SDUM   = LINE//','
   20 CONTINUE
      I2     = INDEX(SDUM,',')
      SDUM(I2:I2) = ';'
      NVAR   = NVAR + 1
      IF (NVAR.GT.9) THEN
          ISTAT  = ERRPAR
          GO TO 70
      END IF
      COLUMN(NVAR) = LINE(I1:I2-1)
      I1     = I2 + 1
      IF (I1.LT.NN) GO TO 20
      K      = IDEG(1)
      L      = IDEG(2)
      IF (K.EQ.0 .AND. L.EQ.0) THEN
          IPOL   = .FALSE.
          N      = NVAR
          NM1    = N - 1
      ELSE
          IPOL   = .TRUE.
          N      = (K+1)* (L+1)
      END IF
      IF (N.GT.NOELM) THEN
          ISTAT  = ERRPAR
          GO TO 70
      END IF
C
C ... SEARCH FOR COLUMNS AND CHECK FORMAT
C
      DO 30 I = 1,NVAR
          CALL TBCSER(TID,COLUMN(I),IVAR(I),ISTAT)
          IF (ISTAT.NE.0) GO TO 70
          IF (IVAR(I).EQ.-1) THEN
              ISTAT  = ERRPAR
              GO TO 70
          END IF
          CALL TBFGET(TID,IVAR(I),FORM,LL,TYPE,ISTAT)
          IF (TYPE.EQ.D_C_FORMAT) THEN
              ISTAT  = ERRFMT
              GO TO 70
          END IF
   30 CONTINUE
      CALL STTPUT(TABLE,ISTAT)
C
C ... SOLVE LEAST SQUARE PROBLEM
C
      XMIN   = 1.D30
      YMIN   = 1.D30
      XMAX   = -XMIN
      YMAX   = -YMIN
      LL     = 0
      N1     = 0
      DO 50 IROW = 1,NROW
          CALL TBSGET(TID,IROW,SEL,ISTAT)
          IF (SEL) THEN
              CALL TBRRDD(TID,IROW,NVAR,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
                  IF (IPOL) THEN
                      XMIN   = MIN(XMIN,VAL(2))
                      XMAX   = MAX(XMAX,VAL(2))
                      YMIN   = MIN(YMIN,VAL(3))
                      YMAX   = MAX(YMAX,VAL(3))
                      CALL TDSET2(LL+1,VAL(2),VAL(3),VAL(1),K,L,G,X,N,
     +                            MAXPAR)
                  ELSE
                      CALL TDSET1(LL+1,VAL(2),VAL(1),NM1,G,X,N,MAXPAR)
                  END IF
                  IF (LL.NE.0) THEN
                      KMIN   = MIN(LL,N+1)
                      DO 40 K1 = 1,KMIN
                          CALL TDHHTR(K1,LL+1,G,X,N,MAXPAR)
   40                 CONTINUE
                  END IF
                  LL     = MIN(LL+1,N+1)
              END IF
          END IF
   50 CONTINUE
      CALL TDSOLV(G,X,N,MAXPAR)
      ERROR  = ABS(G(N+1,N+1))
C
C ... COPY RESULT OF REGRESSION
C
      CPAR(9:16) = TABLE
      IPAR(1) = N1
      IPAR(2) = NVAR - 1
      IPAR(3) = IVAR(1)
      IF (IPOL) THEN
          CPAR(17:20) = 'MULT'
          IPAR(4) = IVAR(2)
          IPAR(5) = IVAR(3)
          IPAR(6) = IDEG(1)
          IPAR(7) = IDEG(2)
      ELSE
          CPAR(17:20) = 'LINE'
          DO 60 I = 1,NVAR - 1
              IPAR(I+3) = IVAR(I+1)
   60     CONTINUE
      END IF
      RPAR(1) = XMIN
      RPAR(2) = XMAX
      RPAR(3) = YMIN
      RPAR(4) = YMAX
      RPAR(5) = SQRT((ERROR*ERROR)/IPAR(1))

      INPNAM = KEYNAM(1:INDKEY)//'C'
      CALL STKWRC(INPNAM,1,CPAR,1,20,KUN,ISTAT)
      INPNAM = KEYNAM(1:INDKEY)//'I'
      CALL STKWRI(INPNAM,IPAR,1,13,KUN,ISTAT)
      INPNAM = KEYNAM(1:INDKEY)//'R'
      CALL STKWRR(INPNAM,RPAR,1,5,KUN,ISTAT)
      INPNAM = KEYNAM(1:INDKEY)//'D'
      CALL STKWRD(INPNAM,X,1,N,KUN,ISTAT)
C
C ... PRINT RESULT
C
      IF (SELE(1:1).NE.'-') THEN
          LINE   = ' SELECT  '//SELE
          CALL STTPUT(LINE,ISTAT)
      END IF
      RPAR(5) = ERROR
      CALL TDRDIS(CPAR,IPAR,RPAR,X,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
C
C ... END
C
      CALL DSCUPT(TID,TID,' ',ISTAT)
      CALL TBTCLO(TID,ISTAT)
   70 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,STATUS)
      END IF
      RETURN
 9000 FORMAT (I4)
      END
