C @(#)rfotcbase.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 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