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