C @(#)rfotctrans.for	19.1 (ES0-DMD) 02/25/03 13:30:14
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
