C @(#)smtpsf.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15:45
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 SMTPSF(LPXL, LSBP, LL, CPSF, DPSF,
     &                  IPSF, INZ, NN, SIGMA, I,
     &                  J, K, L, AMAG, GRI,
     &                  GRJ)
C
      IMPLICIT  NONE
C
      INTEGER   LPXL
      INTEGER   LSBP
      INTEGER   LL
      REAL      CPSF((-LL):LL,(-LL):LL)
      REAL      DPSF((-LL):LL,(-LL):LL)
      INTEGER   IPSF((-LSBP):LSBP,(-LSBP):LSBP)
      INTEGER   INZ
      INTEGER   NN
      REAL      SIGMA
      INTEGER   I
      INTEGER   J
      INTEGER   K
      INTEGER   L
      REAL      AMAG
      REAL      GRI
      REAL      GRJ
C
      INTEGER   MSBP, MSBP2
      INTEGER   MPXL, MPXL2
      INTEGER   II, II1, II2, III
      INTEGER   IC, ISTEP, ITMP
      INTEGER   NJ, NI
      INTEGER   LLL
      INTEGER   KKK	
      INTEGER   JJ, JJ1, JJ2, JJJ
      REAL      DI, DJ
C
      DOUBLE PRECISION   COEF(6,6) , S(6) , TEMP , TMP(6) , X(6)
C
      IF ( NN .LT. 7 ) THEN
          AMAG = 0.0
          GRI = 0.0
          GRJ = 0.0
          RETURN
      ENDIF
      MSBP = 2 * LSBP + 1
      MPXL = 2 * LPXL + 1
      MSBP2 = MSBP * MSBP
      MPXL2 = MPXL * MPXL
      ISTEP = MAX( 2 , NINT( ( SIGMA / 0.04 ) *
     &                       SQRT(FLOAT(MSBP2*MPXL2)/FLOAT(NN)) ) )
      II = I * MSBP - K
      JJ = J * MSBP - L
      IF ( DPSF(II,JJ) .LT. 0.01 ) THEN
          AMAG = CPSF(II,JJ)
          RETURN
      ENDIF
   71 CONTINUE
      DO 10 NJ = 1 , 6
          DO 20 NI = 1 , 6
              COEF(NI,NJ) = 0.0
   20     CONTINUE
   10 CONTINUE
   70 CONTINUE
          II1 = MAX( -LL , II-ISTEP )
          II2 = MIN( LL , II+ISTEP )
          IF ( II1 .EQ. -LL .AND. II2-II1 .LE. 2 ) THEN
              II2 = II1 + 3
          ENDIF
          IF ( II2 .EQ. LL .AND. II2-II1 .LE. 2 ) THEN
              II1 = II2 - 3
          ENDIF
          JJ1 = MAX( -LL , JJ-ISTEP )
          JJ2 = MIN( LL , JJ+ISTEP )
          IF ( JJ1 .EQ. -LL .AND. JJ2-JJ1 .LE. 2 ) THEN
              JJ2 = JJ1 + 3
          ENDIF
          IF ( JJ2 .EQ. LL .AND. JJ2-JJ1 .LE. 2 ) THEN
              JJ1 = JJ2 - 3
          ENDIF
          IC = 0
          DO 30 JJJ = JJ1 , JJ2
              IF ( JJJ .LT. 0 ) THEN
                  LLL = MOD( JJJ-LSBP , MSBP ) + LSBP
              ELSE IF ( JJJ .EQ. 0 ) THEN
                  LLL = 0
              ELSE
                  LLL = MOD( JJJ+LSBP , MSBP ) - LSBP
              ENDIF
              DO 40 III = II1 , II2
                  IF ( III .LT. 0 ) THEN
                      KKK = MOD( III-LSBP , MSBP ) + LSBP
                  ELSE IF ( III .EQ. 0 ) THEN
                      KKK = 0
                  ELSE
                      KKK = MOD( III+LSBP , MSBP ) - LSBP
                  ENDIF
                  ITMP = IPSF(KKK,LLL)
                  IF ( ITMP .GT. 0 ) THEN
                      TEMP = SQRT( DBLE( ITMP ) )
                      IC = IC + 1
                      DI = DBLE(III-II)
                      DJ = DBLE(JJJ-JJ)
                      TMP(1) = 1.0 * TEMP
                      TMP(2) = DI * TEMP
                      TMP(3) = DJ * TEMP
                      TMP(4) = DI * DI * TEMP
                      TMP(5) = DJ * DJ * TEMP
                      TMP(6) = DBLE( CPSF(III,JJJ) ) * TEMP
                      DO 50 NJ = 1 , 6
                          DO 60 NI = NJ , 6
                              COEF(NI,NJ) = COEF(NI,NJ) +
     &                                      TMP(NI) * TMP(NJ)
   60                     CONTINUE
   50                 CONTINUE
                  ENDIF
   40         CONTINUE
   30     CONTINUE
          IF ( IC .LT. 7 .AND. ISTEP .LE. LL/4 ) THEN
              ISTEP = ISTEP + 1
              GOTO 70
          ENDIF
      CONTINUE
      IF ( IC .GE. 7 ) THEN
          CALL LSQSOL( 6 , COEF , IC , X , S )
          IF ( SNGL(S(1)) .GT. 0.05 .AND. ISTEP .LE. LL/4 ) THEN
              ISTEP = ISTEP + 1
              GOTO 71
          ELSE
              AMAG = SNGL(X(1))
              GRI = SNGL(X(2)) * FLOAT(MSBP)
              GRJ = SNGL(X(3)) * FLOAT(MSBP)
          ENDIF
      ELSE
          AMAG = CPSF(II,JJ)
          GRI = 0.0
          GRJ = 0.0
      ENDIF
C
      RETURN
C
      END
C