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
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  15:55 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C  program TDCMPR
C
C.PURPOSE
C  creates a new column with values computed using coefficients
C  generated in a regression command
C  COMP/REGR table output = coefs[(col1,...)]  [R*8]
C
C.KEYWORDS
C  REGRESSION, TABLES
C
C.ALGORITHM
C use table interface routines
C
C.VERSION
C 051017	last modif
C 
C-----------------------------------------------------------
C
      SUBROUTINE TDCMPR
      IMPLICIT NONE
C
      LOGICAL NULL(9),NEXT,IPOL,NEWCOL
C
      INTEGER PARVAL,IPAR(13),ICOL(10)
      INTEGER ISTAT, NPAR, I, TID, NCOL, NROW, NSC, NAC, NAR
      INTEGER OCOL, LL, INDX, IAV, NIND, NIND1, NA, K, L
      INTEGER NN, K1, IP, I1, I2, DUNIT, DNULL
      INTEGER INDVAR, L1, J
      INTEGER INDEX, TYPE, TYPE1
      INTEGER EC, EL, ED, N1, N2
C
      DOUBLE PRECISION DPAR(100),VAL(9),WORK(100),RESULT,DC
C
      CHARACTER*80 TABLE
      CHARACTER*80 LINE1, LINE
      CHARACTER*20 DESNAM, NAME
      CHARACTER*8  FORM, OUTFMT
      CHARACTER*18 OUTCOL
      CHARACTER*80 MSG
      CHARACTER*16 UNIT, LABEL, CTYP
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C 
      DATA PARVAL/5/
      DATA MSG/'ERR:TREGRCOMxxxx'/
      DATA UNIT/' '/,LABEL/' '/
C
C ... GET PARAMETERS
C
      K = 0
      L = 0
      CALL TDPGET(PARVAL,NPAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 80
      TABLE  = TPARBF(1)
      OUTCOL = TPARBF(2)
      LINE   = ' '
      LINE   = TPARBF(4)
      OUTFMT = TPARBF(5)
      DO 10 I = 1,9
         NULL(I) = .FALSE.
         VAL(I) = 0.D0
   10 CONTINUE
      TID = -1
C
C ... read table and check output column
C
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      CALL TBCSER(TID,OUTCOL,OCOL,ISTAT)
      IF (OCOL.EQ.-1) THEN
          LABEL  = OUTCOL(2:)
          IF (OUTFMT(3:3) .EQ. '8') THEN
             TYPE   = D_R8_FORMAT
          ELSE
             TYPE   = D_R4_FORMAT
          ENDIF
          FORM   = 'G14.7'
          CALL TBCINI(TID,TYPE,1,FORM,UNIT,LABEL,OCOL,ISTAT)
      ELSE
          CALL TBFGET(TID,OCOL,FORM,LL,TYPE,ISTAT)
          IF (TYPE.EQ.D_C_FORMAT) THEN
              ISTAT  = ERRFMT
              GO TO 80
          END IF
      END IF
C
C ... read regression parameters
C
      INDX   = INDEX(LINE,'(')
      IF (INDX.EQ.0) THEN
          NEWCOL = .FALSE.
          INDX   = INDEX(LINE,' ') - 1
      ELSE
          NEWCOL = .TRUE.
          INDX   = INDX - 1
      END IF
      DESNAM = LINE(1:INDX)//'C '
      CALL STECNT('GET',EC,EL,ED)
      CALL STECNT('PUT',1,0,0)            !disable errors ...
      CALL STDRDC(TID,DESNAM,1,17,4,IAV,NAME,DUNIT,DNULL,ISTAT)
      CALL STECNT('PUT',EC,EL,ED)
      IF (ISTAT.NE.0) THEN
         CALL STDFND(TID,DESNAM,CTYP,N1,N2,ISTAT)
         IF ((CTYP(1:1) .NE. 'C') .OR. (N1 .LT. 4)) THEN
            CALL STDRDC(TID,DESNAM,1,1,N2,IAV,NAME,DUNIT,DNULL,ISTAT)
         ELSE
            CALL STDRDC(TID,DESNAM,1,N1-3,N1,IAV,NAME,DUNIT,DNULL,ISTAT)
            WRITE(MSG,123) DESNAM,N1
            CALL STTPUT('!OJO! compute/regression: ',ISTAT)
123         FORMAT('!OJO! descr. ',A,' has size',i3,' (not 20)')
            CALL STTPUT(MSG,ISTAT)
         ENDIF
      ENDIF
      DESNAM = LINE(1:INDX)//'I '
      CALL STDRDI(TID,DESNAM,1,13,IAV,IPAR,DUNIT,DNULL,ISTAT)
      NIND   = IPAR(2)
      DO 20, I=1,NIND
          ICOL(I) = IPAR(3+I)
   20 CONTINUE
      IF (NAME(1:4).EQ.'LINE') THEN
          NA     = NIND + 1
          IPOL   = .FALSE.
      ELSE
          K      = IPAR(6)
          L      = IPAR(7)
          NA     = (K+1)* (L+1)
          IPOL   = .TRUE.
      END IF
      DESNAM = LINE(1:INDX)//'D '
      CALL STDRDD(TID,DESNAM,1,NA,IAV,DPAR,DUNIT,DNULL,ISTAT)
C
C ... modify independent variables if required in the syntax
C
      IF (NEWCOL) THEN
          NIND1  = 0
          I1     = INDX + 2
          NN     = INDEX(LINE,')')
          IF (NN.EQ.0) THEN
              ISTAT  = ERRPAR
              GO TO 80
          END IF
          LINE(NN:NN) = ','
          NN     = NN + 1
   30     I2     = INDEX(LINE,',') - 1
          LINE(I2+1:I2+1) = ';'
          OUTCOL = LINE(I1:I2)
          CALL TBCSER(TID,OUTCOL,INDVAR,ISTAT)
          IF (INDVAR.EQ.-1) THEN
              ISTAT  = ERRPAR
              GO TO 80
          END IF
          CALL TBFGET(TID,INDVAR,FORM,LL,TYPE1,ISTAT)
          IF (TYPE1.EQ.D_C_FORMAT) THEN
              ISTAT  = ERRFMT
              GO TO 80
          END IF
          NIND1  = NIND1 + 1
          ICOL(NIND1) = INDVAR
          I1     = I2 + 2
          IF (I1.LT.NN) GO TO 30
          IF (NIND1.NE.NIND) THEN
              ISTAT  = ERRPAR
              GO TO 80
          END IF
      END IF
C
C ... iterate on rows
C
      DO 70 I = 1,NROW
          CALL TBRRDD(TID,I,NIND,ICOL,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
              IF (IPOL) THEN
                  IP     = 0
                  DC     = 1.D0
                  RESULT = 0.D0
                  DO 50 L1 = 0,L
                      IP     = IP + 1
                      WORK(IP) = DC
                      RESULT = RESULT + WORK(IP)*DPAR(IP)
                      DO 40 K1 = 1,K
                          IP     = IP + 1
                          WORK(IP) = WORK(IP-1)*VAL(1)
                          RESULT = RESULT + WORK(IP)*DPAR(IP)
   40                 CONTINUE
                      DC     = DC*VAL(2)
   50             CONTINUE
              ELSE
CC                RESULT = VAL(1)
CC Bug corrected by M.-F. Landrea, Meudon - Paris    23 Oct 1987
                  RESULT = DPAR(1)
                  DO 60 J = 2,NA
                      RESULT = RESULT + VAL(J-1)*DPAR(J)
   60             CONTINUE
              END IF
              CALL TBEWRD(TID,I,OCOL,RESULT,ISTAT)
          ELSE
              CALL TBEDEL(TID,I,OCOL,ISTAT)
          END IF
   70 CONTINUE
C
C ... END
C
      LINE1 = '-'
      CALL TDWSEL(TID,LINE1,ISTAT)
      CALL TDHSTR(TID,ISTAT)
      CALL TBSINI(TID,ISTAT)
      CALL DSCUPT(TID,TID,' ',ISTAT)
      CALL TBTCLO(TID,ISTAT)
   80 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL STTPUT(MSG,ISTAT)
      END IF
      RETURN
C 
 9000 FORMAT (I4)
      END
