C @(#)twodim.for 17.1.1.1 (ES0-DMD) 01/25/02 17:15:46 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 C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C C C----------------------------------------------------------------------- SUBROUTINE TWODIM(AS, MAS, JAPYS, IBUFS, LZ, & APSF, FPSF, IPSF, NCAT, PMTR, & PRCT, NCT, PMT, PRC, IARR, & RARR, SIGMA) C IMPLICIT NONE INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST' C REAL AS(1) INTEGER MAS((-MAXSUB):MAXSUB,(-MAXSUB):MAXSUB) INTEGER JAPYS(1) INTEGER IBUFS(4) INTEGER LZ REAL APSF(0:MAXSUB) REAL FPSF(1) INTEGER IPSF(1) INTEGER NCAT(NIPAR,MAXCNT) REAL PMTR(NRPAR,MAXCNT) REAL PRCT(0:MAXSUB,MAXCNT) INTEGER NCT(NIPAR) REAL PMT(NRPAR) REAL PRC(0:MAXSUB) INTEGER IARR(32) REAL RARR(64) REAL SIGMA C INTEGER NP INTEGER LPXL INTEGER LSBP INTEGER IC INTEGER ICM INTEGER IDF INTEGER ISC INTEGER IOFF, IOFF1, IOFF0, IOFFI, IOFF2, IOFFC INTEGER IS, ISOLD INTEGER ISTART, IACT INTEGER JCM INTEGER JDF INTEGER JSTART, JACT INTEGER I0, I0OLD INTEGER J0, J0OLD INTEGER JSC INTEGER JS, JSOLD INTEGER KCLN INTEGER KSAT INTEGER LPLIM INTEGER MSBP2, MSBP INTEGER MPXL, MPXL2 INTEGER MPRF INTEGER NNOR INTEGER NPAS INTEGER NPX, NPXOLD C REAL TMP(5) REAL BGRD REAL CORX, CORY REAL SCA, SCC, SCAA REAL SIGM, SHLM, SIGM1 REAL XDF REAL YDF DOUBLE PRECISION AR(25) , X(5) , S(5) LOGICAL CENTER , COMP , LAST , STAR C LOGICAL JCNTR NP = MAXCNT LPXL = IARR(20) LSBP = IARR(21) BGRD = PMTR(1,LZ) KCLN = NCAT(4,LZ) KSAT = NCAT(6,LZ) IF ( NCAT(10,LZ) .GT. 0 ) THEN COMP = .TRUE. ELSE COMP = .FALSE. ENDIF IF ( KSAT .EQ. -1 .AND. ABS(PMTR(3,LZ)) .LT. RARR(5) ) THEN STAR = .TRUE. ELSE STAR = .FALSE. ENDIF MSBP = 2 * LSBP + 1 MPXL = 2 * LPXL + 1 MSBP2 = MSBP * MSBP MPXL2 = MPXL * MPXL IOFF = MSBP2 * MPXL2 MPRF = NOSP + MSBP2 C C *** Determine scaling parameter SCA for object and SCC for component. C CALL SCLDET( LZ, NCAT , PMTR , PRCT , NCT , & PMT , PRC , APSF , NPAS , SCA , & SCC ) IF ( SCA .LE. 5.0*RARR(3) ) THEN STAR = .FALSE. ENDIF C C *** Save starting object coordinates. C ISTART = NCAT(1,LZ) JSTART = NCAT(2,LZ) C C *** Start loop. C I0 = 0 J0 = 0 I0OLD = 0 J0OLD = 0 IS = 0 ! addition for initialization JS = 0 ! addition for initialization ISOLD = IS JSOLD = JS NPX = 0 NPXOLD = -1 IC = 1 LAST = .FALSE. SIGM = BGRD C 10 CONTINUE IF ( IC .GT. 6 ) GOTO 11 C IF ( IC .EQ. 1 ) THEN SHLM = 0.6 ELSE SHLM = MIN( 0.6 , 1.0 / FLOAT(MSBP) ) ENDIF C C *** Identify object coordinates and subpixel. C CALL SUBPXL( NCAT(1,LZ) , PMTR(1,LZ) , LSBP , IACT , JACT , & IS , JS ) I0 = IACT - ISTART J0 = JACT - JSTART IF ( IC .EQ. 6 .OR. ( IC .GT. 1 .AND. I0 .EQ. I0OLD .AND. & J0 .EQ. J0OLD .AND. IS .EQ. ISOLD .AND. & JS .EQ. JSOLD .AND. NPX .EQ. NPXOLD ) ) THEN LAST = .TRUE. IC = 6 ENDIF NPXOLD = NPX C C *** Check component's distance. C IF ( COMP ) THEN XDF = PMT(10) - FLOAT(NCAT(1,LZ)) YDF = PMT(11) - FLOAT(NCAT(2,LZ)) IDF = NCT(1) - NCAT(1,LZ) JDF = NCT(2) - NCAT(2,LZ) IF ( ABS(IDF) .LE. LPXL .AND. ABS(JDF) .LE. LPXL ) THEN CENTER = .TRUE. C C *** Identify component's subpixel. C CALL SUBPXL( NCT , PMT , LSBP , ICM , JCM , & ISC , JSC ) ELSE CENTER = .FALSE. ENDIF ELSE CENTER = .FALSE. ENDIF C C *** IOFF0 - adress of psf for central pixel at given IS and JS. C IOFF0 = ( (LSBP+JS) * MSBP + LSBP + IS ) * MPXL2 + & LPXL * ( MPXL + 1 ) + 1 C C *** IOFFC - difference of adress of central C *** pixel at ISC and JSC minus IOFF0. C IF ( CENTER ) THEN IOFFC = ( (JSC-JS) * MSBP + ISC - IS ) * MPXL2 ENDIF C C *** IOFF1 and IOFF2 - corresponding adresses of x and y gradients. C IOFF1 = IOFF0 + IOFF IOFF2 = IOFF1 + IOFF C C *** Check if this object can serve as standard star. C LPLIM = INT( FLOAT(LPXL) * SQRT(2.0) ) IF ( STAR .AND. LAST .AND. KCLN .GT. LPLIM & .AND. NNOR .EQ. 4 ) THEN IF ( SIGM1 .GT. 0.0 ) THEN SCAA = SCA / SIGM1 ELSE SCAA = 0.0 ENDIF CALL IFSTAR( IBUFS , SCAA , LPXL , LSBP , IS , & JS , MPRF , IPSF , FPSF , IOFF0 , & IOFF , STAR , IOFFI ) IC = 7 ENDIF C C *** Calculate normal equations. C CALL NRMLEQ( AS , MAS , JAPYS , IBUFS , LPXL , & APSF , FPSF , IOFF0 , IOFF1 , IOFF2 , & IOFFI , BGRD , I0 , J0 , SCA , & COMP , SCC , CENTER , IDF , JDF , & XDF , YDF , IOFFC , IC , SIGM , & NPX , AR ) C C *** Solve normal equations if there is enough of good points. C IF ( COMP ) THEN NNOR = 5 ELSE NNOR = 4 ENDIF IF ( NPX .GT. NNOR .AND. ( .NOT. LAST ) ) THEN CALL LSQSOL ( NNOR , AR , NPX , X , S ) SIGM = SNGL(S(NNOR)) SCA = SCA + SNGL(X(1)) CORX = SNGL(X(2)) C C *** Save sigma from the first, non trimmed C *** determination of object intensity. C IF ( IC .EQ. 1 ) THEN SIGM1 = SNGL(S(1)) ENDIF IF ( ABS(CORX) .GT. SHLM ) THEN CORX = SHLM * CORX / ABS(CORX) ENDIF PMTR(10,LZ) = FLOAT(NCAT(1,LZ)) + & FLOAT(IS) / FLOAT(MSBP) - CORX PMTR(22,LZ) = SNGL(S(2)) NCAT(1,LZ) = NINT(PMTR(10,LZ)) CORY = SNGL(X(3)) IF ( ABS(CORY) .GT. SHLM ) THEN CORY = SHLM * CORY / ABS(CORY) ENDIF PMTR(11,LZ) = FLOAT(NCAT(2,LZ)) + & FLOAT(JS) / FLOAT(MSBP) - CORY PMTR(23,LZ) = SNGL(S(3)) NCAT(2,LZ) = NINT(PMTR(11,LZ)) IF ( COMP ) THEN SCC = SCC + SNGL(X(4)) ENDIF ENDIF C ISOLD = IS JSOLD = JS I0OLD = I0 J0OLD = J0 C IC = IC + 1 GOTO 10 11 CONTINUE C C *** Apply corrections. C IF ( NPXOLD .GT. NNOR ) THEN PMTR(12,LZ) = SCA IF ( SCA .GT. 0.0 ) THEN PMTR(24,LZ) = 1.086 * SNGL(S(1)) / SCA ELSE PMTR(24,LZ) = 1.0 ENDIF NCAT(9,LZ) = 1 IF ( NCAT(10,LZ) .GT. 0 ) THEN PMT(2) = SCC / ( 9.0 / ( 1.0 + 8.0 * APSF(1) ) ) PMT(12) = SCC ENDIF ELSE NCAT(9,LZ) = 0 ENDIF SIGMA = SNGL(S(NNOR)) C RETURN C END