C
C Single-purpose program generating file 'cv-src.dat' containing the C initial values of hypocentral coordinates and times for linearized C simultaneous inversion C INTEGER LU1,LU2,LU3,LU4,LU5 PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5) CHARACTER*4 RECNAM(14),NAME CHARACTER*12 EVENT REAL REC(3,14),TP(13),TS(13),A(10),E(16) C VSVP0=0. VP0=0. VSVP0=SQRT(1./3.) VP0=18034. C C Receiver coordinates: OPEN(LU1,FILE='cv-stn.dat',STATUS='OLD') NREC=0 11 CONTINUE READ(LU1,'(A4,F7.1,F8.1,F8.1)',END=19) * RECNAM(NREC+1),(REC(I,NREC+1),I=1,3) NREC=NREC+1 GO TO 11 19 CONTINUE CLOSE(LU1) C DPMIN= 999999999. DPMAX=-999999999. DPRMS= 0. NPRMS= 0 RPMIN= 999999999. RPMAX=-999999999. DSMIN= 999999999. DSMAX=-999999999. DSRMS= 0. NSRMS= 0 RSMIN= 999999999. RSMAX=-999999999. X1MIN= 999999999. X1MAX=-999999999. X2MIN= 999999999. X2MAX=-999999999. X3MIN= 999999999. X3MAX=-999999999. TMIN = 999999999. TMAX =-999999999. VPMIN= 999999999. VPMAX=-999999999. VPAVE= 0. VSVPMI= 999999999. VSVPMA=-999999999. VSVPAV= 0. NAVE=0 C C Loop over events: OPEN(LU1,FILE='cv-arr.dat',STATUS='OLD') OPEN(LU2,FILE='cv-abs.out') OPEN(LU3,FILE='cv-rel.out') OPEN(LU4,FILE='cv-loc.out') OPEN(LU5,FILE='cv-eigen.out') 20 CONTINUE C C Reading the event: READ(LU1,'(A12)',END=88) EVENT DO 21 IREC=1,NREC TP(IREC)=0. TS(IREC)=0. 21 CONTINUE DO 24 I=1,NREC READ(LU1,*) NAME,TP0,TS0 DO 22 IREC=1,NREC IF(RECNAM(IREC).EQ.NAME) THEN TP(IREC)=TP0 TS(IREC)=TS0 GO TO 23 END IF 22 CONTINUE 23 CONTINUE 24 CONTINUE C C Ratio VS/VP: TP2=0. TS2=0. DO 32 IREC=1,NREC IF(TP(IREC).NE.0..AND.TS(IREC).NE.0.) THEN DO 31 I=1,IREC-1 IF(TP(I).NE.0..AND.TS(I).NE.0.) THEN IF((TP(IREC)-TP(I))*(TS(IREC)-TS(I)).GT.0.) THEN TP2=TP2+(TP(IREC)-TP(I))**2 TS2=TS2+(TS(IREC)-TS(I))**2 END IF END IF 31 CONTINUE END IF 32 CONTINUE IF(TS2.LE.0.) THEN WRITE(*,'(9A)') ' Event ',EVENT,': Less than 2 P-S time diff.' GO TO 81 END IF VSVP=SQRT(TP2/TS2) IF(VSVP0.NE.0.) THEN VSVP=VSVP0 END IF C C Hypocentral time: N=0 T=0. DO 41 IREC=1,NREC IF(TP(IREC).NE.0..AND.TS(IREC).NE.0.) THEN N=N+1 T=TP(IREC)-VSVP*TS(IREC)+T END IF 41 CONTINUE T=T/(1.-VSVP) T=T/FLOAT(N) C C Average quadratic equation C C0-C1*X1-C2*X2-C3*X3-C4*X4+0.5*(X1*X1+X2*X2+X3*X3): N=0 C0=0. C1=0. C2=0. C3=0. C4=0. DO 51 IREC=1,NREC IF(TP(IREC).NE.0..OR.TS(IREC).NE.0.) THEN N=N+1 C0=C0+REC(1,IREC)**2+REC(2,IREC)**2+REC(3,IREC)**2 C1=C1+REC(1,IREC) C2=C2+REC(2,IREC) C3=C3+REC(3,IREC) IF(TP(IREC).NE.0..AND.TS(IREC).NE.0.) THEN C4=C4+0.5*((TP(IREC)-T)**2+(VSVP*(TS(IREC)-T))**2) ELSE IF(TP(IREC).NE.0.) THEN C4=C4+(TP(IREC)-T)**2 ELSE IF(TS(IREC).NE.0.) THEN C4=C4+(VSVP*(TS(IREC)-T))**2 ELSE WRITE(*,'(9A)') ' ERROR C4' STOP END IF END IF 51 CONTINUE IF(N.LT.4) THEN WRITE(*,'(9A)') ' Event ',EVENT,': Less than four receivers' GO TO 81 END IF C0=C0/2. C4=C4/2. C0=C0/FLOAT(N) C1=C1/FLOAT(N) C2=C2/FLOAT(N) C3=C3/FLOAT(N) C4=C4/FLOAT(N) C C Linear equations (differences of quadratic equations): N=0 DO 61 I=1,10 A(I)=0. 61 CONTINUE R1=0. R2=0. R3=0. R4=0. VP2=0. DO 62 IREC=1,NREC IF(TP(IREC).NE.0..OR.TS(IREC).NE.0.) THEN C One equation: A1*X1+A2*X2+A3*X3+A4*X4=A0 N=N+1 A0=REC(1,IREC)**2+REC(2,IREC)**2+REC(3,IREC)**2 A1=REC(1,IREC)-C1 A2=REC(2,IREC)-C2 A3=REC(3,IREC)-C3 IF(TP(IREC).NE.0..AND.TS(IREC).NE.0.) THEN A4=0.5*((TP(IREC)-T)**2+(VSVP*(TS(IREC)-T))**2) ELSE IF(TP(IREC).NE.0.) THEN A4=(TP(IREC)-T)**2 ELSE IF(TS(IREC).NE.0.) THEN A4=(VSVP*(TS(IREC)-T))**2 ELSE WRITE(*,'(9A)') ' ERROR A4' STOP END IF A0=0.5*A0-C0 A4=0.5*A4-C4 AUX=SQRT(A1*A1+A2*A2+A3*A3) A0=A0/AUX A1=A1/AUX A2=A2/AUX A3=A3/AUX A4=A4/AUX C System of equations: Aij*Xj=Ri A( 1)=A( 1)+A1*A1 A( 2)=A( 2)+A1*A2 A( 3)=A( 3)+A2*A2 A( 4)=A( 4)+A1*A3 A( 5)=A( 5)+A2*A3 A( 6)=A( 6)+A3*A3 A( 7)=A( 7)+A1*A4 A( 8)=A( 8)+A2*A4 A( 9)=A( 9)+A3*A4 A(10)=A(10)+A4*A4 R1=R1+A0*A1 R2=R2+A0*A2 R3=R3+A0*A3 R4=R4+A0*A4 VP2=VP2+A4 END IF 62 CONTINUE VP2=FLOAT(N)/VP2 VP2=VP2/SQRT(3.) IF(VP0.NE.0.) THEN VP2=VP0**2 END IF C Substitution: X4=VP2*X4 A( 7)=A( 7)*VP2 A( 8)=A( 8)*VP2 A( 9)=A( 9)*VP2 A(10)=A(10)*VP2*VP2 R4=R4*VP2 C4=C4*VP2 C C Unknown velocity: IF(VP0.EQ.0.) THEN CALL EIGEN(A,E,4,0) WRITE(LU5,'(3A,4F11.6)') 'Event ',EVENT,':',A(1),A(3),A(6),A(10) * WRITE(LU5,'(19X,4F11.6)') E( 1),E( 5),E( 9),E(13) * WRITE(LU5,'(19X,4F11.6)') E( 2),E( 6),E(10),E(14) * WRITE(LU5,'(19X,4F11.6)') E( 3),E( 7),E(11),E(15) * WRITE(LU5,'(19X,4F11.6)') E( 4),E( 8),E(12),E(16) C Substitution: X1=E( 1)*B1+E( 5)*B2+E( 9)*B3+E(13)*B4 C X2=E( 2)*B1+E( 6)*B2+E(10)*B3+E(14)*B4 C X3=E( 3)*B1+E( 7)*B2+E(11)*B3+E(15)*B4 C X4=E( 4)*B1+E( 8)*B2+E(12)*B3+E(16)*B4 B1=E( 1)*R1+E( 2)*R2+E( 3)*R3+E( 4)*R4 B2=E( 5)*R1+E( 6)*R2+E( 7)*R3+E( 8)*R4 B3=E( 9)*R1+E(10)*R2+E(11)*R3+E(12)*R4 B4=E(13)*R1+E(14)*R2+E(15)*R3+E(16)*R4 B1=B1/A( 1) B2=B2/A( 3) B3=B3/A( 6) C B4=B4/A(10) --- not resolved C Converting C0-C1*X1-C2*X2-C3*X3-C4*X4+0.5*(X1*X1+X2*X2+X3*X3)=0 C into C0-Y1*B1-Y2*B2-Y3*B3-Y4*B4+0.5*(X1*X1+X2*X2+X3*X3)=0 Y1=E( 1)*C1+E( 2)*C2+E( 3)*C3+E( 4)*C4 Y2=E( 5)*C1+E( 6)*C2+E( 7)*C3+E( 8)*C4 Y3=E( 9)*C1+E(10)*C2+E(11)*C3+E(12)*C4 Y4=E(13)*C1+E(14)*C2+E(15)*C3+E(16)*C4 C Converting C0-Y1*B1-Y2*B2-Y3*B3-Y4*B4+0.5*(X1*X1+X2*X2+X3*X3)=0 C into YY-Y4*B4+0.5*B4*B4=0: X1=E( 1)*B1+E( 5)*B2+E( 9)*B3 X2=E( 2)*B1+E( 6)*B2+E(10)*B3 X3=E( 3)*B1+E( 7)*B2+E(11)*B3 YY=C0-Y1*B1-Y2*B2-Y3*B3+0.5*(X1*X1+X2*X2+X3*X3) Y4=Y4-E(13)*X1-E(14)*X2-E(15)*X3 AUX=E(13)**2+E(14)**2+E(15)**2 YY=YY/AUX Y4=Y4/AUX C Solving equation YY-Y4*B4+0.5*B4*B4=0 C (Y4-B4)**2=Y4*Y4-2*YY IF(Y4*Y4-2.*YY.GE.0.) THEN B4=Y4-SQRT(Y4*Y4-2.*YY) * WRITE(LU5,'(19X,4F11.6)') B1,B2,B3,B4 X1=E( 1)*B1+E( 5)*B2+E( 9)*B3+E(13)*B4 X2=E( 2)*B1+E( 6)*B2+E(10)*B3+E(14)*B4 X3=E( 3)*B1+E( 7)*B2+E(11)*B3+E(15)*B4 X4=E( 4)*B1+E( 8)*B2+E(12)*B3+E(16)*B4 VP=SQRT(X4*VP2) VS=VSVP*VP WRITE(LU4,'(3A,6F11.3)') * 'Event ',EVENT,':',X1,X2,X3,T,VP,1./VSVP ELSE B4=Y4 * WRITE(LU5,'(19X,4F11.6)') B1,B2,B3,B4 X1=E( 1)*B1+E( 5)*B2+E( 9)*B3+E(13)*B4 X2=E( 2)*B1+E( 6)*B2+E(10)*B3+E(14)*B4 X3=E( 3)*B1+E( 7)*B2+E(11)*B3+E(15)*B4 X4=E( 4)*B1+E( 8)*B2+E(12)*B3+E(16)*B4 VP=SQRT(X4*VP2) VS=VSVP*VP WRITE(LU4,'(3A,6F11.3)') * 'Event*',EVENT,':',X1,X2,X3,T,VP,1./VSVP END IF C C Known velocity: ELSE R1=R1-A(7) R2=R2-A(8) R3=R3-A(9) CALL EIGEN(A,E,3,0) WRITE(LU5,'(3A,4F11.6)') 'Event ',EVENT,':',A(1),A(3),A(6) C Substitution: X1=E( 1)*B1+E( 4)*B2+E( 7)*B3 C X2=E( 2)*B1+E( 5)*B2+E( 8)*B3 C X3=E( 3)*B1+E( 6)*B2+E( 9)*B3 C X4=1. B1=E( 1)*R1+E( 2)*R2+E( 3)*R3 B2=E( 4)*R1+E( 5)*R2+E( 6)*R3 B3=E( 7)*R1+E( 8)*R2+E( 9)*R3 B1=B1/A( 1) B2=B2/A( 3) C B3=B3/A( 6) --- not resolved C Converting C0-C1*X1-C2*X2-C3*X3-C4*X4+0.5*(X1*X1+X2*X2+X3*X3)= C into C0-Y1*B1-Y2*B2-Y3*B3+0.5*(X1*X1+X2*X2+X3*X3)=0 C0=C0-C4 Y1=E( 1)*C1+E( 2)*C2+E( 3)*C3 Y2=E( 4)*C1+E( 5)*C2+E( 6)*C3 Y3=E( 7)*C1+E( 8)*C2+E( 9)*C3 C Converting C0-Y1*B1-Y2*B2-Y3*B3+0.5*(X1*X1+X2*X2+X3*X3)=0 C into YY-Y3*B3+0.5*B3*B3=0: X1=E( 1)*B1+E( 4)*B2 X2=E( 2)*B1+E( 5)*B2 X3=E( 3)*B1+E( 6)*B2 YY=C0-Y1*B1-Y2*B2+0.5*(X1*X1+X2*X2+X3*X3) Y3=Y3-E( 7)*X1-E( 8)*X2-E( 9)*X3 AUX=E( 7)**2+E( 8)**2+E( 9)**2 YY=YY/AUX Y3=Y3/AUX C Solving equation YY-Y3*B3+0.5*B3*B3=0 C (Y3-B3)**2=Y3*Y3-2*YY IF(Y3*Y3-2.*YY.GE.0.) THEN B3=Y3-SQRT(Y3*Y3-2.*YY) * WRITE(LU5,'(19X,4F11.6)') B1,B2,B3 X1=E( 1)*B1+E( 4)*B2+E( 7)*B3 X2=E( 2)*B1+E( 5)*B2+E( 8)*B3 X3=E( 3)*B1+E( 6)*B2+E( 9)*B3 X4=1. VP=SQRT(X4*VP2) VS=VSVP*VP WRITE(LU4,'(3A,6F11.3)') * 'Event ',EVENT,':',X1,X2,X3,T,VP,1./VSVP ELSE B4=Y3 * WRITE(LU5,'(19X,4F11.6)') B1,B2,B3 X1=E( 1)*B1+E( 4)*B2+E( 7)*B3 X2=E( 2)*B1+E( 5)*B2+E( 8)*B3 X3=E( 3)*B1+E( 6)*B2+E( 9)*B3 X4=1. VP=SQRT(X4*VP2) VS=VSVP*VP WRITE(LU4,'(3A,6F11.3)') * 'Event*',EVENT,':',X1,X2,X3,T,VP,1./VSVP END IF END IF C C Residuals: WRITE(LU2,'(3A,6F11.3)') 'Event ',EVENT,':' WRITE(LU3,'(3A,6F11.3)') 'Event ',EVENT,':' DO 71 IREC=1,NREC IF(TP(IREC).NE.0..OR.TS(IREC).NE.0.) THEN S=SQRT((REC(1,IREC)-X1)**2+(REC(2,IREC)-X2)**2 * +(REC(3,IREC)-X3)**2) IF(TP(IREC).NE.0.) THEN DP=TP(IREC)-T-S/VP RP=DP/(TP(IREC)-T) WRITE(LU2,'(19X,F11.6)') DP WRITE(LU3,'(19X,F11.6)') RP DPMIN=AMIN1(DPMIN,DP) DPMAX=AMAX1(DPMAX,DP) DPRMS=DPRMS+DP**2 NPRMS=NPRMS+1 RPMIN=AMIN1(RPMIN,RP) RPMAX=AMAX1(RPMAX,RP) END IF IF(TS(IREC).NE.0.) THEN DS=TS(IREC)-T-S/VS RS=DS/(TS(IREC)-T) WRITE(LU2,'(30X,F11.6)') DS WRITE(LU3,'(30X,F11.6)') RS DSMIN=AMIN1(DSMIN,DS) DSMAX=AMAX1(DSMAX,DS) DSRMS=DSRMS+DS**2 NSRMS=NSRMS+1 RSMIN=AMIN1(RSMIN,RS) RSMAX=AMAX1(RSMAX,RS) END IF END IF 71 CONTINUE C X1MIN=AMIN1(X1MIN,X1) X1MAX=AMAX1(X1MAX,X1) X2MIN=AMIN1(X2MIN,X2) X2MAX=AMAX1(X2MAX,X2) X3MIN=AMIN1(X3MIN,X3) X3MAX=AMAX1(X3MAX,X3) TMIN =AMIN1(TMIN,T) TMAX =AMAX1(TMAX,T) VPMIN=AMIN1(VPMIN,VP) VPMAX=AMAX1(VPMAX,VP) VPAVE=VPAVE+VP VSVPMI=AMIN1(VSVPMI,VSVP) VSVPMA=AMAX1(VSVPMA,VSVP) VSVPAV=VSVPAV+VSVP NAVE=NAVE+1 C 81 CONTINUE GO TO 20 C 88 CONTINUE WRITE(LU2,'(A,6F11.6)') WRITE(LU2,'(A,6F11.6)') 'Minimum: ',DPMIN,DSMIN WRITE(LU2,'(A,6F11.6)') 'Maximum: ',DPMAX,DSMAX DPRMS=SQRT(DPRMS/FLOAT(NPRMS)) DSRMS=SQRT(DSRMS/FLOAT(NSRMS)) WRITE(LU2,'(A,6F11.6)') 'RMS: ',DPRMS,DSRMS WRITE(LU3,'(A,6F11.6)') WRITE(LU3,'(A,6F11.6)') 'Minimum: ',RPMIN,RSMIN WRITE(LU3,'(A,6F11.6)') 'Maximum: ',RPMAX,RSMAX WRITE(LU4,'(A,6F11.3)') WRITE(LU4,'(A,6F11.3)') 'Minimum: ',X1MIN,X2MIN,X3MIN * ,TMIN,VPMIN,1./VSVPMA WRITE(LU4,'(A,6F11.3)') 'Maximum: ',X1MAX,X2MAX,X3MAX * ,TMAX,VPMAX,1./VSVPMI VPAVE=VPAVE/FLOAT(NAVE) VSVPAV=VSVPAV/FLOAT(NAVE) WRITE(LU4,'(A,44X,6F11.3)') 'Average: ' * ,VPAVE,1./VSVPAV CLOSE(LU1) CLOSE(LU2) CLOSE(LU3) CLOSE(LU4) CLOSE(LU5) STOP END C C======================================================================= C INCLUDE 'eigen.for' C eigen.for C C======================================================================= C