C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LFN_PLFIT       **
C***************************
C
C
      SUBROUTINE LFN_PLFIT ( IPASS)
C
      INCLUDE 'implicit_none.finc'
C
C PURPOSE:       Fits plane A*X+B*Y+C through background
C
C **** Service Routine Only ****
C
C     Minimize sum[(I - B1*X - B2*Y - B3)**2]
C     where I is intensity at X,Y and
C     B1,B2,B3 are constants to be determined.
C
C     The normal matrix PLFMAT is:               RHO is:
C
C     sum(X*X)  sum(X*Y)  sum(X)       sum(I*X)
C     sum(X*Y)  sum(Y*Y)  sum(Y)       sum(I*Y)
C     sum(X)    sum(Y)    sum(1) = n   sum(I)
C
C     and we have PLFMAT * B = RHO
C                 ======   -   ---
C
C                    -1
C     Thus  B = PLFMAT  * RHO
C           -   ======    ---
C
C Author: Q. Hao
C Based on code from MADNES
C     M. Z. Papiz             Daresbury Laboratory
C     J. W. Pflugrath         Cold Spring Harbor Laboratory
C
C====== ARGUMENTS: 
C
      INTEGER IPASS
C     INTEGER IPASS (R)    =1, fit background plane
C                          =2, correct background plane
C====== Parameter statements
C
      INCLUDE 'saint_limits.finc'
C
      INTEGER    QX, QY, QXX, QYY, QXY, QN
      PARAMETER (QXX=1, QXY=2, QX=3, QYY=4, QY=5, QN=6)
C
C====== GLOBALS:
C 
      INCLUDE 'lfn_saint.finc' 
C
C====== LOCALS:
C
      INTEGER I
      INTEGER J
      INTEGER IMVEC(6)
      REAL PLFMAT(3,3)
      REAL PLFINV(3,3)
      REAL RHO(3)
      SAVE RHO
C
C     INTEGER I         !Temp/loop variable
C     INTEGER J         !Temp/loop variable
C     INTEGER IMVEC(6)	!Array to hold elements of normal matrix
C     REAL PLFMAT(3,3)	!Normal matrix
C     REAL PLFINV(3,3)	!Inverted normal matrix
C     REAL RHO(3)	!Array to hold sum(I*X), sum(I*Y), sum(I)
C
C------------------------------------------------------------------------------ 
C
C====== Clear least-squares sums (ONLY IN PASS 1)
C
      IF (IPASS .EQ. 1) THEN
         DO 100 I = 1, 6
100         NMVEC(I) = 0
         CONTINUE
         DO 110 I = 1 , 3
110         RHO(I) = 0.0
         CONTINUE
         IMVEC(QN) = 1
      ENDIF
C
      IF (IPASS .EQ. 1)  THEN
C
C====== Collect background sums for least squares refinement : PASS 1
C
            DO 300 I = 1, NUM_IN_SPOT
                  IF (SHAPE(I).LT.0) GO TO 300
                  IF (SHAPE(I) .LE. 1) THEN
C
C====== This is a background pixel
C
                     IMVEC(QXX) = IXP(I)**2
                     IMVEC(QXY) = IXP(I)*IYP(I)
                     IMVEC(QYY) = IYP(I)**2
                     IMVEC(QY) = IYP(I)
                     IMVEC(QX) = IXP(I)
                     RHO(1) = RHO(1)  +  PIXV(I) * IXP(I)
                     RHO(2) = RHO(2)  +  PIXV(I) * IYP(I)
                     RHO(3) = RHO(3)  +  PIXV(I)
                     DO 200 J = 1, 6
                        NMVEC(J) = NMVEC(J) + IMVEC(J)
200                  CONTINUE
                  ENDIF
300            CONTINUE
      ELSE IF ( (IPASS .EQ. 2) .AND. (NREJ .GT. 0) ) THEN
C
C====== Correct background sums for least squares refinement : PASS 2
C
         DO 350 J = 1, NREJ
            I  = IREJ(J)
            NMVEC(QYY) = NMVEC(QYY) - IYP(I)**2
            NMVEC(QY)  = NMVEC(QY)  - IYP(I)
            NMVEC(QXY) = NMVEC(QXY) - IYP(I)*IXP(I)
            NMVEC(QX)  = NMVEC(QX)  - IXP(I)
            NMVEC(QXX) = NMVEC(QXX) - IXP(I)**2
            RHO(1) = RHO(1)  -  PIXV(I) * IXP(I)
            RHO(2) = RHO(2)  -  PIXV(I) * IYP(I)
            RHO(3) = RHO(3)  -  PIXV(I)
350      CONTINUE
         NMVEC(QN) = NMVEC(QN) - NREJ
      ENDIF
      PLFMAT(1,1) = NMVEC(QXX)
      PLFMAT(1,2) = NMVEC(QXY)
      PLFMAT(1,3) = NMVEC(QX)
      PLFMAT(2,1) = NMVEC(QXY)
      PLFMAT(2,2) = NMVEC(QYY)
      PLFMAT(2,3) = NMVEC(QY)
      PLFMAT(3,1) = NMVEC(QX)
      PLFMAT(3,2) = NMVEC(QY)
      PLFMAT(3,3) = NMVEC(QN)
C
C====== Invert normal matrix and solve equations
C
      CALL MTX_INV(PLFMAT, PLFINV, 3)
C
C====== Calculate plane fit constants A,B,C in BACK(1) ->
C
      DO 400 I = 1, 3
         BACK(I) = 0.0
         DO 400 J = 1, 3
            BACK(I) = BACK(I) + PLFINV(I, J) * RHO(J)
400      CONTINUE
      CONTINUE
      RETURN
      END
