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