C @(#)mtablv.for	19.1 (ESO-DMD) 02/25/03 13:24:28
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.IDENT:    program MTABLV                             A. Lauberts  ESO  890713
C.KEYWORDS: ESOLV data base evaluation, table file, diameter
C.PURPOSE:  get semidiameter for specified magnitude interpolated
C           from photometric parameters stored in ESOLV table file
C.USAGE:    MIDAS application program
C           execute as @@ MTABLV P1 P2 P3 P4 P5
C.PARAMETERS:
C           P1  PCAT   C*12  name of input table file; default PCAT
C           P2  MCOL   I     col no of 1st mag value in array of max 13 elems
C           P3  TCOL   I     col no of reference mag (normally total mag)
C           P4  FRAC   R     fractional (total) light wanted
C           P5  FCOL   I     col no of output semidiameter corresp to FRAC
C.OUTPUT:   semidiameter in column no FCOL
C.COMMENTS: Undefined semidiameters are set to NULL.
C           Algorithm tested against ESOLV values at 50% total light;
C           systematic difference about 1% in radius is due to different
C           interpolation techniques.
C.NOTE:     Uses MATHLIB rotines E01BAF and E02BCF.
C.VERSION:  910207 RHW IMPLICIT NONE added
C-----------------------------------------------------------------------------
       PROGRAM MTABLV
C
       IMPLICIT     NONE
       INTEGER      ACOL, AROW
       INTEGER      COL(256)
       INTEGER      FCOL
       INTEGER      I, IACT, ISTAT
       INTEGER      IFAIL
       INTEGER      J
       INTEGER      KNUL,KUN
       INTEGER      MADRID(1)
       INTEGER      M,MCOL
       INTEGER      NCOL, NROW, NSC, NSEQ
       INTEGER      TCOL
       INTEGER      TID
       INTEGER      TINULL
C
       DOUBLE PRECISION TDNULL
       REAL         FRAC
       REAL         MAG
       REAL         LMD     
       REAL         ROW(256)
       REAL         SMD
       REAL         TRNULL
C
       LOGICAL      NULL(256)
C
       CHARACTER*60 PCAT
       CHARACTER*80 STRING
       CHARACTER*60 TABLE
C
       INCLUDE      'MID_INCLUDE:ST_DEF.INC/NOLIST'
       COMMON       /VMR/MADRID
       INCLUDE      'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
C *** start 
       CALL STSPRO('MTABLV')
C
       DO I=1,256
         COL(I) = I
       ENDDO
C
C  get keywords
       CALL STKRDC('IN_A',1,1,60,IACT,PCAT,KNUL,KUN,ISTAT)
       CALL STKRDI('INPUTI',1,1,IACT,MCOL,KNUL,KUN,ISTAT)
       CALL STKRDI('INPUTI',2,1,IACT,TCOL,KNUL,KUN,ISTAT)
       CALL STKRDI('INPUTI',3,1,IACT,FCOL,KNUL,KUN,ISTAT)
       CALL STKRDR('INPUTR',1,1,IACT,FRAC,KNUL,KUN,ISTAT)
       CALL TBMNUL(TINULL,TRNULL,TDNULL)
C
       FRAC = 2.5*ALOG10(FRAC)
       IF (FCOL.LT.183.OR.FCOL.GT.189) THEN
          STRING = '*** FATAL: Output column not in scratch '//
     2             'interval 183-189'
          CALL STTPUT(STRING,ISTAT)
          CALL STSEPI
       ENDIF
       M = 0
C
C open PCAT table file; don't count last (ascii) column
       CALL TBTOPN(TABLE,F_IO_MODE,TID,ISTAT)
       CALL TBIGET(TID,NCOL,NROW,NSC,ACOL,AROW,ISTAT)
       NCOL = NCOL-1
C
C read PCAT
       DO NSEQ=1,NROW
          M=M+1
          IF (M.EQ.1000) THEN
             M=0
C             WRITE(6,*) NSEQ
          ENDIF
          CALL TBRRDR(TID,NSEQ,NCOL,COL,ROW,NULL,ISTAT)
          IF (.NOT.NULL(TCOL)) THEN
             MAG = ROW(TCOL)-FRAC
             DO I=1,13
                J = MCOL+I-1
                IF (NULL(J)) ROW(J)=-99.
             ENDDO
             CALL IPOL(ROW(MCOL),MAG,LMD,IFAIL)
             IF (IFAIL.EQ.0) THEN
                SMD = 0.5*10.**LMD
             ELSE
                SMD = TRNULL
             ENDIF
          ELSE
             SMD = TRNULL
          ENDIF
          CALL TBEWRR(TID,NSEQ,FCOL,SMD,ISTAT)
       ENDDO
C
C *** finish
       CALL TBTCLO(TID,ISTAT)
       CALL STSEPI
       END



       SUBROUTINE IPOL(AC,MAG,LMD,IFAIL)
C+
C.IDENT:    subroutine IPOL                            A. Lauberts  ESO  881216
C.KEYWORDS  galaxy photometry, database
C.PURPOSE:  extracts interpolated photom data from table file
C.USAGE:    MIDAS application subroutine called from MTABLV
C.PARAMETERS (arguments)
C           input:
C           AC    R   array    magnitudes in 1 to 13 element
C           MAG   R            magnitude
C           output:
C           LMD   R            log major diam corresp to MAG
C           other:
C           DL    R   array    log major diameters
C           BM    R   array    magnitudes
C-
      IMPLICIT         NONE
      REAL             AC(1)
      REAL             MAG
      REAL             LMD
      INTEGER          IFAIL
C
      INTEGER          K
      INTEGER          LCK
      INTEGER          LWRK
      INTEGER          N
      INTEGER          NB
C
      REAL             A
      REAL             BM(13)
      REAL             DL(13)
      REAL             X, X0, X1
      DOUBLE PRECISION CIB(17)
      DOUBLE PRECISION KIB(17)
      DOUBLE PRECISION KLA(13)
      DOUBLE PRECISION KBM(13)
      DOUBLE PRECISION SDIF(4)
      DOUBLE PRECISION WRK(94)
      DOUBLE PRECISION XX
C
      X=1.
      N=0
      DO K=1,13
         IF (AC(K).NE.-99.) THEN
            N=N+1
            DL(N)=X
            BM(N)=AC(K)
         ENDIF
         X=X+0.2
      ENDDO
C
C fill KLA and KBM with data; skip if MAG outside interval or N < 2
C
      IF (MAG.LT.BM(N).OR.MAG.GT.BM(1)) THEN
         IFAIL=1
         CALL FAIL(IFAIL)
         RETURN
      ENDIF
      IF (N.LT.2) THEN
         IFAIL=2
         CALL FAIL(IFAIL)
         RETURN
      ENDIF
      NB=N
      CALL PREPARE(NB,DL,BM,KLA,KBM)
C+
C interpolate with cubic splines on B mags
C (if N<4 then linear or parabolic interpol only,
C  if N=1 then skip interpolation)
C the independent variable is log a, a = major diam
C arr  dim      contains
C input:
C KLA  13       log a
C KBM  13       B integr mag
C output:
C KIB  17       =KBM + 4 extra points
C CIB  17       cubic spline coeffs for KBM
C SDIF  4       value + first three derivatives
C i/o:
C WRK  94       work area for cubic splines
C-
C interpol coeffs
C
      LCK=NB+4
      LWRK=6*NB+16
      IFAIL=0
      IF (NB.GT.3) THEN
         CALL E01BAF(NB,KLA,KBM,KIB,CIB,LCK,WRK,LWRK,IFAIL)
      ELSE
         CALL NEWTON1(NB,KLA,KBM,CIB,IFAIL)
      ENDIF
      IF (IFAIL.NE.0) THEN
        CALL FAIL(IFAIL)
        RETURN
      ENDIF
C
C  value (derivatives not used) found by bisection method
C  return if mag error < 0.01 or no of iterations > 20 (failed)
C
      X0=KLA(1)
      X1=KLA(NB)
      IFAIL=0
      N=0
10    N=N+1
      X=(X0+X1)/2.
      XX=X
      IF (NB.GT.3) THEN
         CALL E02BCF(NB+4,KIB,CIB,XX,0,SDIF,IFAIL)
      ELSE
         CALL NEWTON2(NB,KLA,CIB,XX,SDIF)
      ENDIF
      IF (IFAIL.NE.0.) THEN
         CALL FAIL(IFAIL)
         RETURN
      ENDIF
      A=SDIF(1)
      IF (ABS(A-MAG).LT.0.001) THEN
         LMD=X
         RETURN
      ELSE IF(N.GT.20) THEN
         IFAIL=20
         CALL FAIL(IFAIL)
         RETURN
      ENDIF
      IF (A.LT.MAG) THEN
         X1=X
      ELSE
         X0=X
      ENDIF
      GOTO 10
      END



      SUBROUTINE NEWTON1(N,KLA,KBM,CIB,IFAIL)
C+
C  Newtons interpolation formulas for 2 and 3 points (1st and 2nd degree)
C     P1(X) = FX0 + (X-X0)*(FX1-FX0)/(X1-X0)
C           = CIB(1) + (X-X0)*CIB(2)
C       P1' = CIB(2)
C
C     P2(X) = FX0 + (X-X0)*(FX1-FX0)/(X1-X0) +
C             (X-X0)*(X-X1)*((FX2-FX1)/(X2-X1)-(FX1-FX0)/(X1-X0))/(X2-X0)
C           = CIB(1) + (X-X0)*CIB(2) + (X-X0)*(X-X1)*CIB(3)
C       P2' = CIB(2) + (2X-X0-X1)*CIB(3)
C       P2''= 2*CIB(3) 
C
C determine coefficients CIB for later use in NEWTON2
C case N=2: linear interpolation
C-
      IMPLICIT NONE
      INTEGER  N
      DOUBLE PRECISION  KLA(13)
      DOUBLE PRECISION  KBM(13)
      DOUBLE PRECISION  CIB(3)
      DOUBLE PRECISION  X0,X1,X2,FX0,FX1,FX2
      INTEGER  IFAIL
C
      IF (N.EQ.2) THEN
         X0=KLA(1)
         X1=KLA(2)
         FX0=KBM(1)
         FX1=KBM(2)
         IF (FX0.LT.FX1) THEN
            IFAIL=-2
            RETURN
         ENDIF
         CIB(1)=FX0
         CIB(2)=(FX1-FX0)/(X1-X0)
         RETURN
      ENDIF
C case N=3: parabolic interpolation
      IF (N.GE.3) THEN
         X0=KLA(1)
         X1=KLA(2)
         X2=KLA(3)
         FX0=KBM(1)
         FX1=KBM(2)
         FX2=KBM(3)
         IF (FX0.LT.DMIN1(FX1,FX2)) THEN
            IFAIL=-3
            RETURN
         ENDIF
         CIB(1)=FX0
         CIB(2)=(FX1-FX0)/(X1-X0)
         CIB(3)=((FX2-FX1)/(X2-X1)-CIB(2))/(X2-X0)
         RETURN
      ENDIF
      END



      SUBROUTINE NEWTON2(N,KLA,CIB,X,SDIF)
C+
C  using coeffs CIB from NEWTON1
C  output value + derivatives to SDIF array
C-
      INTEGER           N
      DOUBLE PRECISION  KLA(13)
      DOUBLE PRECISION  CIB(3)
      DOUBLE PRECISION  SDIF(4)
      DOUBLE PRECISION  X,X0,X1
C
      IF (N.EQ.2) THEN
         X0=KLA(1)
         SDIF(1)=CIB(1)+(X-X0)*CIB(2)
         SDIF(2)=CIB(2)
         SDIF(3)=0.
         SDIF(4)=0.
         RETURN
      ENDIF
      IF (N.GE.3) THEN
         X0=KLA(1)
         X1=KLA(2)
         SDIF(1)=CIB(1)+(X-X0)*(CIB(2)+(X-X1)*CIB(3))
         SDIF(2)=CIB(2)+(X*2.-X0-X1)*CIB(3)
         SDIF(3)=CIB(3)*2.
         SDIF(4)=0.
         RETURN
      ENDIF
      END



      SUBROUTINE PREPARE(NB,DL,BM,KLA,KBM)
C
      IMPLICIT   NONE
      INTEGER    NB
      REAL       DL(1)
      REAL       BM(1)
      DOUBLE PRECISION KLA(13)
      DOUBLE PRECISION KBM(13)
C
      INTEGER    I
C
      DO I=1,NB
         KBM(I)=BM(I)
         KLA(I)=DL(I)
      ENDDO
      END



      SUBROUTINE FAIL(IFAIL)
      IMPLICIT      NONE
      INTEGER       IFAIL
C
      INTEGER       ISTAT
      CHARACTER*80  STRING
C
  100 FORMAT('*** WARNING: IPOL--->IFAIL=',I2)

      WRITE(STRING,100) IFAIL
      CALL STTPUT(STRING,ISTAT)
      END
