C @(#)tdregt.for 17.1.1.1 (ES0-DMD) 01/25/02 17:47:16 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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 16:14 - 19 NOV 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.IDENTIFICATION C C program TDREGT C C.PURPOSE C C POLYNOMIAL FIT OF ONE OR TWO VARIABLES BETWEEN TWO TABLES C THE PROGRAM CAN BE USED TO FIND CORRESPONDENCES BETWEEN C ENTRIES IN THE TWO TABLES C C REGRESSION/TABLE TABLE X1[,X2] CATAL Y1[,Y2] DEGREE TOL [COEFFS] C C THE METHOD IS NOT YET IMPLEMENTED FOR TWO VARIABLES C C.KEYWORDS C C TABLES, REGRESSION, IDENTIFICATIONS C C.ALGORITHM C C THERE ARE TWO MODES OF OPERATION : C 1. INITIAL GUESS FOR THE CORRESPONDENCE IS DEFINED BY THE C REFERENCE COLUMNS OF BOTH TABLES C IN THIS CASE REFERNCE COLUMN IN SECOND OUT HAS TO BE SORTED C 2. INITIAL GUESS IS DEFINED BY A SET OF COEFFS. WITH THE C GEOMETRIC TRANSFORMATION C THE RESULT IS C 1. THE SET OF COEFFS. OF THE TRANSFORMATION C X1 <- Y1 OR C X1 <- Y1,Y2 C C 2. THE CORRESPONDENCE BETWEEN THE TWO TABLES IS DEFINED BY C UPDATING THE REFERENCE COLUMN OF THE FIRST INPUT TABLE C C 3. If two columns of the problem table are specified, the C second column defines the iteration variable. C C C----------------------------------------------------------- C SUBROUTINE TDREGT IMPLICIT NONE C C ... PARAMETER DEFINITION C CHARACTER LINE1*64,KEY*64 CHARACTER*17 COLX1,COLX2,COLY1,COLY2 CHARACTER*8 FORM CHARACTER*8 FORM1,FORM2 CHARACTER*16 LABEL1,LABEL2,UNIT1,UNIT2 CHARACTER*64 TABLE,CATAL CHARACTER*80 OUTPUT1, OUTPUT2, OUTPUT3 C INTEGER MADRID INTEGER DTYPE1, DTYPE2, DTYPE, DTYPR1, DTYPR2 INTEGER IDEG(2),ICOL1(3),ICOL2(3),ICOL3(2) INTEGER NPIX(2),INDX1(400),INDX2(400) INTEGER IACT,STATUS,INDX,NINC,NINC1,NINC2 INTEGER NOUTC, TID, TID1, LEN, NACT, NRCOL, IST, I, J INTEGER NSTART, NORD, NCOL2, NROW2, NSC2, NAC2, NAR2 INTEGER*8 IBUF INTEGER NCOL1, NROW1, NSC1, NAC, NAR, NBYTES, IS INTEGER INDEX, KUN, KNUL C REAL TOL REAL START(2) REAL STEP(2) C DOUBLE PRECISION ERROR,COEFF(30),COEFFR(30) C INCLUDE 'MID_INCLUDE:TABLES.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:TABLED.INC' C DATA STEP/1.,1./ DATA START/0.,0./ DATA OUTPUT1 ./' scan total unid. ident. abs(res) '/ DATA OUTPUT2 ./' no. iterat. lines lines lines >3*std.dev. stdev.'/ DATA OUTPUT3 ./' ---- ------ ----- ----- ----- ----------- ------'/ C C ... GET INPUT PARAMETERS + DEFAULT C CALL STKRDC('P1',1,1,64,IACT,TABLE,KUN,KNUL,STATUS) CALL STKRDC('P2',1,1,64,IACT,LINE1,KUN,KNUL,STATUS) INDX = INDEX(LINE1,',') IF (INDX.EQ.0) THEN NINC = 1 COLX1 = LINE1 ELSE NINC = 2 COLX1 = LINE1(1:INDX-1) COLX2 = LINE1(INDX+1:) END IF CALL STKRDC('P3',1,1,64,IACT,CATAL,KUN,KNUL,STATUS) CALL STKRDC('P4',1,1,64,IACT,LINE1,KUN,KNUL,STATUS) INDX = INDEX(LINE1,',') IF (INDX.EQ.0) THEN NOUTC = 1 COLY1 = LINE1 ELSE NOUTC = 2 COLY1 = LINE1(1:INDX-1) COLY2 = LINE1(INDX+1:) CALL STTPUT(' Not yet implemented ',STATUS) GO TO 50 END IF CALL STKRDI('INPUTI',1,2,IACT,IDEG,KUN,KNUL,STATUS) CALL STKRDR('INPUTR',1,1,IACT,TOL,KUN,KNUL,STATUS) CALL STKRDC('P7',1,1,64,IACT,KEY,KUN,KNUL,STATUS) C C ... INIT INPUT TABLE C CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS) CALL TBIGET(TID,NCOL1,NROW1,NSC1,NAC,NAR,STATUS) CALL TBKGET(TID,ICOL1(1),STATUS) IF (ICOL1(1).LE.0) THEN CALL STTPUT(' There is no REFERENCE column in first table', + STATUS) GO TO 50 END IF CALL TBFGET(TID,ICOL1(1),FORM,LEN,DTYPR1,STATUS) CALL TBCSER(TID,COLX1,ICOL1(2),STATUS) CALL TBFGET(TID,ICOL1(2),FORM,LEN,DTYPE,STATUS) IF (DTYPE.EQ.D_C_FORMAT) THEN CALL STTPUT(' Illegal column type in first table',STATUS) GO TO 50 END IF IF (NINC.EQ.2) THEN CALL TBCSER(TID,COLX2,ICOL1(3),STATUS) CALL TBFGET(TID,ICOL1(3),FORM,LEN,DTYPE,STATUS) IF (DTYPE.EQ.D_C_FORMAT) THEN CALL STTPUT(' Illegal column type in first table',STATUS) GO TO 50 END IF END IF IF (ICOL1(2).EQ.-1 .OR. ICOL1(3).EQ.-1) THEN CALL STTPUT(' Column in first table not found',STATUS) GO TO 50 END IF C C ... INIT REFERENCE TABLE C CALL TBTOPN(CATAL,F_I_MODE,TID1,STATUS) CALL TBIGET(TID1,NCOL2,NROW2,NSC2,NAC2,NAR2,STATUS) CALL TBKGET(TID1,ICOL2(1),STATUS) IF (ICOL2(1).LE.0) THEN CALL STTPUT(' There is no REFERENCE column in second table', + STATUS) GO TO 50 END IF IF (NSC2 .NE. ICOL2(1)) THEN CALL STTPUT(' Second table is not sorted by REFERENCE column', + STATUS) GO TO 50 END IF CALL TBFGET(TID1,ICOL2(1),FORM,LEN,DTYPR2,STATUS) IF (DTYPR1.NE.DTYPR2) THEN CALL STTPUT(' Warning: columns with different types',STATUS) C GO TO 50 END IF CALL TBCSER(TID1,COLY1,ICOL2(2),STATUS) CALL TBFGET(TID1,ICOL2(2),FORM1,LEN,DTYPE1,STATUS) IF (DTYPE1.EQ.D_C_FORMAT) THEN CALL STTPUT(' Illegal column type in second table',STATUS) GO TO 50 END IF IF (NOUTC.EQ.2) THEN CALL TBCSER(TID1,COLY2,ICOL2(3),STATUS) CALL TBFGET(TID1,ICOL2(3),FORM2,LEN,DTYPE2,STATUS) IF (DTYPE2.EQ.D_C_FORMAT) THEN CALL STTPUT(' Illegal column type in second table',STATUS) GO TO 50 END IF END IF IF (ICOL2(2).EQ.-1 .OR. ICOL2(3).EQ.-1) THEN CALL STTPUT(' Column in second table not found',STATUS) GO TO 50 END IF C C ... IDENTIFY FEATURES(iteration on equal order numbers ) C NINC1 = NINC + 1 NINC2 = NOUTC + 1 NPIX(1) = NROW1 NPIX(2) = 6*2 C C ... PREPARE WORK SPACE C NBYTES = NPIX(1)*NPIX(2)*4 CALL TDMGET(NBYTES,IBUF,IS) CALL SPCOP1(NROW1,TID,NINC1,ICOL1,TID1,NINC2,ICOL2, + MADRID(IBUF),NACT,NORD,INDX1,INDX2,NSTART) C C ... HEADER C CALL STTPUT(OUTPUT1, STATUS) CALL STTPUT(OUTPUT2, STATUS) CALL STTPUT(OUTPUT3, STATUS) C C ... IDENTIFICATION OF MAIN REFERENCE ORDER C NRCOL = NINC2 - 1 ERROR = TOL CALL IDENTB(NROW1,MADRID(IBUF),INDX1(NSTART),INDX2(NSTART), + TID1,NROW2,NRCOL,ICOL2(2),ERROR,COEFF,IDEG,NSTART,IST) DO 10 I = 1,30 COEFFR(I) = COEFF(I) 10 CONTINUE C C ... IDENTIFICATION OF PREVIOUS ORDERS C IF (NSTART.NE.1) THEN DO 20 I = NSTART - 1,1,-1 DO 15 J = 1,30 COEFF(J) = COEFFR(J) 15 CONTINUE CALL IDENTA(NROW1,MADRID(IBUF),INDX1(I),INDX2(I), + TID1,NROW2,NRCOL,ICOL2(2),ERROR,COEFF,IDEG,I,IST) 20 CONTINUE END IF C C ... IDENTIFICATION OF NEXT ORDERS C IF (NSTART.NE.NORD) THEN DO 40 I = NSTART + 1,NORD DO 30 J = 1,30 COEFF(J) = COEFFR(J) 30 CONTINUE CALL IDENTA(NROW1,MADRID(IBUF),INDX1(I),INDX2(I), + TID1,NROW2,NRCOL,ICOL2(2),ERROR,COEFF,IDEG,I,IST) 40 CONTINUE END IF C C ... COPY FINAL RESULTS WITH THE IDENTIFICATIONS IN OUT TABLE C CALL TBCSER(TID,COLY1,ICOL3(1),STATUS) IF (ICOL3(1).LE.0) THEN CALL TBFGET(TID1,ICOL2(2),FORM1,LEN,DTYPE1,STATUS) CALL TBUGET(TID1,ICOL2(2),UNIT1,STATUS) CALL TBLGET(TID1,ICOL2(2),LABEL1,STATUS) CALL TBCINI(TID,DTYPE1,1,FORM1,UNIT1,LABEL1,ICOL3(1),STATUS) ENDIF IF (NOUTC.EQ.2) THEN CALL TBCSER(TID,COLY2,ICOL3(2),STATUS) IF (ICOL3(1).LE.0) CALL TBCINI(TID,DTYPE2,1,FORM2,UNIT2, + LABEL2,ICOL3(2),STATUS) END IF CALL SPCOP2(NROW1,MADRID(IBUF),TID,NOUTC,ICOL3) C C ... end C 50 CALL TBTCLO(TID,STATUS) CALL TBTCLO(TID1,STATUS) CALL TDMFRE(NBYTES,IBUF,IS) RETURN END SUBROUTINE IDENTA(NROW1,BUF,I1,I2,TID,NROW2,NC,IC,ERROR, + COEFF,IDEG,INDX,IST) C C IDENTIFY DETECTED LINES C INITIAL GUESSES DEFINED BY THE COEFFS C INPUT STORED IN BUF WITH THE FOLLOWING STRUCTURE C BUF(I,1) - X1 (SAMPLE) C BUF(I,2) - X2 (ORDER) C BUF(I,3) - Y1 (IDENTIFIED WAVELENGTH) C BUF(I,4) - Y2 (SECOND VARIABLE) C BUF(I,5) - LINE WEIGHT AS 1 MANUAL IDENTIFICATION C 0 NOT IDENTIFIED C -1 AUTOMATIC IDENTIFICATION C BUF(I,6) - RESIDUALS C IMPLICIT NONE C INTEGER NROW1, I1, I2, NROW2, NC, IST, NSTEP, NV, NTOT INTEGER NTOTAL, NTOTU, NTOTI, NTOTR1,NTOTR2, I, TID, ISTAT INTEGER IC(2),IDEG(2),INDX C DOUBLE PRECISION ERROR,BUF(NROW1,6),COEFF(30) DOUBLE PRECISION STDEV,STDEV1,STDEV3 C CHARACTER*80 OUT C STDEV1 = -1.D20 IST = 0 C C ITERATION C DO 20 NSTEP = 1,10 NV = (IDEG(1)+1)* (IDEG(2)+1) C C BUILD L - S MATRIX AND SOLVE SYSTEM C IF (NSTEP.GT.1) CALL LSOLVE(NROW1,BUF,I1,I2,IDEG,NV,COEFF) C C IDENTIFY FEATURES C CALL LSIDEN(IDEG,NV,COEFF,NROW1,BUF,I1,I2,TID,NROW2,NC, + IC,ERROR,NTOT,STDEV) C C SELECT LINES WITH LOW RESIDUALS C STDEV3 = STDEV*3.D0 IF (NTOT.EQ.0) THEN OUT = ' ' CALL STTPUT(' Wrong identifications in input table',ISTAT) IST = 1 RETURN END IF C C STATISTICS ON THE IDENTIFICATIONS C ! TOTAL NUMBER OF LINES NTOTAL = 0 ! TOTAL NUMBER OF UNIDENTIFIED LINES NTOTU = 0 ! TOTAL NUMBER OF IDENTIFIED LINES NTOTI = 0 ! TOTAL NUMBER OF LINES WITH RESIDUALS .LT. 0.1 PIXEL NTOTR1 = 0 ! TOTAL NUMBER OF LINES WITH RESIDUALS .GT.3*STDEV NTOTR2 = 0 DO 10 I = I1,I2 NTOTAL = NTOTAL + 1 IF (DABS(BUF(I,5)).GT.0.5D0) THEN NTOTI = NTOTI + 1 IF (DABS(BUF(I,6)).GT.STDEV3) THEN BUF(I,4) = 0.0D0 BUF(I,5) = 0.0D0 BUF(I,6) = 0.0D0 NTOTR2 = NTOTR2 + 1 END IF END IF 10 CONTINUE NTOTU = NTOTAL - NTOTI IF (NTOTI.LT.NV) THEN CALL STTPUT(' Error : not enough identified entries', + ISTAT) IST = 1 RETURN END IF IF (DABS((STDEV-STDEV1)/STDEV).LE.0.005D0) GOTO 100 STDEV1 = STDEV 20 CONTINUE RETURN 100 WRITE(OUT,300) INDX,NSTEP,NTOTAL,NTOTU,NTOTI,NTOTR2,STDEV CALL STTPUT(OUT,ISTAT) RETURN 300 FORMAT(1X,I5,I6,2X,I5,2X,I5,2X,I6,2X,I8,5X,E10.3) END SUBROUTINE IDENTB(NROW1,BUF,I1,I2,TID,NROW2,NC,IC,ERROR, + COEFF,IDEG,INDX,IST) C C IDENTIFY DETECTED LINES C INITIAL GUESSES DEFINED BY THE IDENTIFICATIONS IN BUF(I,3) C INPUT STORED IN BUF WITH THE FOLLOWING STRUCTURE C BUF(I,1) - X1 (SAMPLE) C BUF(I,2) - X2 (ORDER) C BUF(I,3) - Y1 (IDENTIFIED WAVELENGTH) C BUF(I,4) - Y2 (SECOND VARIABLE) C BUF(I,5) - LINE WEIGHT AS 1 MANUAL IDENTIFICATION C 0 NOT IDENTIFIED C -1 AUTOMATIC IDENTIFICATION C BUF(I,6) - RESIDUALS C IMPLICIT NONE C INTEGER NROW1, I1, I2, NROW2, NC, IST, NSTEP, NV, NTOT INTEGER NTOTAL, NTOTU, NTOTI, NTOTR1,NTOTR2, I, TID INTEGER ISTAT,INDX INTEGER IC(2),IDEG(2) C DOUBLE PRECISION ERROR,BUF(NROW1,6),COEFF(30) DOUBLE PRECISION STDEV,STDEV1,STDEV3 C CHARACTER*80 OUT C STDEV1 = -1.D20 IST = 0 C C ITERATION C DO 20 NSTEP = 1,10 NV = (IDEG(1)+1)* (IDEG(2)+1) C C BUILD L - S MATRIX AND SOLVE SYSTEM C CALL LSOLVE(NROW1,BUF,I1,I2,IDEG,NV,COEFF) C C IDENTIFY FEATURES C CALL LSIDEN(IDEG,NV,COEFF,NROW1,BUF,I1,I2,TID,NROW2,NC, + IC,ERROR,NTOT,STDEV) C C SELECT LINES WITH LOW RESIDUALS C STDEV3 = STDEV*3.D0 IF (NTOT.EQ.0) THEN OUT = ' ' CALL STTPUT(' Wrong identifications in input table',ISTAT) IST = 1 RETURN END IF C C STATISTICS ON THE IDENTIFICATIONS C ! TOTAL NUMBER OF LINES NTOTAL = 0 ! TOTAL NUMBER OF UNIDENTIFIED LINES NTOTU = 0 ! TOTAL NUMBER OF IDENTIFIED LINES NTOTI = 0 ! TOTAL NUMBER OF LINES WITH RESIDUALS .LT. 0.1 PIXEL NTOTR1 = 0 ! TOTAL NUMBER OF LINES WITH RESIDUALS .GT.3*STDEV NTOTR2 = 0 DO 10 I = I1,I2 NTOTAL = NTOTAL + 1 IF (DABS(BUF(I,5)).GT.0.5D0) THEN NTOTI = NTOTI + 1 IF (DABS(BUF(I,6)).GT.STDEV3) THEN BUF(I,4) = 0.0D0 BUF(I,5) = 0.0D0 BUF(I,6) = 0.0D0 NTOTR2 = NTOTR2 + 1 END IF END IF 10 CONTINUE NTOTU = NTOTAL - NTOTI IF (NTOTI.LT.NV) THEN CALL STTPUT(' Error : not enough identified entries', + ISTAT) IST = 1 RETURN END IF IF (DABS((STDEV-STDEV1)/STDEV).LE.0.005D0) GOTO 100 STDEV1 = STDEV 20 CONTINUE RETURN 100 WRITE(OUT,300) INDX,NSTEP,NTOTAL,NTOTU,NTOTI,NTOTR2,STDEV CALL STTPUT(OUT,ISTAT) RETURN 300 FORMAT(1X,I5,I6,2X,I5,2X,I5,2X,I6,2X,I8,5X,E10.3) END SUBROUTINE LSIDEN(IDEG,NV,COEFF,NROW1,BUF,I1,I2,TID,NROW2, + NC,IC,ERROR,NTOT,STDEV) C C IDENTIFY DETECTED LINES USING DISPERSION COEFFS C IMPLICIT NONE C INTEGER NTOT, NV, NROW1, I1, I2, TID, NROW2, NC INTEGER K, L, NACT, I, IPOINT, STATUS, ISTART INTEGER IC(NC),IDEG(2) C LOGICAL NULL(2) C DOUBLE PRECISION COEFF(NV),BUF(NROW1,6),ERROR,SAMPLE DOUBLE PRECISION DD(2),STDEV, POLEV C DATA DD/0.D0,0.D0/ DATA NULL/.FALSE.,.FALSE./ C C LOOP ON THE LINES IN THE CATAL C NTOT = 0 STDEV = 0.D0 C K = IDEG(1) + 1 C L = IDEG(2) + 1 K = IDEG(1) L = IDEG(2) NACT = I2 - I1 + 1 ISTART = 1 DO 10 I = 1,NROW2 C C ... USE ONLY SUITABLE RANGE !! C CALL TBRRDD(TID,I,NC,IC,DD,NULL,STATUS) IF ( .NOT. NULL(1) .AND. .NOT. NULL(2)) THEN SAMPLE = POLEV(K,L,NV,COEFF,DD) C CALL SEABAD(BUF(I1,1),NACT,SAMPLE,ERROR,ISTART,IPOINT) C C ... LINE IDENTIFIED C IF (IPOINT.GT.0) THEN IPOINT = IPOINT + I1 - 1 BUF(IPOINT,3) = DD(1) BUF(IPOINT,4) = DD(1) BUF(IPOINT,5) = -1.D0 BUF(IPOINT,6) = BUF(IPOINT,1) - SAMPLE STDEV = STDEV + BUF(IPOINT,6)*BUF(IPOINT,6) NTOT = NTOT + 1 END IF END IF 10 CONTINUE IF (NTOT.GT.0) STDEV = DSQRT(STDEV/NTOT) RETURN END C SUBROUTINE LSOLVE(NROW1,BUF,I1,I2,IDEG,NV,COEFF) C C SOLVE THE L-S FIT TO THE DATA IMPLICIT NONE C INTEGER NROW1, I1, I2, NV, K, L, NTOTAL, LL, N, I INTEGER KMIN, N1, K1, ISTAT, LL1 INTEGER IDEG(2) C DOUBLE PRECISION BUF(NROW1,6),COEFF(NV) DOUBLE PRECISION G(20,20) C C BUILD INPUT ARRAYS C K = IDEG(1) L = IDEG(2) NTOTAL = 0 LL = 0 N = NV N1 = N + 1 DO 20 I = I1,I2 IF (DABS(BUF(I,5)).GT.0.5D0) THEN NTOTAL = NTOTAL + 1 LL1 = LL+1 CALL TDSET2(LL1,BUF(I,3),BUF(I,4),BUF(I,1),K,L,G,COEFF,N,20) IF (LL.NE.0) THEN KMIN = MIN(LL,N1) DO 10 K1 = 1,KMIN CALL TDHHTR(K1,LL1,G,COEFF,N,20) 10 CONTINUE END IF LL = MIN(LL1,N1) END IF 20 CONTINUE C C SOLVE THE LINEAR SYSTEM C IF (NTOTAL.LT.NV) THEN CALL STTPUT(' Not enough identified features',ISTAT) RETURN END IF CALL TDSOLV(G,COEFF,N,20) RETURN END DOUBLE PRECISION FUNCTION POLEV(K,L,NV,DPAR,VAL) C C COMPUTE FUNCTION VALUE C IMPLICIT NONE C INTEGER K, L, NV, IP, L1, K1 C DOUBLE PRECISION DPAR(NV),VAL(2) DOUBLE PRECISION DC,RESULT,WORK(30) C IP = 0 DC = 1.D0 RESULT = 0.D0 DO 20 L1 = 0, L IP = IP + 1 WORK(IP) = DC RESULT = RESULT + WORK(IP)*DPAR(IP) DO 10 K1 = 1,K IP = IP + 1 WORK(IP) = WORK(IP-1)*VAL(1) RESULT = RESULT + WORK(IP)*DPAR(IP) 10 CONTINUE DC = DC*VAL(2) 20 CONTINUE POLEV = RESULT RETURN END SUBROUTINE SPCOP1(N,TID,NC1,IC1,TID1,NC2,IC2,BUF,NACT, + NORD,INDX1,INDX2,NSTART) C C PREPARE WORK SPACE C IMPLICIT NONE C INTEGER NACT, TID, TID1, N, NC1, NC2, NORD INTEGER NSTART, NC, LEN, STATUS, NIDENT, NMAX, I, I1 INTEGER IC1(NC1),IC2(NC2) INTEGER INDX1(1),INDX2(1) INTEGER DTYPE, NEXT C LOGICAL SEL,NULL1(3),NULL2(3), NULL3 C REAL RVAL1,RVAL2,RERR C DOUBLE PRECISION BUF(N,6) DOUBLE PRECISION DERR DOUBLE PRECISION VAL1(3) DOUBLE PRECISION VAL2(3) DOUBLE PRECISION REF C CHARACTER FORM*10 C INCLUDE 'MID_INCLUDE:TABLES.INC' INCLUDE 'MID_INCLUDE:TABLED.INC' C DATA VAL1/0.D0,0.D0,0.D0/ DATA VAL2/0.D0,0.D0,0.D0/ DATA REF/1.E30/ DATA NULL1/.FALSE.,.FALSE.,.FALSE./ DATA NULL2/.FALSE.,.FALSE.,.FALSE./ C I1 = 1 NACT = 0 DERR = .1D0 NC = NC2 - 1 NORD = 0 C C ... CHECK FORMAT IN CATAL C CALL TBFGET(TID1,IC2(1),FORM,LEN,DTYPE,STATUS) IF (DTYPE.EQ.D_R4_FORMAT) RERR = DERR C NSTART = 1 NIDENT = 0 NMAX = 0 DO 10 I = 1,N CALL TBSGET(TID,I,SEL,STATUS) IF (SEL) THEN CALL TBRRDD(TID,I,NC1,IC1,VAL1,NULL1,STATUS) IF ( .NOT. NULL1(2) .AND. .NOT. NULL1(3)) THEN NACT = NACT + 1 ! ASSIGN X1 BUF(NACT,1) = VAL1(2) ! ASSIGN X2 (OPTIONAL) BUF(NACT,2) = VAL1(3) ! NEW ORDER IF (VAL1(3).NE.REF) THEN REF = VAL1(3) NORD = NORD + 1 INDX1(NORD) = NACT ! CHECK NO. OF IDENTIFICATIONS IF (NIDENT.GT.NMAX) THEN NMAX = NIDENT NSTART = NORD END IF NIDENT = 0 IF (NORD.GT.1) INDX2(NORD-1) = NACT - 1 END IF ! ASSIGN Y1 BUF(NACT,3) = 0.D0 ! ASSIGN Y2 (OPTIONAL) BUF(NACT,4) = 0.D0 ! FLAG AS NON-EXISTING BUF(NACT,5) = 0.D0 ! RESIDUAL (X1-FIT(Y1,Y2)) BUF(NACT,6) = 0.D0 IF ( .NOT. NULL1(1)) THEN NIDENT = NIDENT + 1 IF (DTYPE.EQ.D_R4_FORMAT) THEN RVAL1 = VAL1(1) CALL TBESRR(TID1,IC2,RVAL1,RERR,I1,NEXT, + STATUS) CALL TBERDR(TID1,NEXT,IC2,RVAL2, + NULL3,STATUS) VAL2(3) = RVAL2 ELSE CALL TBESRD(TID1,IC2,VAL1(1),DERR,I1,NEXT, + STATUS) CALL TBERDD(TID1,NEXT,IC2,VAL2(3),NULL3, + STATUS) END IF IF (NEXT.GT.0) THEN CALL TBSGET(TID1,NEXT,SEL,STATUS) IF (SEL) THEN CALL TBRRDD(TID1,NEXT,NC,IC2(2), . VAL2,NULL2,STATUS) IF (.NOT.NULL2(1).AND..NOT.NULL2(2)) THEN !ASSIGN Y1 BUF(NACT,3) = VAL2(1) ! ASSIGN Y2 (OPTIONAL) BUF(NACT,4) = VAL2(2) ! FLAG AS EXISTING BUF(NACT,5) = 1.D0 I1 = NEXT ! ASSUMES SORTED TABLES END IF END IF END IF END IF END IF END IF 10 CONTINUE IF (NIDENT.GT.NMAX) NSTART = NORD NSTART = NSTART-1 INDX2(NORD) = NACT RETURN END SUBROUTINE SPCOP2(N,BUF,TID,NC,IC) C C COPY FINAL IDENTIFICATIONS C IMPLICIT NONE C INTEGER N, TID, NC, I, NACT, STATUS INTEGER IC(2) C LOGICAL SEL C DOUBLE PRECISION BUF(N,6),VAL(3) C INCLUDE 'MID_INCLUDE:TABLES.INC' INCLUDE 'MID_INCLUDE:TABLED.INC' C NACT = 0 DO 10 I = 1,N CALL TBSGET(TID,I,SEL,STATUS) IF (SEL) THEN NACT = NACT + 1 IF (DABS(BUF(NACT,5)).GT.0.5D0) THEN VAL(1) = BUF(NACT,3) VAL(2) = BUF(NACT,4) CALL TBRWRD(TID,I,NC,IC,VAL,STATUS) ELSE CALL TBEDEL(TID,I,IC(1),STATUS) IF (NC.EQ.2) CALL TBEDEL(TID,I,IC(2),STATUS) END IF END IF 10 CONTINUE RETURN END SUBROUTINE SEABAD(ARRAY, N, VALUE, ERROR, IFIRST,IPOS) C C BINARY SEARCH OF ELEMENTS IN THE ARRAY C DOUBLE PRECISION VERSION C ASCENDING ORDER C IMPLICIT NONE INTEGER N DOUBLE PRECISION ARRAY(N) DOUBLE PRECISION VALUE DOUBLE PRECISION ERROR INTEGER IFIRST INTEGER IPOS C INTEGER I, LOWER, UPPER C C ... INITIALIZE COUNTERS C LOWER = IFIRST UPPER = N IPOS = 0 C C ... LOOP C 5 CONTINUE I = (LOWER + UPPER)/2 C C ... COMPARE KEYS C IF (DABS(VALUE-ARRAY(I)).LE.ERROR) THEN IPOS = I GO TO 10 ELSE IF (VALUE.LT.ARRAY(I)) THEN UPPER = I - 1 ELSE LOWER = I + 1 ENDIF ENDIF IF (LOWER.LE.UPPER) GO TO 5 C C ... FIND FIRST ENTRY IN A SET OF EQUAL KEYS C 10 IF (IPOS.LE.IFIRST) RETURN IF (DABS(VALUE-ARRAY(IPOS-1)).GT.ERROR) RETURN 20 CONTINUE IPOS = IPOS - 1 IF (IPOS.EQ.IFIRST) RETURN IF (DABS(VALUE-ARRAY(IPOS-1)).LE.ERROR) GOTO 20 RETURN END