C 
C Program TCCOMP to read three files with propagator matrices U1, U2, U3
C written in the form of file GREENTC,
C to compute matrix DU=inv(U1)*(U2-U3), and to compute its norm
C ||DU||=square root(trace(DUh*DU)), where DUh is Hermitian adjoint
C matrix to the matrix DU.
C
C Version: 5.40
C Date: 2000, February 17
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of a SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Name of the input file with the propagation matrices:
C     GREENTC1='string'... Name of the file with matrix U1.
C         Description of file GREENTC1
C             Default: GREENTC1=' '
C     GREENTC2='string'... Name of the file with matrix U2.
C         Description of file GREENTC2
C             Default: GREENTC2=' '
C     GREENTC='string'... Name of the file with matrix U3.
C         Description of file GREENTC3
C             Default: GREENTC=' '
C Name of the output file:
C     TCCOMP='string'... Name of the output formatted file with
C             matrix DU and its norm.
C             Default: TCCOMP='tccomp.out'
C
C Output formatted file TCCOMP:
C   For each frequency following line:
C     F,REDU11,IMDU11,REDU21,IMDU21,REDU12,IMDU12,REDU22,IMDU22,ADU,TEXT
C   where F is the frequency, REDU11,IMDU11,REDU21,IMDU21,REDU12,IMDU12,
C   REDU22,IMDU22 are real and imaginary parts of the 2x2 matrix DU,
C   ADU is its norm as defined above, and text is a character enabling
C   the output file to be included into PostScript file. Thus the first
C   line terminates with character M, and other lines terminate with
C   character L.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T
C     ERROR ... File error.for.
C     RSEP1,RSEP3T ... File sep.for.
C
C-----------------------------------------------------------------------
C
C     Input and output data files:
      CHARACTER*80 FSEP,FIN1,FIN2,FIN3,FOUT
      CHARACTER*2 TEXT
      INTEGER LU,LU1,LU2,LU3,LU4
      PARAMETER (LU=5,LU1=1,LU2=2,LU3=3,LU4=4)
C     Auxiliary storage locations:
      REAL RDU11,RDU21,RDU12,RDU22,IDU11,IDU21,IDU12,IDU22
      REAL RU11,RU21,RU12,RU22,IU11,IU21,IU12,IU22
      REAL RV11,RV21,RV12,RV22,IV11,IV21,IV12,IV22
      REAL RW11,RW21,RW12,RW22,IW11,IW21,IW12,IW22
      REAL RX11,RX21,RX12,RX22,IX11,IX21,IX12,IX22
      REAL RY11,RY21,RY12,RY22,IY11,IY21,IY12,IY22
      REAL F,RD,ID,ADU
C
C.......................................................................
C
C     Main input data:
      FSEP=' '
      WRITE(*,'(A)') '+TCCOMP: Enter input filename: '
      READ(*,*) FSEP
      IF (FSEP.EQ.' ') THEN
C       TCCOMP-01
        CALL ERROR('TCCOMP-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
      WRITE(*,'(A)') '+TCCOMP: Working...            '
C
C     Reading all the data from file FSEP to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
C
C     Recalling the data:
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3T('GREENTC1',FIN1,' ')
      CALL RSEP3T('GREENTC2',FIN2,' ')
      CALL RSEP3T('GREENTC' ,FIN3,' ')
      CALL RSEP3T('TCCOMP',FOUT,'tccomp.out')
      OPEN(LU1,FILE=FIN1,STATUS='OLD')
      OPEN(LU2,FILE=FIN2,STATUS='OLD')
      OPEN(LU3,FILE=FIN3,STATUS='OLD')
      OPEN(LU4,FILE=FOUT)
      TEXT=' M'
  10  CONTINUE
        READ(LU1,*,END=100) F,RU11,IU11,RU21,IU21,RU12,IU12,RU22,IU22
        READ(LU2,*,END=100) F,RV11,IV11,RV21,IV21,RV12,IV12,RV22,IV22
        READ(LU3,*,END=100) F,RW11,IW11,RW21,IW21,RW12,IW12,RW22,IW22
        RX11=RV11-RW11
        IX11=IV11-IW11
        RX21=RV21-RW21
        IX21=IV21-IW21
        RX12=RV12-RW12
        IX12=IV12-IW12
        RX22=RV22-RW22
        IX22=IV22-IW22
        RD=RU11*RU22-IU11*IU22 - RU21*RU12+IU21*IU12
        ID=RU11*IU22+IU11*RU22 - RU21*IU12-IU21*RU12
        ADU=SQRT(RD*RD+ID*ID)
        RD= RD/ADU
        ID=-ID/ADU
        RY11=RU22
        IY11=IU22
        RY21=-RU21
        IY21=-IU21
        RY12=-RU12
        IY12=-IU12
        RY22=RU11
        IY22=IU11
        RU11=RY11*RD-IY11*ID
        IU11=RY11*ID+IY11*RD
        RU21=RY21*RD-IY21*ID
        IU21=RY21*ID+IY21*RD
        RU12=RY12*RD-IY12*ID
        IU12=RY12*ID+IY12*RD
        RU22=RY22*RD-IY22*ID
        IU22=RY22*ID+IY22*RD
        RDU11=RU11*RX11-IU11*IX11 + RU12*RX21-IU12*IX21
        IDU11=RU11*IX11+IU11*RX11 + RU12*IX21+IU12*RX21
        RDU21=RU21*RX11-IU21*IX11 + RU22*RX21-IU22*IX21
        IDU21=RU21*IX11+IU21*RX11 + RU22*IX21+IU22*RX21
        RDU12=RU11*RX12-IU11*IX12 + RU12*RX22-IU12*IX22
        IDU12=RU11*IX12+IU11*RX12 + RU12*IX22+IU12*RX22
        RDU22=RU21*RX12-IU21*IX12 + RU22*RX22-IU22*IX22
        IDU22=RU21*IX12+IU21*RX12 + RU22*IX22+IU22*RX22
        ADU=RDU11**2+IDU11**2+RDU21**2+IDU21**2+RDU12**2+IDU12**2+
     *      RDU22**2+IDU22**2
        ADU=SQRT(0.5*ADU)
        WRITE(LU4,'(9(F9.6,1X),F9.6,A)')
     *    F,RDU11,IDU11,RDU21,IDU21,RDU12,IDU12,RDU22,IDU22,ADU,TEXT
        TEXT=' L'
      GOTO 10
  100 CONTINUE
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      WRITE(*,'(A)') '+TCCOMP: Done.                  '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
C=======================================================================
C