C @(#)nrmleq.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15:43
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-----------------------------------------------------------------------
      SUBROUTINE 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
      IMPLICIT           NONE
      INCLUDE            'MID_REL_INCL:INVENT.INC/NOLIST'
C
      REAL               AS(1)             ! IN: Data array
      INTEGER            MAS((-MAXSUB):MAXSUB,(-MAXSUB):MAXSUB) ! Mask array
      INTEGER            JAPYS(1)          ! IN: Pointers to data lines
      INTEGER            IBUFS(4)          ! IN: Limits of image buffer
      INTEGER            LPXL              ! IN: Pixel's extend
      REAL               APSF(0:MAXSUB)    ! IN: One dimensional p.s.f.
      REAL               FPSF(1)           ! MOD: Two dimensional p.s.f.
      INTEGER            IOFF0             ! IN: P.s.f. offset
      INTEGER            IOFF1             ! IN: X gradient offset
      INTEGER            IOFF2             ! IN: Y gradient offset
      INTEGER            IOFFI             ! IN: Corrections offset
      REAL               BGRD              ! IN: Sky background
      INTEGER            I0                ! IN: Object x shift
      INTEGER            J0                ! IN: Object y shift
      REAL               SCA               ! IN: Object intensity
      LOGICAL            COMP              ! IN: TRUE if there is a component
      REAL               SCC               ! IN: Component intensity
      LOGICAL            CENTER            ! IN: TRUE if close component
      INTEGER            IDF               ! IN: X pixel offset of component
      INTEGER            JDF               ! IN: Y pixel offset of component
      REAL               XDF               ! IN: X real offset of component
      REAL               YDF               ! IN: Y real offset of component
      INTEGER            IOFFC             ! IN: Adress offset of component
      INTEGER            IC                ! IN: Iteration number
      REAL               SIGM              ! IN: Sigma of pixel data
      INTEGER            NPX               ! IN: Number of good pixels
      DOUBLE PRECISION   AR(25)            ! OUT: Normal equations coefficients
C
      INTEGER   I , I1 , I2 , IDIF , IOF
      INTEGER   J , J1 , J2 , JARG , JOF , JOFC
      INTEGER   K , L , MPXL , NNOR , NR
      REAL      ADIF , DIFF , TMP(5)
      LOGICAL   JCNTR
C
      MPXL = 2 * LPXL + 1
C
C *** Prepare for calculations.
C
      DO 50 J = 1 , 25
          AR(J) = 0.0
   50 CONTINUE
      NPX = 0
      I1 = MAX( (-LPXL) , IBUFS(1)-I0 )
      I2 = MIN( LPXL , IBUFS(3)-I0 )
      J1 = MAX( (-LPXL) , IBUFS(2)-J0 )
      J2 = MIN( LPXL , IBUFS(4)-J0 )
C
C *** Form normal equations for least squares solution.
C
      DO 10 J = J1 , J2
          JARG = JAPYS(J+J0-IBUFS(2)+1) + I0
          JOF = MPXL * J
          IF ( CENTER .AND. ABS(J+JDF) .LE. LPXL ) THEN
              JCNTR = .TRUE.
              JOFC = MPXL * (J+JDF)
          ELSE
              JCNTR = .FALSE.
          ENDIF
          DO 20 I = I1 , I2
              IF ( MAS(I,J) .GE. 1 ) THEN
                  IOF = JOF + I
                  TMP(1) = FPSF(IOFF0+IOF)
                  TMP(2) = FPSF(IOFF1+IOF) * SCA
                  TMP(3) = FPSF(IOFF2+IOF) * SCA
                  IF ( COMP ) THEN
                      IF ( JCNTR .AND. ABS(I+IDF) .LE. LPXL ) THEN
                          TMP(4) = FPSF(IOFF0+IOFFC+JOFC+I+IDF)
                      ELSE
                          DIFF = SQRT( (XDF-FLOAT(I))**2.0 +
     &                                 (YDF-FLOAT(J))**2.0 )
                          IDIF = INT(DIFF)
                          IF ( IDIF .LT. MAXSUB ) THEN
                              ADIF = DIFF - IDIF
                              TMP(4) = APSF(IDIF) * (1.0-ADIF) +
     &                                 APSF(IDIF+1) * ADIF
                          ELSE
                              TMP(4) = 0.0
                          ENDIF
                      ENDIF
                      TMP(5) = AS(JARG+I) - BGRD - TMP(1) * SCA -
     &                                             TMP(4) * SCC
                      NNOR = 5
                  ELSE
                      TMP(4) = AS(JARG+I) - BGRD - TMP(1) * SCA
                      NNOR = 4
                  ENDIF
                  IF ( ABS(TMP(NNOR)) .LT. 2.0 * SIGM ) THEN
                      NPX = NPX + 1
                      DO 70 L = 1 , NNOR
                          DO 80 K = L , NNOR
                              NR = (L-1)*NNOR + K
                              AR(NR) = AR(NR) + DBLE(TMP(K)*TMP(L))
   80                     CONTINUE
   70                 CONTINUE
                  ENDIF
                  IF ( IC .EQ. 7 ) THEN
                      FPSF(IOFFI+IOF) = TMP(NNOR) / SCA
                  ENDIF
              ENDIF
   20     CONTINUE
   10 CONTINUE
C
      RETURN
C
      END
C