C @(#)modpsf.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.INPUT/OUTPUT
C  Input parameters
C  RARR      real*4 array         real parameters
C  MAXSUB    integer*4            maximal subarray half-edge
C  IPSF      integer*4            pointers to FPSF
C  LPXL      integer*4            extend of 2-d point spread function
C  LSBP      integer*4            number of subpixels
C  Output parameters
C  FPSF      real*4 array         2-d point spread function
C-----------------------------------------------------------------------
      SUBROUTINE MODPSF(RARR, FPSF, NPRF, IPSF, LPXL, LSBP)
C
      IMPLICIT    NONE
      INCLUDE     'MID_REL_INCL:INVENT.INC/NOLIST'
C
      REAL      RARR(64)
      REAL      FPSF(1)
      INTEGER   NPRF
      INTEGER   IPSF(0:NPRF)
      INTEGER   LPXL
      INTEGER   LSBP
C
      INTEGER   I, II
      INTEGER   IT1, IT2, IT3, IT4
      INTEGER   IDM, IDN, IDK1, IDK2, INZ
      INTEGER   IOFN, IOF4, IOF5, IOF42, IOF6
      INTEGER   K
      INTEGER   J
      INTEGER   L, LL
      INTEGER   MPXL, MPXL2
      INTEGER   MSBP, MSBP2
      INTEGER   NN, NI
      INTEGER   NOSP1
      REAL      AMEAN	
      REAL      GRI, GRJ
      REAL      SIGMA
      REAL      TEMP
      REAL      X
C
C ***
      MPXL  = 2 * LPXL + 1
      MSBP  = 2 * LSBP + 1
      MPXL2 = MPXL * MPXL
      MSBP2 = MSBP * MSBP
      IDM   = MPXL2 * MSBP2
      IOF4  = 3 * IDM + 1
      IOF5  = IOF4 + IDM
      IOF6  = IOF5 + IDM
      IDN   = 2 * IDM
C
C *** Calculate corrections to point spread function.
      LL = ( MPXL * MSBP ) / 2
      CALL CORPSF( LPXL , LSBP , LL , FPSF(IOF6) , IPSF(NOSP+1) ,
     &             RARR(39) , FPSF(IOF4) , FPSF(IOF5) )
C
C *** Check data accuracy.
      CALL MEAN( FPSF(IOF4) , IDM , AMEAN , SIGMA )
C
C *** Correct for zeroes in data.
      IDK1 = NOSP + 1
      IDK2 = NOSP + MSBP2
      INZ = 0
      DO 10 II = IDK1 , IDK2
          IF ( IPSF(II) .GT. 0 ) THEN
              INZ = INZ + 1
          ENDIF
   10 CONTINUE
      IF ( INZ .GT. 0 ) THEN
          SIGMA = SIGMA * SQRT( FLOAT(IDM-1) / FLOAT(INZ*MPXL2-1) )
      ELSE
          RETURN
      ENDIF
C
C *** Count non zero values of corrections.
      NN = 0
      IOF42 = IOF5 - 1
      DO 80 NI = IOF4 , IOF42
          IF ( FPSF(NI) .NE. 0.0 ) THEN
              NN = NN + 1
          ENDIF
   80 CONTINUE
C
      IF ( NN .LT. 7 ) RETURN
C
C     CALL MEAN( FPSF(IOF4) , IDM , AMEAN , SIGMA )
C
C     SIGMA = SIGMA * SQRT( FLOAT(IDM-1) / FLOAT(NN-1) )
C
C *** Smooth corrections to point spread function
C *** and calculate corrections to gradients.
C
      IF ( INZ .LT. 1 .OR. INZ*MPXL2 .LT. 9 ) RETURN
      DO 20 L = -LSBP , LSBP
          IT1 = (L+LSBP) * MSBP + LSBP
          DO 30 K = -LSBP , LSBP
              IT2 = ( IT1 + K ) * MPXL2
              DO 40 J = -LPXL , LPXL
                  IT3 = IT2 + (J+LPXL) * MPXL + LPXL + 1
                  DO 50 I = -LPXL , LPXL
                      IT4 = IT3 + I
                      CALL SMTPSF( LPXL, LSBP, LL, FPSF(IOF4),
     &                             FPSF(IOF5), IPSF(NOSP+1), INZ,
     &                             NN, SIGMA, I, J, K, L, X, GRI, GRJ)
                      FPSF(IT4) = FPSF(IT4) + X
                      TEMP = FPSF(IT4+IDM) + GRI
                      IF ( TEMP*FLOAT(I) .GT. 0.0 ) THEN
                          TEMP = 0.0
                      ENDIF
                      FPSF(IT4+IDM) = TEMP
                      TEMP = FPSF(IT4+IDN) + GRJ
                      IF ( TEMP*FLOAT(J) .GT. 0.0 ) THEN
                          TEMP = 0.0
                      ENDIF
                      FPSF(IT4+IDN) = TEMP
   50             CONTINUE
   40         CONTINUE
   30     CONTINUE
   20 CONTINUE
C
C *** Adjust one- and two-dimensional point spread functions.
      CALL ADJPSF( RARR , LPXL , LSBP , FPSF )
C
      IOFN = IOF6 + IDM * (NOSP+2) - 1
      DO 60 II = IOF6 , IOFN
          FPSF(II) = 0.0
   60 CONTINUE
      NOSP1 = NOSP+1
      DO 70 II = NOSP1 , NPRF
          IPSF(II) = 0
   70 CONTINUE
C
      RETURN
      END