C @(#)smoot.for	19.1 (ES0-DMD) 02/25/03 13:29: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
      SUBROUTINE SMOOT (RMAT,NY,NX)
C
C      11.X.88
C
      IMPLICIT  NONE
      INTEGER   NX
      INTEGER   NY
      REAL      RMAT(NY,NX)
      INTEGER   MX1
      INTEGER   MY1
      PARAMETER (MX1=55)
      PARAMETER (MY1=55)
C
      INTEGER   N35, N36
      INTEGER   I, J, L, K, INC, INR
      REAL      SP, S
      REAL      TMT(3,3)
      REAL      RMAT2(MY1,MX1)
C
      DATA      TMT/1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0/
C
      N35 = 3
      N36 = N35/2+1
      SP  = 10.
      DO 100 K = N36,NY-N36+1
         INC = K-N36
         DO 110 I = N36,NX-N36+1
            S   = 0.
            INR = I-N36
            DO 120 J = 1,N35
               DO 130 L = 1,N35
                  S = S+RMAT(INR+L,INC+J)*TMT(L,J)
  130          CONTINUE
  120       CONTINUE
            RMAT2(I,K) = S/SP
  110    CONTINUE
  100 CONTINUE
C
      DO 200 K = N36,NY-N36+1
         DO 210 I = N36,NX-N36+1
            RMAT(I,K) = RMAT2(I,K)
  210    CONTINUE
  200 CONTINUE

      RETURN
      END


      REAL FUNCTION GALE2(NPX,NPY,H,SIG,BET,XA,YA,X0,Y0)
C++++++++++
C     CALCOLA L' INTEGRALE BIDIMENSIONALE DI UNA
C     CURVA DI MOFFAT IN UN PIXEL CON IL METODO
C     DI GAUSS-LEGENDRE
C
C     H       = ALTEZZA       CURVA DI MOFFAT
C     SIG     = R               "    "    "
C     BET     = BETA            "    "    "
C     X0 - Y0 = COORD. CENTRO   "    "    "
C     XA - YA = COORD. INIZIALI DI INTEGRAZIONE
C     NPX     = NUMERO DI PUNTI INTERMEDI SULL' ASSE X
C               SU CUI CALCOLARE L' INTEGRALE NUMERICO (MAX=6)
C     NPY     = NUMERO DI PUNTI INTERMEDI SULL' ASSE Y
C               SU CUI CALCOLARE L' INTEGRALE NUMERICO (MAX=6)
C-------------
      IMPLICIT NONE
      INTEGER  NPX
      INTEGER  NPY
      REAL     H
      REAL     SIG
      REAL     BET
      REAL     XA
      REAL     YA
      REAL     X0
      REAL     Y0
C
      REAL     X, Y
      REAL     XB, YB
      REAL     XT, YT
      REAL     VT
      REAL     PSF
      INTEGER  I, J
      REAL     HE(6,6),HX(6,6)
C
      DATA     HX/ 0.0    ,0.0     ,0.0     ,0.0    ,0.0    ,0.0,
     *            -.577350,.577350 ,0.0     ,0.0    ,0.0    ,0.0,
     *            -.774597, 0.     ,.774597 ,0.0    ,0.0    ,0.0,
     *            -.861136,-.339981,.339981 ,.861136,0.0    ,0.0,
     *            -.906179,-.538469,0.      ,.538469,.906179,0.0,
     *            -.932469,-.661209,-.238619,.238619,.661209,.932469/
      DATA     HE/2.    ,0.0    ,0.0    ,0.0    ,0.0    ,0.0,
     *            1.    ,1.     ,0.0    ,0.0    ,0.0    ,0.0,
     *           .555555,.888889,.555555,0.0    ,0.0    ,0.0,
     *           .347855,.652145,.652145,.347855,0.0    ,0.0,
     *           .236927,.478629,.568889,.478629,.236927,0.0,
     *           .171324,.360762,.467914,.467914,.360762,.171324/
C
      PSF(H,X,Y,X0,Y0,SIG,BET) = H*(1+((X-X0)**2+(Y-Y0)**2)
     *                            /SIG**2)**(-BET)
      XB = XA+1
      YB = YA+1
      VT = 0.
C
      DO 10  I = 1,MIN0(6,NPX)
         XT = .5*(HX(I,NPX)+XB+XA)
         DO 20  J = 1,MIN0(6,NPY)
            YT = .5*(HX(J,NPY)+YB+YA)
            VT = VT+HE(I,NPX)*HE(J,NPY)*PSF(H,XT,YT,X0,Y0,SIG,BET)
  20     CONTINUE
  10  CONTINUE
      GALE2 = VT/4.

      RETURN
      END


      REAL FUNCTION GALEG(NPX,NPY,H,SIG,XA,YA,X0,Y0)
C
C     CALCOLA L' INTEGRALE BIDIMENSIONALE DI UNA
C     CURVA DI MOFFAT IN UN PIXEL CON IL METODO
C     DI GAUSS-LEGENDRE
C
C     H       = ALTEZZA       CURVA DI MOFFAT
C     SIG     = R               "    "    "
C     X0 - Y0 = COORD. CENTRO   "    "    "
C     XA - YA = COORD. INIZIALI DI INTEGRAZIONE
C     NPX     = NUMERO DI PUNTI INTERMEDI SULL' ASSE X
C               SU CUI CALCOLARE L' INTEGRALE NUMERICO (MAX=6)
C     NPY     = NUMERO DI PUNTI INTERMEDI SULL' ASSE Y
C               SU CUI CALCOLARE L' INTEGRALE NUMERICO (MAX=6)
C
      IMPLICIT NONE
      INTEGER  NPX
      INTEGER  NPY
      REAL     H
      REAL     SIG
      REAL     XA
      REAL     YA
      REAL     X0
      REAL     Y0
C
      REAL     X, Y
      REAL     XB, YB
      REAL     XT, YT
      REAL     VT
      REAL     PSF
      INTEGER  I, J
      REAL     HE(6,6),HX(6,6)
C
      DATA     HX/ 0.0    ,0.0     ,0.0     ,0.0    ,0.0    ,0.0,
     *            -.577350,.577350 ,0.0     ,0.0    ,0.0    ,0.0,
     *            -.774597, 0.     ,.774597 ,0.0    ,0.0    ,0.0,
     *            -.861136,-.339981,.339981 ,.861136,0.0    ,0.0,
     *            -.906179,-.538469,0.      ,.538469,.906179,0.0,
     *            -.932469,-.661209,-.238619,.238619,.661209,.932469/
      DATA     HE/2.    ,0.0    ,0.0    ,0.0    ,0.0    ,0.0,
     *            1.    ,1.     ,0.0    ,0.0    ,0.0    ,0.0,
     *           .555555,.888889,.555555,0.0    ,0.0    ,0.0,
     *           .347855,.652145,.652145,.347855,0.0    ,0.0,
     *           .236927,.478629,.568889,.478629,.236927,0.0,
     *           .171324,.360762,.467914,.467914,.360762,.171324/
C
      PSF(H,X,Y,X0,Y0,SIG) = H*EXP(-4*ALOG(2.)*((X-X0)**2+(Y-Y0)**2)
     *                        /SIG**2)
      XB = XA+1
      YB = YA+1
      VT = 0.
C
      DO 10  I=1,MIN0(6,NPX)
         XT = .5*(HX(I,NPX)+XB+XA)
         DO 20 J=1,MIN0(6,NPY)
            YT = .5*(HX(J,NPY)+YB+YA)
            VT = VT+HE(I,NPX)*HE(J,NPY)*PSF(H,XT,YT,X0,Y0,SIG)
   20    CONTINUE
   10 CONTINUE
C
      GALEG=VT/4.
      RETURN
      END
