C @(#)rfotctrans.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18:15 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 PROGRAM CTRANS C+++ C.IDENTIFICATION: RFOTCTRANS C.PURPOSE: Create two tables for coordinate transformation between two frames C.AUTHOR: R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola C Osservatorio Astronomico di Roma, and C. R.H. Warmels, ESO-IPG C.VERSION: 881119 RB First installation in MIDAS C 890914 RHW Rewritten for MIDAS tables C---- IMPLICIT NONE INCLUDE 'MID_REL_INCL:RFOTDECL.INC' C INTEGER LC INTEGER LT PARAMETER (LC=10000) PARAMETER (LT=41) C INTEGER LF(LC) INTEGER MADRID(1) INTEGER ICOL(12) INTEGER ICIDN, IDN1, IDN2 INTEGER IAV, KUN, KNUL INTEGER ISTAT INTEGER EC, ED, EL INTEGER IGP INTEGER KK INTEGER NCP, NHL INTEGER I INTEGER J INTEGER NST, NPA INTEGER LL, LLK, LF1 INTEGER IROW INTEGER NCINT,NRINT,NSINT,NCAINT,NRAINT INTEGER NCCB1,NRCB1,NSCB1,NCACB1,NRACB1 INTEGER NCCB2,NRCB2,NSCB2,NCACB2,NRACB2 INTEGER TIDINT,TIDCB1,TIDCB2 C DOUBLE PRECISION X(LC),Y(LC),Z(LC),ZX(LC),ZY(LC),XD,YD DOUBLE PRECISION C(LT),CX(LT),CY(LT) REAL TAB1(12),TAB2(12) REAL DEX,DEY,DEX2,DEY2 REAL BIP REAL VZ REAL SQST, SQM C CHARACTER STRING*80 CHARACTER YH*1 CHARACTER TABINT*60,TABCB1*60,TABCB2*60 C LOGICAL NUL(12) C INCLUDE 'MID_INCLUDE:TABLES.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:TABLED.INC' C DATA LF/LC*0/ DATA ICOL/2,3,4,5,6,7,8,9,10,11,12,13/ DATA ICIDN/1/ C 901 FORMAT(' Seq IDENT Difference') 902 FORMAT(' ',I5,2X,I5,2X,E12.5) 903 FORMAT('*** INFO: SQM = ',E12.5) C C *** Is MIDAS out there CALL STSPRO('CTRANS') C C *** start input procedures intermediate file CALL STKRDC('IN_A',1,1,60,IAV,TABINT,KUN,KNUL,ISTAT) !intermediate file CALL STECNT('GET',EC,ED,EL) CALL STECNT('PUT',1,0,0) CALL TBTOPN(TABINT,F_IO_MODE,TIDINT,ISTAT) !open it IF (ISTAT.NE.0) THEN STRING = '*** FATAL: No intermediate file available, '// 2 'sorry ...' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C CALL TBIGET(TIDINT,NCINT,NRINT,NSINT,NCAINT,NRAINT,ISTAT) !get table info IF (NRINT.EQ.0 .OR. ISTAT.NE.0) THEN STRING = '*** FATAL: Something wrong with the intermediate '// 2 'file, sorry ...' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C C *** get coordinate tables CALL STKRDC('INPUTC',1,1,30,IAV,TABCB1,KUN,KNUL,ISTAT) !first table CALL TBTOPN(TABCB1,F_I_MODE,TIDCB1,ISTAT) !open it IF (ISTAT.NE.0) THEN STRING = '*** FATAL: First coordinate table not available, '// 2 'sorry ...' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C CALL TBIGET(TIDCB1,NCCB1,NRCB1,NSCB1,NCACB1,NRACB1,ISTAT) !get table info IF (NRCB1.EQ.0 .OR. ISTAT.NE.0) THEN STRING = '*** FATAL: Something wrong with the first '// 2 'coordinate table, sorry ...' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C CALL STKRDC('INPUTC',1,31,30,IAV,TABCB2,KUN,KNUL,ISTAT) !second table CALL TBTOPN(TABCB2,F_I_MODE,TIDCB2,ISTAT) !open it IF (ISTAT.NE.0) THEN STRING = '*** FATAL: Second coordinate table not available, '// 2 'sorry ...' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C CALL TBIGET(TIDCB2,NCCB2,NRCB2,NSCB2,NCACB2,NRACB2,ISTAT) !get table info IF (NRCB2.EQ.0 .OR. ISTAT.NE.0) THEN STRING = '*** FATAL: Something wrong with the seond '// 2 'coordinate table, sorry ...' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF CALL STECNT('PUT',EC,ED,EL) C C *** get the degrees of the fit CALL STKRDI('INPUTI',1,1,IAV,IGP,KUN,KNUL,ISTAT) C C *** get the data from the table *************************************** DEX = 0. DEY = 0. DEX2 = 0. DEY2 = 0. KK = 0 C DO 10 I=1,NRCB1 CALL TBERDI(TIDCB1,I,ICIDN,IDN1,NUL,ISTAT) CALL TBRRDR(TIDCB1,I,12,ICOL,TAB1,NUL,ISTAT) J = 0 LF1 = 0 C 101 CONTINUE IF (LF1.NE.0 .OR. J.GE.NRCB2) GO TO 102 LF1 = 0 J = J+1 CALL TBERDI(TIDCB2,J,ICIDN,IDN2,NUL,ISTAT) CALL TBRRDR(TIDCB2,J,12,ICOL,TAB2,NUL,ISTAT) IF (IDN1.EQ.IDN2) THEN LF1 = 1 ENDIF GO TO 101 C 102 CONTINUE IF (LF1.EQ.1) THEN LF(I) = 1 KK = KK+1 IF (IGP.GT.0) THEN X(KK) = TAB1(1) Y(KK) = TAB1(2) ZX(KK) = TAB2(1) ZY(KK) = TAB2(2) ELSE DEX = DEX+TAB2(1)-TAB1(1) DEY = DEY+TAB2(2)-TAB1(2) DEX2 = DEX2+(TAB2(1)-TAB1(1))**2 DEY2 = DEY2+(TAB2(2)-TAB1(2))**2 END IF END IF 10 CONTINUE C IF (IGP.EQ.0) THEN DEX = DEX/KK DEY = DEY/KK END IF NST = KK C C *** do the fitting ***************************************************** IF (IGP.GT.0) THEN DO 201 I = 1,NST Z(I) = ZX(I) 201 CONTINUE C CALL IPERD(X,Y,Z,C,NST,IGP,SQM,NPA) WRITE(STRING,901) CALL STTPUT(STRING,ISTAT) LLK = 0 DO 202 LL=1,NST IF (LF(LL).EQ.1) THEN LLK = LLK+1 VZ = BIP(X(LLK),Y(LLK),C,IGP) SQST = VZ-Z(LLK) CALL TBERDI(TIDCB1,LL,ICIDN,IDN1,NUL,ISTAT) CALL TBRRDR(TIDCB1,LL,12,ICOL,TAB1,NUL,ISTAT) WRITE(STRING,902) LL,IDN1,SQST CALL STTPUT(STRING,ISTAT) END IF 202 CONTINUE C DO 203 I=1,NPA CX(I) = C(I) 203 CONTINUE C CALL STTPUT(' ',ISTAT) C DO 204 I = 1,NST Z(I) = ZY(I) 204 CONTINUE C C *** inizio proc. calcola scarto CALL IPERD(X,Y,Z,C,NST,IGP,SQM,NPA) WRITE(STRING,901) LLK = 0 DO 301 LL=1,NST IF (LF(LL).EQ.1) THEN LLK = LLK+1 VZ = BIP(X(LLK),Y(LLK),C,IGP) SQST = VZ-Z(LLK) CALL TBERDI(TIDCB1,LL,ICIDN,IDN1,NUL,ISTAT) CALL TBRRDR(TIDCB1,LL,12,ICOL,TAB1,NUL,ISTAT) WRITE(STRING,902) LL,IDN1,SQST CALL STTPUT(STRING,ISTAT) END IF 301 CONTINUE C DO 302 I = 1,NPA CY(I) = C(I) 302 CONTINUE ELSE SQM = SQRT(DEX2/KK-DEX**2) WRITE(STRING,903) SQM CALL STTPUT(STRING,ISTAT) SQM = SQRT(DEY2/KK-DEY**2) WRITE(STRING,903) SQM CALL STTPUT(STRING,ISTAT) END IF C C *** do we want this or not CALL STTPUT('*** Do you want this transformation [N]?',ISTAT) CALL STKPRC('*** WARNING: If "Y" the intermediate table will '// 2 'be modified: ','INPUTC',1,1,1,IAV,YH, 3 KUN,KNUL,ISTAT) CALL UPCAS(YH,YH) IF (YH.EQ.'Y') THEN C C *** read table descriptor IROW = 1 401 CONTINUE CALL INTWRD(TIDINT,IROW,NCP,NHL) IF (IGP.GT.0) THEN XD = PARINT(1) YD = PARINT(2) VZ = BIP(XD,YD,CX,IGP) PARINT(1) = VZ+.5 VZ = BIP(XD,YD,CY,IGP) PARINT(2) = VZ+.5 ELSE PARINT(1) = PARINT(1)+DEX PARINT(2) = PARINT(2)+DEY END IF CALL INTWWR(TIDINT,IROW,NCP,NHL) IROW = IROW + NCP + NHL IF (IROW.LE.NRINT) THEN GO TO 401 ENDIF ELSE STRING = '*** INFO: User veto: no conversion done' CALL STTPUT(STRING,ISTAT) ENDIF C CALL STSEPI END