C @(#)mtablv.for 17.1.1.1 (ESO-DMD) 01/25/02 17:14:43 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