C @(#)rfotcbase.for	19.1 (ES0-DMD) 02/25/03 13:30:13
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 CBASE
C+++
C.IDENTIFICATION: RFOTCBASE
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          890913 RHW  Rewritten on the basis on MCL and MIDAS tables
C----
      IMPLICIT       NONE
      INCLUDE        'MID_REL_INCL:RFOTDECL.INC'
C
      INTEGER        NCPAR
      PARAMETER      (NCPAR=12)
      INTEGER        ICBASE(2)
      INTEGER        IVX(9)
      INTEGER        NAXIS1,NPIX1(2)
      INTEGER        NAXIS2,NPIX2(2)
      INTEGER        MADRID(1)
      INTEGER        LX(2),LY(2)
      INTEGER        ED, EC, EL
      INTEGER        ICIDN
      INTEGER*8      IPNR
      INTEGER        IAV, ISTAT
      INTEGER        I, IR, IROW, INX, IT
      INTEGER        J
      INTEGER        KUN, KNUL
      INTEGER        KK, KL
      INTEGER        NCTAB, NRTAB, NSTAB, NACTAB, NARTAB
      INTEGER        NPL0, NL0
      INTEGER        NLIV, NC
      
      INTEGER        IPNTR1,IPNTR2,IMF1,IMF2
      INTEGER        TIDI,TIDO1,TIDO2
      INTEGER        NCOL, NROW
      INTEGER        REGTYP,REGCOL
      INTEGER        ICPAR(NCPAR)
      INTEGER        FIT
C
      DOUBLE PRECISION BEGIN1(2),STEP1(2)
      DOUBLE PRECISION BEGIN2(2),STEP2(2)
C
      REAL           COO(2)
      REAL           FO1(1000,2)
      REAL           H1(1000,2)
      REAL           PAR(4),PARS(4)
      REAL           RX(100),RY(100)
      REAL           TAB(NCPAR)
      REAL           VET(100)
      REAL           X1(1000,2),Y1(1000,2)
      REAL           FG
      REAL           RME, VME
      REAL           SQR, SL, SQM
C
      CHARACTER*60   IDENT1,CUNIT1,IDENT2,CUNIT2
      CHARACTER      TABI*60,TABO1*60,TABO2*60
      CHARACTER*60   IMAG1, IMAG2
      CHARACTER*16   LABEL(NCPAR),REGLAB
      CHARACTER*16   UNIT(NCPAR),REGUNI
      CHARACTER*80   STRING
      CHARACTER*16   REGFOR
      CHARACTER*16   FORMR4,FORMI4
C
      LOGICAL        NUL(2)
C
      INCLUDE        'MID_INCLUDE:ST_DEF.INC'
      COMMON         /VMR/MADRID
      INCLUDE        'MID_INCLUDE:ST_DAT.INC'
C
      DATA           ICIDN/1/
      DATA           ICBASE/1,2/
      DATA           FORMI4/'I6'/
      DATA           FORMR4/'E12.4'/
      DATA           LABEL /'X      ', 'Y      ', 'INT    ', 'LOC_BGD ', 
     2                      'MAG1   ', 'MAG2   ', 'MAG3   ', 'MAG_CNV ', 
     3                      'SIGMA  ', 'BETA   ', 'SIQ    ', 'CHI_SQ  '/
      DATA           UNIT  /'PIXEL  ', 'PIXEL  ', '       ', '        ',
     2                      'MAG.   ', 'MAG.   ', 'MAG.   ', 'MAG.    ',
     3                      '       ', '       ', '       ', '        '/
      DATA           IVX/1,2,3,4,5,6,7,8,9/
C
C *** Is MIDAS out there
      CALL STSPRO('CBASE')
C
C *** get the frame names
      CALL STKRDC('IN_A',1,1,60,IAV,IMAG1,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',1,1,60,IAV,IMAG2,KUN,KNUL,ISTAT)
      CALL STIGET(IMAG1,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,NAXIS1,
     2            NPIX1,BEGIN1,STEP1,IDENT1,CUNIT1,IPNTR1,IMF1,ISTAT)
      CALL STIGET(IMAG2,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,NAXIS2,
     2            NPIX2,BEGIN2,STEP2,IDENT2,CUNIT2,IPNTR2,IMF2,ISTAT)
C
C *** get the input coordinate and output 
      CALL STKRDC('INPUTC',1,1,60,IAV,TABI,KUN,KNUL,ISTAT)
      CALL STKRDC('P3',1,1,60,IAV,TABO1,KUN,KNUL,ISTAT)
      CALL STKRDC('P4',1,1,60,IAV,TABO2,KUN,KNUL,ISTAT)
C
C *** read the input table 
      CALL TBTOPN(TABI,F_I_MODE,TIDI,ISTAT)
      CALL STECNT('GET',EC,ED,EL)
      CALL STECNT('PUT',1,0,0)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: No input coordinate table available'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
      CALL TBIGET(TIDI,NCTAB,NRTAB,NSTAB,NACTAB,NARTAB,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: Problems with input coordinate table'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
      IF (NRTAB.EQ.0) THEN
         STRING = '*** FATAL: No points in coordinate table'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
      CALL STECNT('PUT',EC,ED,EL)
C
C *** create the output tables
      NROW  = NRTAB/2
      NCOL  = NCPAR + 1
C
      CALL TBTINI(TABO1,0,F_O_MODE,NCOL,NROW,TIDO1,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL:  Problems with opening first'//
     2            'transformation table'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
      CALL TBTINI(TABO2,0,F_O_MODE,NCOL,NROW,TIDO2,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL:  Problems with opening second '//
     2            'tranformation table'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
      REGTYP = D_I4_FORMAT
      REGFOR = FORMI4
      REGUNI = '   '
      REGLAB = 'IDENT'
      CALL TBCINI(TIDO1,REGTYP,1,REGFOR,REGUNI,
     2                  REGLAB,REGCOL,ISTAT)           ! create the ident column
      CALL TBCINI(TIDO2,REGTYP,1,REGFOR,REGUNI,
     2                  REGLAB,REGCOL,ISTAT)           ! create the ident column
C
      DO 101 I = 1,NCPAR
         ICPAR(I) = I + 1
         REGTYP   = D_R4_FORMAT
         REGFOR   = FORMR4
         REGUNI   = UNIT(I)
         REGLAB   = LABEL(I)
         CALL TBCINI(TIDO1,REGTYP,1,REGFOR,REGUNI,
     2                     REGLAB,REGCOL,ISTAT)        ! create the data columns
         CALL TBCINI(TIDO2,REGTYP,1,REGFOR,REGUNI,
     2                     REGLAB,REGCOL,ISTAT)        ! create the data columns
  101 CONTINUE
C
C *** do the fitting
      DO 10 IR = 1,NROW
         IROW  = IR
         CALL TBRRDR(TIDI,IROW,2,ICBASE,COO,NUL,ISTAT)
         LX(1) = INT((COO(1)-BEGIN1(1))/STEP1(1)) + 1
         LY(1) = INT((COO(2)-BEGIN1(2))/STEP1(2)) + 1
         IROW  = IR + NROW
         CALL TBRRDR(TIDI,IROW,2,ICBASE,COO,NUL,ISTAT)
         LX(2) = INT((COO(1)-BEGIN2(1))/STEP2(1)) + 1
         LY(2) = INT((COO(2)-BEGIN2(2))/STEP2(2)) + 1
C
C *** fit the data; get first image
         KK   = 1
         IPNR = IPNTR1
         NPL0 = NPIX1(1)
         NL0  = NPIX1(2)
         ASSIGN 9001 TO FIT
         GO TO 9000
 9001    CONTINUE
C
C *** get the input of the second image
         KK   = 2 
         IPNR = IPNTR2
         NPL0 = NPIX2(1)
         NL0  = NPIX2(2)
         ASSIGN 9002 TO FIT
         GO TO 9000
 9002    CONTINUE
   10 CONTINUE
C
C *** write tables
      DO 40 IR = 1,NROW
         TAB(1)  = X1(IR,1)
         TAB(2)  = Y1(IR,1)
         TAB(3)  = H1(IR,1)
         TAB(4)  = FO1(IR,1)
         CALL TBEWRI(TIDO1,IR,ICIDN,IR,ISTAT)
         CALL TBRWRR(TIDO1,IR,12,ICPAR,TAB,ISTAT)
         TAB(1)  = X1(IR,2)
         TAB(2)  = Y1(IR,2)
         TAB(3)  = H1(IR,2)
         TAB(4)  = FO1(IR,2)
         CALL TBEWRI(TIDO2,IR,ICIDN,IR,ISTAT)
         CALL TBRWRR(TIDO2,IR,12,ICPAR,TAB,ISTAT)
   40 CONTINUE
C
      CALL TBTCLO(TIDO1,ISTAT)
      CALL TBTCLO(TIDO2,ISTAT)
      CALL STSEPI
C
C+++++
 9000 CONTINUE                                                     ! FIT DATA
C----
      DO 9110 J=1,9
         RX(J) = 0
         RY(J) = 0
 9110 CONTINUE
C
      NC = 0
      H1(IR,KK) = -10.E15
      DO 9120 I = LY(KK)-4,LY(KK)+4
         NC  = NC+1
         INX = LX(KK)-4
         CALL REALIN(NPL0,NL0,I,INX,9,MADRID(IPNR),VET)
         RME = 0
         DO 9121 J = 1,9
            H1(IR,KK)  = AMAX1(H1(IR,KK),VET(J))
            RME        = RME + VET(J)
            RX(J)      = RX(J) + VET(J)
 9121    CONTINUE
         RY(NC) = RME/9.
 9120 CONTINUE
C
      PAR(1) = 0.
      PAR(4) = 10.E15
      DO 9130 J = 1,9
         RX(J) = RX(J)/9.
         IF (RX(J).GT.PAR(1)) THEN
            PAR(1) = RX(J)
            PAR(2) = J
         END IF
         IF (RX(J).LT.PAR(4)) THEN
            PAR(4) = RY(J)
         ENDIF
 9130 CONTINUE
C
      PAR(1) = PAR(1)-PAR(4)
      KL     = PAR(2)
      VME    = (RX(KL-1)+RX(KL+1))/2.
      IF (VME.GE.PAR(1)) THEN
         VME = PAR(1)/2.
      ENDIF
      PAR(3) = SQRT(4*ALOG(2.)/ALOG(PAR(1)/VME))
      SQR    = 10.**15
      SL     = SQR
      IT     = 0
      NLIV   = 9
C
      DO 9140 J=1,4
         PARS(J) = PAR(J)
 9140 CONTINUE
C
 9201 CONTINUE
      IF (SL.LE.0.0001 .OR. IT.GT.20) GO TO 9200
      IT = IT+1
      CALL MONO4(IVX,RX,NLIV,PAR,.7)
      IF (PAR(1).LT..1 .OR. ABS(PAR(2)-PARS(2)).GT.4 .OR.
     2    PAR(3).LE.0. .OR. PAR(3).GT.20.) THEN
         DO 9141 J=1,4
            PAR(J) = PARS(J)
 9141    CONTINUE
         GO TO 9200
      END IF
C
      SQM = 0.
      DO 9150 I=1,NLIV
         FG  = PAR(1)*EXP(-4*ALOG(2.)*((I-PAR(2))/PAR(3))**2)
         SQM = SQM+(RX(I)-FG-PAR(4))**2
 9150 CONTINUE
      SQM = SQRT(SQM/NLIV)
      SL  = ABS(SQR-SQM)/SQR
      SQR = SQM
      GO TO 9201
C
 9200 CONTINUE
      X1(IR,KK)  = PAR(2)+LX(KK)-5
      FO1(IR,KK) = PAR(4)
      PAR(1)     = 0
      PAR(4)     = 10.E15
      DO 9160 J = 1,9
         IF (RY(J).GT.PAR(1)) THEN
            PAR(1) = RY(J)
            PAR(2) = J
         END IF
         IF (RY(J).LT.PAR(4)) THEN
            PAR(4) = RY(J)
         ENDIF
 9160 CONTINUE
C
      PAR(1) = PAR(1)-PAR(4)
      KL     = PAR(2)
      VME    = (RY(KL-1)+RY(KL+1))/2.
      IF (VME.GE.PAR(1)) VME = PAR(1)/2.
      PAR(3) = SQRT(4*ALOG(2.)/ALOG(PAR(1)/VME))
      SQR    = 10.**15
      SL     = SQR
      IT     = 0
      NLIV   = 9
C
      DO 9170 J=1,4
         PARS(J) = PAR(J)
 9170 CONTINUE
C
 9301 CONTINUE
      IF (SL.LE.0.0001.OR.IT.GT.20) GO TO 9300
      IT = IT+1
      CALL MONO4(IVX,RY,NLIV,PAR,.7)
      IF (PAR(1).LT..1 .OR. ABS(PAR(2)-PARS(2)).GT.4 .OR. 
     2    PAR(3).LE.0. .OR. PAR(3).GT.20.) THEN
         DO 9171 J=1,4
            PAR(J) = PARS(J)
 9171    CONTINUE
         GO TO 9300
      END IF
      SQM = 0.
      DO 9180 I=1,NLIV
         FG  = PAR(1)*EXP(-4*ALOG(2.)*((I-PAR(2))/PAR(3))**2)
         SQM = SQM+(RY(I)-FG-PAR(4))**2
 9180 CONTINUE
      SQM = SQRT(SQM/NLIV)
      SL  = ABS(SQR-SQM)/SQR
      SQR = SQM
      GO TO 9301
C
 9300 CONTINUE
      Y1(IR,KK)  = PAR(2)+LY(KK)-5
      FO1(IR,KK) = (FO1(IR,KK)+PAR(4))/2.
      H1(IR,KK)  = H1(IR,KK)-FO1(IR,KK)
C 
      GO TO FIT
      END
