C @(#)elmoi.for	19.1 (ESO-DMD) 02/25/03 13:29:42
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 ELMIF(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano fisso
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VMO, VG, VTP
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK, INK1
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      NIN = NC*3
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VPI     = VP(THREE)
      VTN(1)  = 1.
      DO 30  K=1,NP
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N  = L*3-2
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = .5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = .5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N)   = EP*HH+VTN(N)
                  VTN(N+1) = CO*DIX*HH+VTN(N+1)
                  VTN(N+2) = CO*DIY*HH+VTN(N+2)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N)   = VTN(N)/4
            VTN(N+1) = VTN(N+1)/4
            VTN(N+2) = VTN(N+2)/4
   32    CONTINUE
   30 CONTINUE
C
      VFU = VFU/4+VPI
      DO 41 I = 1,NIN
         VTP   = VTN(I)*WEI(K)
         VC(I) = VC(I)+(VZ(K)-VFU)*VTP
         DO 42 J = 1,I
            RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42    CONTINUE
   41 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         DO 71 I=1,NC
            DO 72 J=4,6
               INK     = J+4*(I-1)
               INK1    = J+3*(I-2)
               VMO     = VC(INK1)*VPES(J)
               VP(INK) =VMO+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM=0.
            DO 81 I = 1,NP
               VG=0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = .5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                          EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                          (1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
C
               VG  = VG/4.+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIFV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano orizzontale
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VMO, VG, VTP
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK, INK1
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      NIN = NC*4
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VPI     = VP(THREE)
      VTN(1)  = 1.
      DO 30  K=1,NP
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N  = L*4-3
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = 0.5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N)   = EP*HH+VTN(N)
                  VTN(N+1) = CO*DIX*HH+VTN(N+1)
                  VTN(N+2) = CO*DIY*HH+VTN(N+2)
                  VTN(N+3) = CO*EPES*HH/VP(N1+3)+VTN(N+3)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N)   = VTN(N)/4
            VTN(N+1) = VTN(N+1)/4
            VTN(N+2) = VTN(N+2)/4
            VTN(N+3) = VTN(N+3)/4
   32    CONTINUE
C
         VFU = VFU/4+VPI
         DO 41 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
            DO 42 J = 1,I
               RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42       CONTINUE
   41    CONTINUE
   30 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)*VPES(3)*VP(THREE)
         DO 71 I=1,NC
            DO 72 J=4,7
               INK     = J+4*(I-1)
               INK1    = INK-3
               VMO     = VC(INK1)*VPES(J)
               VP(INK) = VMO+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM=0.
            DO 81 I = 1,NP
               VG=0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = .5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                          EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                          (1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
C
               VG  = VG/4.+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIR(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano fisso
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VMO, VG, VTP, VVV
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK, INK1
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      NIN = NC*3+1
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VPI     = VP(THREE)
      VTN(1)  = 1.
      DO 30  K=1,NP
         VFU = 0.0
         VVV = 0.0
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N  = L*3-1
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = .5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = .5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N)   = EP*HH+VTN(N)
                  VTN(N+1) = CO*DIX*HH+VTN(N+1)
                  VTN(N+2) = CO*DIY*HH+VTN(N+2)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N)   = VTN(N)/4
            VTN(N+1) = VTN(N+1)/4
            VTN(N+2) = VTN(N+2)/4
   32    CONTINUE
C
         VFU = VFU/4+VPI
         DO 41 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
            DO 42 J = 1,I
               RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42       CONTINUE
   41    CONTINUE
   30 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD=NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
         DO 71 I=1,NC
            DO 72 J=4,6
               INK     = J+4*(I-1)
               INK1    = J+3*(I-2)-2
               VMO     = VC(INK1)*VPES(J)
               VP(INK) = VMO+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 81 I = 1,NP
               VG = 0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = .5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*
     *                          VP(N)*EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*
     2                          VP(N)*(1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
               VG  = VG/4.+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIRV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano orizzontale
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VMO, VG, VTP
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK, INK1
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      NIN = NC*4+1
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VPI     = VP(THREE)
      VTN(1)  = 1.
      DO 30  K=1,NP
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N  = L*4-2
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = 0.5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N)   = EP*HH+VTN(N)
                  VTN(N+1) = CO*DIX*HH+VTN(N+1)
                  VTN(N+2) = CO*DIY*HH+VTN(N+2)
                  VTN(N+3) = CO*EPES*HH/VP(N1+3)+VTN(N+3)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N)   = VTN(N)/4
            VTN(N+1) = VTN(N+1)/4
            VTN(N+2) = VTN(N+2)/4
            VTN(N+3) = VTN(N+3)/4
   32    CONTINUE
   30 CONTINUE
C
      VFU = VFU/4+VPI
      DO 41 I = 1,NIN
         VTP   = VTN(I)*WEI(K)
         VC(I) = VC(I)+(VZ(K)-VFU)*VTP
         DO 42 J = 1,I
            RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42    CONTINUE
   41 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)*VPES(3)*VP(THREE)
         DO 71 I=1,NC
            DO 72 J=4,7
               INK     = J+4*(I-1)
               INK1    = INK-2
               VMO     = VC(INK1)*VPES(J)
               VP(INK) = VMO+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM=0.
            DO 81 I = 1,NP
               VG=0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = .5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                          EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                          (1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
C
               VG  = VG/4.+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIPF(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano fisso
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VFU, VG
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EEE
      REAL         DIX, DIY
      INTEGER      N, NIN, N1, NCD
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      NIN = NC+1
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VTN(1)  = 1.
      DO 30  K=1,NP
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = .5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                  END IF
                  VTN(L+1) = EP*HH+VTN(L+1)
   34          CONTINUE
   33       CONTINUE
            VTN(L+1) = VTN(L+1)/4
   32    CONTINUE
C
         DO 41 I = 1,NIN
            VC(I) = VC(I)+VZ(K)*VTN(I)
            DO 42 J = 1,I
               RMT(I,J)=RMT(I,J)+VTN(I)*VTN(J)
   42       CONTINUE
   41    CONTINUE
   30 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)
         DO 71 I=1,NC
            VP(I*4) = VC(I+1)
   71    CONTINUE
C
         SQM=0.
         DO 81 I = 1,NP
            VG=0
            DO 82 J = 1,NC
               N = J*4
               DO 83 II = 1,MX(I)
                  XT = .5*PAS(II,MX(I))+IVX(I)
                  DO 84 JJ=1,MY(I)
                     YT  = .5*PAS(JJ,MY(I))+IVY(I)
                     EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                     (VP(N+3)**2)
                     IF (BETA.LE.0.) THEN
                        VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                       EXP(-EEE*4*ALOG(2.))
                     ELSE
                        VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                       (1+EEE)**(-BETA)
                     END IF
   84             CONTINUE
   83          CONTINUE
   82       CONTINUE
C
            VG  = VG/4.+VP(THREE)
            SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81    CONTINUE
         SQM = SQM/(NP-NIN)
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIPV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano orizzontale
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE, TWO
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VG, VTP
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK, INK1
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      TWO = 2
      NIN = NC*2+1
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VPI     = VP(THREE)
      VTN(3)  = 1.
      DO 30 K=1,NP
         VTN(1) = IVX(K)
         VTN(2) = IVY(K)
         VPI    = VP(1)*IVX(K)+VP(TWO)*IVY(K)+VP(THREE)
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N  = L*2                              ! addition for initialization
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = 0.5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N)   = EP*HH+VTN(N)
                  VTN(N+1) = CO*EPES*HH+VP(N1+3)+VTN(N+1)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N)   = VTN(N)/4
            VTN(N+1) = VTN(N+1)/4
   32    CONTINUE
C
         VFU = VFU/4+VPI
         DO 41 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
            DO 42 J = 1,I
               RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42       CONTINUE
   41    CONTINUE
   30 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
         DO 71 I=1,NC
            DO 72 J=4,7
               INK     = J+4*(I-1)
               INK1    = INK/2
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM=0.
            DO 81 I = 1,NP
               VG=0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = 0.5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                          EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                          (1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
C
               VG  = VG/4.+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano orizzontale
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE, TWO
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VMO, VG, VTP
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      TWO = 2
      NIN = NC*4+3
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VTN(3)  = 1.
      DO 30 K=1,NP
         VTN(1) = IVX(K)
         VTN(2) = IVY(K)
         VPI    = VP(1)*IVX(K)+VP(TWO)*IVY(K)+VP(THREE)
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = 0.5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N1+1) = CO*DIX*HH+VTN(N1+1)
                  VTN(N1+2) = CO*DIY*HH+VTN(N1+2)
                  VTN(N1+3) = CO*EPES*HH/VP(N1+3)+VTN(N1+3)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N1)   = VTN(N1)/4
            VTN(N1+1) = VTN(N1+1)/4
            VTN(N1+2) = VTN(N1+2)/4
            VTN(N1+3) = VTN(N1+3)/4
   32    CONTINUE
   30 CONTINUE
C
      VFU = VFU/4+VPI
      DO 41 I = 1,NIN
         VTP   = VTN(I)*WEI(K)
         VC(I) = VC(I)+(VZ(K)-VFU)*VTP
         DO 42 J = 1,I
            RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42    CONTINUE
   41 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         DO 70 I = 1,3
            VP(I) = VC(I)*VPES(I)*VP(I)
   70    CONTINUE
         DO 71 I=1,NC
            DO 72 J=4,7
               INK     = J+4*(I-1)
               VMO     = VC(INK)*VPES(J)
               VP(INK) = VMO+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM=0.
            DO 81 I = 1,NP
               VG=0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = 0.5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                          EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                          (1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
C
               VG  = VG/4.+VP(1)*IVX(I)+VP(TWO)*IVY(I)+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIH(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano fisso
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VFU, VG
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EEE
      REAL         DIX, DIY
      INTEGER      N, NIN, N1, NCD
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      NIN = NC
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      DO 30  K=1,NP
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = .5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                  END IF
                  VTN(L) = EP*HH+VTN(L)
   34          CONTINUE
   33       CONTINUE
            VTN(L) = VTN(L)/4
   32    CONTINUE
C
         DO 41 I = 1,NIN
            VC(I) = VC(I)+VZ(K)*VTN(I)*WEI(K)
            DO 42 J = 1,I
               RMT(I,J)=RMT(I,J)+VTN(I)*VTN(J)*WEI(K)
   42       CONTINUE
   41    CONTINUE
   30 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)
         DO 71 I=1,NC
            VP(I*4) = VC(I)
   71    CONTINUE
C
         SQM=0.
         DO 81 I = 1,NP
            VG=0
            DO 82 J = 1,NC
               N = J*4
               DO 83 II = 1,MX(I)
                  XT = .5*PAS(II,MX(I))+IVX(I)
                  DO 84 JJ=1,MY(I)
                     YT  = .5*PAS(JJ,MY(I))+IVY(I)
                     EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                     (VP(N+3)**2)
                     IF (BETA.LE.0.) THEN
                        VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                       EXP(-EEE*4*ALOG(2.))
                     ELSE
                        VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                       (1+EEE)**(-BETA)
                     END IF
   84             CONTINUE
   83          CONTINUE
   82       CONTINUE
C
            VG  = VG/4.+VP(THREE)
            SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81    CONTINUE
         SQM = SQM/(NP-NIN)
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END


      SUBROUTINE ELMIX(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                 SIG,MX,MY)
C+++
C.PURPOSE: MOFFAT o Gauss sigma fisso - piano orizzontale
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1)
      INTEGER      IVY(1)
      REAL         VZ(1)
      INTEGER      NP
      REAL         VP(1)
      REAL         FS
      REAL         VPES(7)
      INTEGER      NC
      REAL         BETA
      REAL         SQM
      INTEGER      LG, THREE, TWO
      REAL         WEI(1)
      REAL         SIG
      INTEGER      MX(1),MY(1)
C
      REAL         RMT
      REAL         VTN, VC, VPI, VFU, VMO, VG, VTP
      REAL         XT, YT
      REAL         HH
      REAL         EPES, EP, EPE, EP2, EEE
      REAL         DIX, DIY
      REAL         CO
      INTEGER      N, NIN, N1, NCD
      INTEGER      INK, INK1
      INTEGER      I, J, K, L, LL
      INTEGER      II, JJ
      REAL         FK(NST)
      REAL         PAS(6,6),H(6,6)
C
      COMMON       /SUFR/RMT(NIND,NIND),VTN(NIND),VC(NIND)
      DATA  PAS/
     2   0.0,       0.0,       0.0,      0.0,      0.0,      0.0,     
     2  -0.577350,  0.577350,  0.0,      0.0,      0.0,      0.0,
     *  -0.774597,  0.0,       0.774597, 0.0,      0.0,      0.0,
     *  -0.861136, -0.339981,  0.339981, 0.861136, 0.0,      0.0,
     *  -0.906179, -0.538469,  0.0,      0.538469, 0.906179, 0.0,
     *  -0.932469, -0.661209, -0.238619, 0.238619, 0.661209, 0.932469/
      DATA   H/
     2  2.0,       0.0,       0.0,      0.0,      0.0,      0.0,
     3  1.0,       1.0,       0.0,      0.0,      0.0,      0.0,
     4  0.555555,  0.888889,  0.555555, 0.0,      0.0,      0.0,
     5  0.347855,  0.652145,  0.652145, 0.347855, 0.0,      0.0,
     6  0.236927,  0.478629,  0.568889, 0.478629, 0.236927, 0.0,
     7  0.171324,  0.360762,  0.467914, 0.467914, 0.360762, 0.171324/
C
C ***
      THREE = 3
      TWO = 2
      NIN = NC*3+3
      DO 10 I=1,NIN
         VC(I)=0
         DO 11 J=1,NIN
            RMT(I,J)=0
  11     CONTINUE
  10  CONTINUE
C
      DO 20 I=1,NC
         IF (BETA.LE.0.) THEN
            FK(I) = -4*ALOG(2.)/(VP(4*I+3)**2)
         ELSE
            FK(I) = 1./(VP(4*I+3)**2)
         END IF
   20 CONTINUE
C
      VTN(3)  = 1.
      DO 30 K=1,NP
         VTN(1) = IVX(K)
         VTN(2) = IVY(K)
         VPI    = VP(1)*IVX(K)+VP(TWO)*IVY(K)+VP(THREE)
         VFU=0.
         DO 31 LL=2,NIN
            VTN(LL)=0.
   31    CONTINUE
C
         DO 32 L=1,NC
            N  = L*3 + 1
            N1 = L*4
            DO 33 I = 1,MX(K)
               XT = 0.5*PAS(I,MX(K))+IVX(K)
               DO 34 J = 1,MY(K)
                  HH   = H(I,MX(K))*H(J,MY(K))
                  YT   = 0.5*PAS(J,MY(K))+IVY(K)
                  DIX  = XT-VP(N1+1)
                  DIY  = YT-VP(N1+2)
                  EPES = DIX*DIX+DIY*DIY
C
                  IF (BETA.LE.0) THEN
                     EP = EXP(EPES*FK(L))
                     CO = -VP(N1)*EP*2.*FK(L)
                  ELSE
                     EPE = 1.+EPES*FK(L)
                     EP  = EPE**(-BETA)
                     EP2 = EPE**(-BETA-1.)
                     CO  = VP(N1)*BETA*EP2*2*FK(L)
                  END IF
                  VTN(N)   = EP*HH+VTN(N)
                  VTN(N+1) = CO*DIX*HH+VTN(N+1)
                  VTN(N+2) = CO*DIY*HH+VTN(N+2)
                  VFU      = VP(N1)*EP*HH+VFU
   34          CONTINUE
   33       CONTINUE
            VTN(N)   = VTN(N)/4
            VTN(N+1) = VTN(N+1)/4
            VTN(N+2) = VTN(N+2)/4
   32    CONTINUE
C
         VFU = VFU/4+VPI
         DO 41 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
            DO 42 J = 1,I
               RMT(I,J)=RMT(I,J)+VTN(J)*VTP
   42      CONTINUE
   41    CONTINUE
   30 CONTINUE
C
      DO 51 I = 2,NIN
         DO 52 J=1,I-1
            RMT(J,I)=RMT(I,J)
   52    CONTINUE
   51 CONTINUE
C
      DO 61 I=1,NIN
         RMT(I,I)=RMT(I,I)*(1+FS**2)
   61 CONTINUE
C
      NCD= NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         DO 70 I = 1,3
            VP(I) = VC(I)*VPES(I)*VP(I)
   70    CONTINUE
         DO 71 I=1,NC
            DO 72 J=4,6
               INK     = J+4*(I-1)
               INK1    = J+3*(I-1)
               VMO     = VC(INK1)*VPES(J)
               VP(INK) = VMO+VP(INK)
               IF (ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   72       CONTINUE
   71    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM=0.
            DO 81 I = 1,NP
               VG=0
               DO 82 J = 1,NC
                  N = J*4
                  DO 83 II = 1,MX(I)
                     XT = 0.5*PAS(II,MX(I))+IVX(I)
                     DO 84 JJ=1,MY(I)
                        YT  = .5*PAS(JJ,MY(I))+IVY(I)
                        EEE = ((VP(N+1)-XT)**2+(VP(N+2)-YT)**2)/
     2                        (VP(N+3)**2)
                        IF (BETA.LE.0.) THEN
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     *                          EXP(-EEE*4*ALOG(2.))
                        ELSE
                           VG = VG+H(II,MX(I))*H(JJ,MY(I))*VP(N)*
     2                          (1+EEE)**(-BETA)
                        END IF
   84                CONTINUE
   83             CONTINUE
   82          CONTINUE
C
               VG  = VG/4.+VP(1)*IVX(I)+VP(TWO)*IVY(I)+VP(THREE)
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   81       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
      IF (NCD.LT.1) LG=1
      RETURN
      END

