C @(#)elmoi.for 17.1.1.1 (ESO-DMD) 01/25/02 17:17:55 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