C @(#)elmor.for	19.1 (ES0-DMD) 02/25/03 13:29: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
      SUBROUTINE ELMRF(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                  SIG)
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), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, INK1, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG	
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
     2                  VC(NIND)
C
C *** start the code
      THREE  = 3
      NIN=NC*3

      DO 10 I=1,NIN
         VFF(I) = 0
         VFZ(I) = 0
         VC(I)  = 0
         DO 11 J = 1,NIN
            RMT(I,J) = 0
   11    CONTINUE
   10 CONTINUE

      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

      VPI    = VP(THREE)
      VTN(1) = 1.

      DO 30 K=1,NP
         VFU = 0.
         DO 31 L = 1,NC
            N    = L*3-2
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY

            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
            VTN(N+1) = CO*DIX
            VTN(N+2) = CO*DIY
            VFU      = VP(N1)*EP+VFU
   31    CONTINUE

         VFU = VFU+VPI
         DO 32 I = 1,NIN
            VTP    = VTN(I)*WEI(K)
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
            DO 36 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   36       CONTINUE
   32    CONTINUE
   30 CONTINUE

      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
 
      DO 50 I=1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 CONTINUE

      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)

      IF (NCD.GT.0) THEN
         DO 60 I = 1,NC
             DO 61 J = 4,6
               INK     = J+4*(I-1)
               INK1    = J+3*(I-1)
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   61        CONTINUE
   60    CONTINUE

         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 63 I = 1,NP
               VG = VP(THREE)
               DO 64 J = 1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
     2                                         (VP(N+3)**2)
                  IF(BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
   64          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   63       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF

      IF(NCD.LT.1) LG=1

      RETURN
      END

      SUBROUTINE ELMRFV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                  SIG)
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), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, INK1, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
     2                  VC(NIND)
C
C *** start the code
      THREE = 3
      NIN=NC*4

      DO 10 I=1,NIN
         VFF(I) = 0
         VFZ(I) = 0
         VC(I)  = 0
         DO 11 J = 1,NIN
            RMT(I,J) = 0
   11    CONTINUE
   10 CONTINUE

      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

      VPI    = VP(THREE)
      VTN(1) = 1.

      DO 30 K=1,NP
         VFU = 0.
         DO 31 L = 1,NC
            N    = L*4-3
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY

            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
            VTN(N+1) = CO*DIX
            VTN(N+2) = CO*DIY
            VTN(N+3) = CO*EPES/VP(N1+3)
            VFU      = VP(N1)*EP+VFU
   31    CONTINUE

         VFU = VFU+VPI
         DO 32 I = 1,NIN
            VTP    = VTN(I)*WEI(K)
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
            DO 36 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   36       CONTINUE
   32    CONTINUE
   30 CONTINUE

      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
 
      DO 50 I=1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 CONTINUE

      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)

      IF (NCD.GT.0) THEN
         DO 60 I = 1,NC
             DO 61 J = 4,7
               INK     = J+4*(I-1)
               INK1    = INK-3
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   61        CONTINUE
   60    CONTINUE

         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 63 I = 1,NP
               VG = VP(THREE)
               DO 64 J = 1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
     2                                         (VP(N+3)**2)
                  IF(BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
   64          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   63       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF

      IF(NCD.LT.1) LG=1

      RETURN
      END

      SUBROUTINE ELMRR(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                  SIG)
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), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, INK1, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
     2                  VC(NIND)
C
C *** start the code
      THREE = 3
      NIN=NC*3+1

      DO 10 I=1,NIN
         VFF(I) = 0
         VFZ(I) = 0
         VC(I)  = 0
         DO 11 J = 1,NIN
            RMT(I,J) = 0
   11    CONTINUE
   10 CONTINUE

      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

      VPI    = VP(THREE)
      VTN(1) = 1.

      DO 30 K=1,NP
         VFU = 0.
         DO 31 L = 1,NC
            N    = L*3-1
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY

            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
            VTN(N+1) = CO*DIX
            VTN(N+2) = CO*DIY
            VFU      = VP(N1)*EP+VFU
   31    CONTINUE

         VFU = VFU+VPI
         DO 32 I = 1,NIN
            VTP    = VTN(I)*WEI(K)
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
            DO 36 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   36       CONTINUE
   32    CONTINUE
   30 CONTINUE

      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
 
      DO 50 I=1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 CONTINUE

      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)

      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
         DO 60 I = 1,NC
             DO 61 J = 4,6
               INK     = J+4*(I-1)
               INK1    = J+3*(I-1)-2
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   61        CONTINUE
   60    CONTINUE

         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 63 I = 1,NP
               VG = VP(THREE)
               DO 64 J = 1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
     2                                         (VP(N+3)**2)
                  IF(BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
   64          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   63       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF

      IF(NCD.LT.1) LG=1

      RETURN
      END

      SUBROUTINE ELMRRV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                  SIG)
C
C     MOFFAT O GAUSS SIGMA VAR. - PIANO ORIZZONTALE
C
C-------------
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER (NST=40)
      PARAMETER (NIND=NST*4+3)
C
      INTEGER      IVX(1), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, INK1, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
     2             VC(NIND)

      THREE = 3
      NIN = NC*4+1
      DO 10 I=1,NIN
         VFF(I)=0
         VFZ(I)=0
         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 L=1,NC
            N    = L*4-2
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY
            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
C
            VTN(N)   = EP
            VTN(N+1) = CO*DIX
            VTN(N+2) = CO*DIY 
            VTN(N+3) = CO*EPES/VP(N1+3)
            VFU      = VP(N1)*EP+VFU
   31    CONTINUE
C
         VFU = VFU+VPI
         DO 32 I=1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
            DO 132 J=1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
  132       CONTINUE
   32    CONTINUE
   30 CONTINUE
C 
      DO 40 I=2,NIN
         DO 41 J=1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
C
      DO 50 I=1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 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 61 I=1,NC
            DO 161 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
  161       CONTINUE
   61    CONTINUE
C
         IF (NCD.GT.0) THEN 
            SQM = 0.
            DO 62 I=1,NP
               VG=VP(THREE)
               DO 162 J=1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
     2                   (VP(N+3)**2)
                  IF (BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
  162          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   62       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
C
      IF (NCD.LT.1) LG=1
      RETURN
      END

      SUBROUTINE ELMRPF(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,SIG)
C+++
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale - posiz. fissa
C---  
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER (NST=40)
      PARAMETER (NIND=4*NST+3)
C
      INTEGER      IVX(1), 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
C
      REAL         VC, VTP
      INTEGER      N1
      INTEGER      NCD
      REAL         EPE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES,EEE
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON    /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),
     2                VFF(NIND),VC(NIND)
C
      THREE = 3
      NIN = NC+1
      DO 10 I = 1,NIN
         VC(I)=0
         DO 20  J = 1,NIN
            RMT(I,J)=0
   20    CONTINUE
   10 CONTINUE
C
      DO 30 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
   30 CONTINUE
C
      VTN(1) = 1.
      DO 40  K = 1,NP
         VFU = 0.
         DO 41 L=1,NC
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY
            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
   41    CONTINUE
C
         DO 42 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+VZ(K)*VTP
            DO 43 J=1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   43       CONTINUE
   42    CONTINUE
   40 CONTINUE
C
      DO 50 I = 2,NIN
         DO 51 J=1,I-1
            RMT(J,I)=RMT(I,J)
   51    CONTINUE
   50 CONTINUE
C
      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)
         DO 60 I=1,NC
            VP(I*4) = VC(I+1)
   60    CONTINUE
         SQM=0.
         DO 70 I=1,NP
            VG = VP(THREE)
            DO 71 J=1,NC
               N   = J*4
               EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
     2                (VP(N+3)**2)
               IF (BETA.LE.0.) THEN
                  VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
               ELSE
                  VG = VG+VP(N)*(1+EEE)**(-BETA)
               END IF
   71       CONTINUE
            SQM=SQM+(VZ(I)-VG)**2*WEI(I)
   70    CONTINUE
         SQM=SQM/(NP-NIN)
      END IF
C
      IF (NCD.LT.1) LG=1
      RETURN
      END

      SUBROUTINE ELMRPV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                  SIG)
C+++
C.PURPOSE: Moffat o Gauss sigma fisso - piano orizzontale - posiz. fissa
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER    (NST=40)
      PARAMETER    (NIND=4*NST+3)
C
      INTEGER      IVX(1), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, INK1, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
     2                  VC(NIND)
C
C *** start the code
      THREE = 3 
      NIN=NC*2+1

      DO 10 I=1,NIN
         VFF(I) = 0
         VFZ(I) = 0
         VC(I)  = 0
         DO 11 J = 1,NIN
            RMT(I,J) = 0
   11    CONTINUE
   10 CONTINUE

      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

      VPI    = VP(THREE)
      VTN(1) = 1.

      DO 30 K=1,NP
         VFU = 0.
         DO 31 L = 1,NC
            N    = L*2
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY

            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
            VTN(N+1) = CO*EPES/VP(N1+3)
            VFU      = VP(N1)*EP+VFU
   31    CONTINUE

         VFU = VFU+VPI
         DO 32 I = 1,NIN
            VTP    = VTN(I)*WEI(K)
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
            DO 36 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   36       CONTINUE
   32    CONTINUE
   30 CONTINUE

      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
 
      DO 50 I=1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 CONTINUE

      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)

      IF (NCD.GT.0) THEN
         VP(THREE) = VC(1)*VPES(3)+VP(THREE)
         DO 60 I = 1,NC
             DO 61 J = 4,7,3
               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
   61        CONTINUE
   60    CONTINUE

         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 63 I = 1,NP
               VG = VP(THREE)
               DO 64 J = 1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
     2                                         (VP(N+3)**2)
                  IF(BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
   64          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   63       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF

      IF(NCD.LT.1) LG=1

      RETURN
      END

      SUBROUTINE ELMRV(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,SIG)
C+++
C.PURPOSE: Moffat o Gauss sigma var. - piano inclinato
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER (NST=40)
      PARAMETER (NIND=4*NST+3)
C
      INTEGER      IVX(1), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON    /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),
     2                 VFF(NIND),VC(NIND)
C
      TWO   = 2
      THREE = 3
      NIN = NC*4+3
      DO 10 I = 1,NIN
         VFF(I) = 0
         VFZ(I) = 0
         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(THREE) = 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 L = 1,NC
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY
            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)   = EP
            VTN(N1+1) = CO*DIX
            VTN(N1+2) = CO*DIY
            VTN(N1+3) = CO*EPES/VP(N1+3)
            VFU       = VP(N1)*EP+VFU
   31    CONTINUE
C
         VFU = VFU+VPI
         DO 32 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+(VZ(K)-VFU)*VTP
            DO 33 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   33       CONTINUE
   32     CONTINUE
   30 CONTINUE
C
      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE

      DO 50 I = 1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 CONTINUE
C
      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         DO 60 I = 1,3
            VP(I) = VC(I)*VPES(I)+VP(I)
   60    CONTINUE
C
         DO 70 I = 1,NC
            DO 71 J = 4,7
              INK     = J+4*(I-1)
               VP(INK) = VC(INK)*VPES(J)+VP(INK)
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   71       CONTINUE
   70    CONTINUE
C
         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 80 I = 1,NP
               VG = VP(1)*IVX(I)+VP(TWO)*IVY(I)+VP(THREE)
               DO 81 J = 1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
     2                   (VP(N+3)**2)
                  IF (BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
   81          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   80       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF
C
      IF (NCD.LT.1) LG=1
      RETURN
      END

      SUBROUTINE ELMRH(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,SIG)
C+++
C.PURPOSE: Moffat o Gauss sigma fisso - piano fisso - posiz. fissa
C---
      IMPLICIT     NONE
      INTEGER      NST
      INTEGER      NIND
      PARAMETER  (NST=40)
      PARAMETER  (NIND=4*NST+3)
C
      INTEGER      IVX(1), 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
C
      REAL         VC, VTP
      INTEGER      N1
      INTEGER      NCD
      REAL         EPE, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON    /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),
     2                VFF(NIND),VC(NIND)
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 L = 1,NC
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY
            IF (BETA.LE.0) THEN
               EP = EXP(EPES*FK(L))
            ELSE
               EPE = 1.+EPES*FK(L)
               EP  = EPE**(-BETA)
            END IF
            VTN(L) = EP
   31    CONTINUE
C
         DO 32 I = 1,NIN
            VTP   = VTN(I)*WEI(K)
            VC(I) = VC(I)+VZ(K)*VTP
            DO 33 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   33       CONTINUE
   32    CONTINUE
   30 CONTINUE
C
      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
C
      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)
      IF (NCD.GT.0) THEN
         DO 50 I = 1,NC
            VP(I*4) = VC(I)
   50    CONTINUE
C
         SQM = 0.
         DO 60 I = 1,NP
            VG = VP(THREE)   
            DO 61 J = 1,NC
               N   = J*4  
               EEE = ((VP(N+1)-IVX(I))**2+(VP(N+2)-IVY(I))**2)/
     2                (VP(N+3)**2)
               IF (BETA.LE.0.) THEN
                  VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
               ELSE
                  VG = VG+VP(N)*(1+EEE)**(-BETA)
               END IF
   61       CONTINUE
            SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   60    CONTINUE
         SQM = SQM/(NP-NIN)
      END IF
C
      IF (NCD.LT.1) LG=1
      RETURN
      END

      SUBROUTINE ELMRX(IVX,IVY,VZ,NP,VP,FS,VPES,NC,BETA,SQM,LG,WEI,
     2                  SIG)
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), 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
C
      REAL         VC, VTP, VPI
      REAL         CO
      INTEGER      N1
      INTEGER      INK, INK1, NCD
      REAL         EPE, EP2, EEE
      REAL         RMT
      REAL         VTN, VFZ, VFF, VFU, VG
      REAL         DIX, DIY
      REAL         EP, EPES
      INTEGER      NIN
      INTEGER      I, J, K, L, N
      REAL         FK(NST)
C
      COMMON      /SUFR/RMT(NIND,NIND),VTN(NIND),VFZ(NIND),VFF(NIND),
     2                  VC(NIND)
C
C *** start the code
      TWO   = 2
      THREE = 3
      NIN=NC*3+3

      DO 10 I=1,NIN
         VFF(I) = 0
         VFZ(I) = 0
         VC(I)  = 0
         DO 11 J = 1,NIN
            RMT(I,J) = 0
   11    CONTINUE
   10 CONTINUE

      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

      VTN(THREE) = 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 L = 1,NC
            N    = L*3+1
            N1   = L*4
            DIX  = IVX(K)-VP(N1+1)
            DIY  = IVY(K)-VP(N1+2)
            EPES = DIX*DIX+DIY*DIY

            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
            VTN(N+1) = CO*DIX
            VTN(N+2) = CO*DIY
            VFU      = VP(N1)*EP+VFU
   31    CONTINUE

         VFU = VFU+VPI
         DO 32 I = 1,NIN
            VTP    = VTN(I)*WEI(K)
            VC(I)  = VC(I)+(VZ(K)-VFU)*VTP
            DO 36 J = 1,I
               RMT(I,J) = RMT(I,J)+VTN(J)*VTP
   36       CONTINUE
   32    CONTINUE
   30 CONTINUE

      DO 40 I = 2,NIN
         DO 41 J = 1,I-1
            RMT(J,I) = RMT(I,J)
   41    CONTINUE
   40 CONTINUE
 
      DO 50 I=1,NIN
         RMT(I,I) = RMT(I,I)*(1+FS**2)
   50 CONTINUE

      NCD = NIND
      CALL LISIB(RMT,VC,NIN,NCD,SIG)

      IF (NCD.GT.0) THEN
         DO 59 I = 1,3
            VP(I) = VC(I)*VPES(I)+VP(I) 
   59    CONTINUE
         DO 60 I = 1,NC
            DO 61 J = 4,6
               INK     = J+4*(I-1)
               INK1    = J+3*(I-1)
               VP(INK) = VC(INK1)*VPES(J)+VP(INK)
               IF(ABS(VP(INK)).GT.1000..AND.J.NE.4) NCD=-1
   61        CONTINUE
   60    CONTINUE

         IF (NCD.GT.0) THEN
            SQM = 0.
            DO 63 I = 1,NP
               VG = VP(1)*IVX(I)+VP(TWO)*IVY(I)+VP(THREE)
               DO 64 J = 1,NC
                  N   = J*4
                  EEE = ((VP(N+1)-IVX(I))**2 + (VP(N+2)-IVY(I))**2)/
     2                                         (VP(N+3)**2)
                  IF(BETA.LE.0.) THEN
                     VG = VG+VP(N)*EXP(-EEE*4*ALOG(2.))
                  ELSE
                     VG = VG+VP(N)*(1+EEE)**(-BETA)
                  END IF
   64          CONTINUE
               SQM = SQM+(VZ(I)-VG)**2*WEI(I)
   63       CONTINUE
            SQM = SQM/(NP-NIN)
         END IF
      END IF

      IF(NCD.LT.1) LG=1
      RETURN
      END
