C @(#)defpsf.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15: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 DEFPSF(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   MPXL, MPXL2
      INTEGER   MSBP, MSBP2
      INTEGER   I
      INTEGER   ISTEP
      INTEGER   IDI, IDJ
      INTEGER   I1, I2
      INTEGER   LMAX
      INTEGER   LADR, JADR, KADR
      INTEGER   L
      INTEGER   K, KK
      INTEGER   J
      INTEGER   IADR
C
      REAL      SCL
      REAL      DR, DL
      REAL      RD, RDD
      REAL      ACN, ASN
      REAL      APSF(0:MAXSUB)
      REAL      COEF(2,MAXSUB)
      REAL      ZERO
      REAL      SCL2
C
C *** Recall one-dimensional point spread function.
C
      CALL INAPSF( RARR , MAXSUB , APSF )
C
C *** Calculate dimensions of arrays.
C
      MPXL = 2 * LPXL + 1
      MSBP = 2 * LSBP + 1
      MPXL2 = MPXL * MPXL
      MSBP2 = MSBP * MSBP
      ISTEP = MSBP2 * MPXL2
      IDI = ISTEP
      IDJ = 2 * ISTEP
C
C *** Initialize array IPSF.
C
      IPSF(0) = NOSP
      I1 = IPSF(0) + 1
      I2 = IPSF(0) + MSBP2 
      DO 60 I = I1 , I2
          IPSF(I) = 0
   60 CONTINUE
C
C *** Approximate to parabola inside of central pixel.
C
      SCL = ( 1.0 - APSF(1) ) / 1.496
      ZERO = 1.0 + 0.168 * SCL
      SCL2 = (-2.0) * SCL
C
C *** Calculate coefficients of parabolic interpolation.
C
      COEF(1,2) = 0.415*(APSF(3)-APSF(2)) - 0.824*(APSF(1)-APSF(2))
      COEF(2,2) = 0.585*(APSF(3)-APSF(2)) + 0.824*(APSF(1)-APSF(2))
      LMAX = MIN( 2*LPXL , MAXSUB-1 )
      DO 15 L = 3 , LMAX
          COEF(1,L) = ( APSF(L+1) - APSF(L-1) ) / 2.0
          COEF(2,L) = ( APSF(L+1) + APSF(L-1) - 2*APSF(L) ) / 2.0
   15 CONTINUE
C
C *** Calculate initial two-dimensional point spread function.
C
      DL = 1.0 / FLOAT( 2 * LSBP + 1 )
      DO 20 L = -LSBP , LSBP
          LADR = 1 + (L+LSBP) * MSBP * MPXL2
          DO 30 K = -LSBP , LSBP
              KADR = LADR + (K+LSBP) * MPXL2
              DO 40 J = -LPXL , LPXL
                  JADR = KADR + (J+LPXL) * MPXL
                  DO 50 I = -LPXL , LPXL
                      IADR = JADR + (I+LPXL)
                      RDD = (J-L*DL)**2.0 + (I-K*DL)**2.0
                      IF ( RDD .EQ. 0.0 ) THEN
                          RDD = 0.168 * DL * DL
                      ENDIF
                      RD = SQRT( RDD )
                      IF ( RD .GT. 0.0 ) THEN
                          ACN = (I-K*DL) / RD
                          ASN = (J-L*DL) / RD
                      ELSE
                          ACN = 0.0
                          ASN = 0.0
                      ENDIF
                      IF ( RD .LT. 1.29 ) THEN
                          FPSF(IADR) = ZERO - SCL * RDD
                          FPSF(IADR+IDI) = SCL2 * RD * ACN
                          FPSF(IADR+IDJ) = SCL2 * RD * ASN
                      ELSE IF ( RD .LT. 2.6 ) THEN
                          DR = RD - 2.0
                          FPSF(IADR) = APSF(2) + ( COEF(1,2)
     &                                 + COEF(2,2) * DR ) * DR
                          FPSF(IADR+IDI) = ( COEF(1,2) +
     &                                     2.0 * COEF(2,2) * DR ) * ACN
                          FPSF(IADR+IDJ) = ( COEF(1,2) +
     &                                     2.0 * COEF(2,2) * DR ) * ASN
                      ELSE
                          KK = NINT(RD)
                          DR = RD - KK
                          FPSF(IADR) = APSF(KK) + ( COEF(1,KK)
     &                                 + COEF(2,KK) * DR ) * DR
                          FPSF(IADR+IDI) = ( COEF(1,KK) +
     &                                     2.0 * COEF(2,KK) * DR ) * ACN
                          FPSF(IADR+IDJ) = ( COEF(1,KK) +
     &                                     2.0 * COEF(2,KK) * DR ) * ASN
                      ENDIF
   50             CONTINUE
   40         CONTINUE
   30     CONTINUE
   20 CONTINUE
C
      RETURN
C
      END
C