C 
C Program TCGREEN to compute propagator matrices U of equation (39) of
C the paper by Klimes (1999).  The propagator matrices are written
C to the output formatted files in the form of file
C GREEN.
C
C References:
C     Klimes L.(1999): 'Analytical one-way plane-wave solution in the
C     1-D anisotropic "twisted crystal" model'.  Report 8, pp.103-118,
C     Charles University, Prague.
C
C     Bulant, P., & Klimes L.(1999): 'Comparison of ray methods
C     with the exact solution in the 1-D anisotropic
C     "twisted crystal" model'.  Report 8, pp.119-126,
C     Charles University, Prague.
C
C Version: 5.80
C Date: 2004, June 11
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 Data describing the twisted-crystal model:
C     SIN2TH=real ... Second power of sinus theta
C              (equation 15 of the paper by Klimes).
C              Used to specify the default values of EPSILON and V0.
C              Has no meaning if both EPSILON and V0 are specified.
C              Default: SIN2TH=0.
C     GAMMA=real ... Gamma of equation 15 of the paper by Klimes.
C              Used to specify the default values of EPSILON and V0.
C              Has no meaning if both EPSILON and V0 are specified.
C              Default: GAMMA=0.
C     TCK=real ... Uppercase K of equation 9 of the paper by Klimes.
C              Default: TCK=0.
C     EPSILON=real ... Epsilon of equation 9 of the paper by Klimes.
C              Default: EPSILON=-(GAMMA*SIN2TH)/(1.+GAMMA*SIN2TH)
C     A44=real ...  a44 of equation 3 of the paper by Bulant & Klimes.
C              Used to specify the default value of V0, and to specify
C              the reference stress for output file TCSTRESS.
C              Default: A44=0.
C     V0=real ... v0 of equation 8 of the paper by Klimes.
C              Default: V0=SQRT(A44*(1.+GAMMA*SIN2TH))
C     VREF2=real ... second power of reference velocity
C              (see equation 111 of the paper by Klimes).
C              Default: VREF2=V0*V0
C Data describing source and receiver:
C     REC='string'... If non-blank, the name of the file with the names
C              and coordinates of the receiver points.  Only the first
C              point in the file is taken into account. Its name is
C              used in output files. If X3 is not defined, x3 coordinate
C              of the point is used instead.
C           Description of file REC
C              Default: REC=' '
C     SRC='string'... If non-blank, the name of the file with the name
C              and coordinates of the source point.  The name is used in
C              the output files, the coordinates are not considered.
C           Description of file SRC
C              Default: SRC=' '
C     X3=real ... Parameter specifying the source-receiver
C              distance X3=X3receiver-X3source in the direction of the
C              X3 axis. If not specified, coordinates of the source and
C              receiver are taken from files SRC and REC (default).
C              Default: X3=X3receiver-X3source
C Data describing the frequency domain:
C     DF=real ... Frequency step.
C              Default: DF=0.
C     OF=real ... The propagator matrix is calculated for NF
C              frequencies OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF.
C              Default: OF=0.
C     NF=integer ... Number of frequencies.
C              Default: NF=1
C Names of the output files:
C              The output files are not generated if the corresponding
C              name equals ' '.
C     TCGREENE='string'... Name of the output formatted file with the
C              exact propagator matrix U (equation 39) written in the
C              form of the file GREEN.
C              Default: TCGREENE='tcgreene.out'
C     TCGREENW='string'... Name of the output formatted file with the
C              coupling ray theory propagator matrix (equation 39 with
C              coefficients F0 to F3 from equation 106) written in the
C              form of the file GREEN.
C              Default: TCGREENW='tcgreenw.out'
C     TCGREENQ='string'... Name of the output formatted file with the
C              quasi-isotropic propagator matrix (equation 39 with
C              coefficients F0 to F3 from equation 111) written in the
C              form of the file GREEN.
C              Default: TCGREENQ='tcgreenq.out'
C     TCGREENA='string'... Name of the output formatted file with the
C              anisotropic propagator matrix (equation 39 with
C              coefficients F0 to F3 from equation 116) written in the
C              form of the file GREEN.
C              Default: TCGREENA='tcgreena.out'
C     TCGREENI='string'... Name of the output formatted file with the
C              isotropic propagator matrix (equation 39 with
C              coefficients F0 to F3 from equation 121) written in the
C              form of the file GREEN.
C              Default: TCGREENI='tcgreeni.out'
C     TCSTRESS='string'... Name of the output formatted file with the
C              stress components corresponding to the exact propagator
C              matrix at the initial point.  The stress is divided by
C              "reference stress"=2*pi*F*SQRT(A44),
C              where F is frequency.
C              Default: TCSTRESS='tcstress.out'
C Rotation of the model:
C     TCE1=real, TCE2=real, TCE3=real ... Three components of
C              the rotation vector describing the rotation of the model.
C              The model described in the papers referenced above may be
C              further rotated by specifying the rotation vector.
C              The model is rotated along the axis given by the rotation
C              vector, the angle of the rotation in radians is given by
C              the Cartesian length of the vector.
C              Only the exact solution may be calculated in the rotated
C              model, all the parameters TCGREENW(Q,A,I) must equal ' '
C              if the rotation vector is specified.
C              Default: TCE1=TCE2=TCE3=0. (no rotation).
C
C Description of the output files TCGREENE(W,Q,A,I):
C     Formatted files containing the computed propagator matrix.
C     The files have the same format as the file GREEN generated by
C     program GREEN. Strings 'Src' and 'Rec' are written in place
C     of names of the source and the receiver, coordinates of the source
C     are [0.,0.,0.], coordinates of the receiver are [0.,0.,X3],
C     travel time is X3/V0, zeros are written in place of derivatives
C     of the travel time with respect to the coordinates of the receiver
C     and coordinates of the source. Components of the 2x2 propagator
C     matrix are written as corresponding amplitudes, zeros are written
C     instead of the remaining amplitudes.
C     Description of file GREEN
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,RSEP3R,FORM2,LENGTH,
     * EQ39,TCGMAT,TCGTRA
      INTEGER LENGTH
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C             File sep.for.
C     FORM2 ... File forms.for.
C     LENGTH ... File length.for.
C     EQ39,TCGMAT,TCGTRA ... This file.
C
C-----------------------------------------------------------------------
C
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      REAL GREEN(MRAM)
      EQUIVALENCE (GREEN,RAM)
C
C     Input and output data files:
      CHARACTER*80 FSEP,FOUTE,FOUTW,FOUTQ,FOUTA,FOUTI,FOUTS,
     *             FILREC,FILSRC,RECNAM,SRCNAM
      CHARACTER*260 FORMAT
      INTEGER LU
      PARAMETER (LU=1)
      REAL PI
      PARAMETER (PI=3.1415926)
C
C     Auxiliary storage locations:
      INTEGER I1,I2,I,NF,NGREEN
      REAL EPS,S2TH,GAMMA,OF,DF,F,VK,K0,A44,V0,X3,K02,VK2,
     * E1,E2,E3,AE,AE2,COSAE,SINAE,AUX1,AUX2,
     * K02VK2,EPS2,JEPS2,SJEPS2,SQJMEP,SQJPEP,
     * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,RECIT,
     * REJM,IMJM,AUX,REAUX,IMAUX,REPZ,IMPZ,REDZ,IMDZ,
     * REF12,IMF12,REF02,IMF02,AF02,FIF02,AF0,FIF0,
     * REU11,IMU11,REU21,IMU21,REU31,IMU31,
     * REU12,IMU12,REU22,IMU22,REU32,IMU32,
     * REU13,IMU13,REU23,IMU23,REU33,IMU33,REU(3,3),IMU(3,3),
     * REE11,IME11,REE21,IME21,REE31,IME31,
     * REE12,IME12,REE22,IME22,REE32,IME32,
     * REE13,IME13,REE23,IME23,REE33,IME33,REE(3,3),IME(3,3),
     * CFIF02,SFIF02,TT,VR,VR2,R1,R2,R3,
     * REA(3,3),IMA(3,3),REET(3,3),IMET(3,3)
      LOGICAL LROT
      EQUIVALENCE (REU11,REU(1,1)),(IMU11,IMU(1,1)),(REU21,REU(2,1)),
     *            (IMU21,IMU(2,1)),(REU31,REU(3,1)),(IMU31,IMU(3,1)),
     *            (REU12,REU(1,2)),(IMU12,IMU(1,2)),(REU22,REU(2,2)),
     *            (IMU22,IMU(2,2)),(REU32,REU(3,2)),(IMU32,IMU(3,2)),
     *            (REU13,REU(1,3)),(IMU13,IMU(1,3)),(REU23,REU(2,3)),
     *            (IMU23,IMU(2,3)),(REU33,REU(3,3)),(IMU33,IMU(3,3)),
     *            (REE11,REE(1,1)),(IME11,IME(1,1)),(REE21,REE(2,1)),
     *            (IME21,IME(2,1)),(REE31,REE(3,1)),(IME31,IME(3,1)),
     *            (REE12,REE(1,2)),(IME12,IME(1,2)),(REE22,REE(2,2)),
     *            (IME22,IME(2,2)),(REE32,REE(3,2)),(IME32,IME(3,2)),
     *            (REE13,REE(1,3)),(IME13,IME(1,3)),(REE23,REE(2,3)),
     *            (IME23,IME(2,3)),(REE33,REE(3,3)),(IME33,IME(3,3))
      DATA REE,IME/18*0./
C
C-----------------------------------------------------------------------
C
C     Main input data:
      FSEP=' '
      WRITE(*,'(A)') '+TCGREEN: Enter input filename: '
      READ(*,*) FSEP
      IF (FSEP.EQ.' ') THEN
C       TCGREEN-01
        CALL ERROR('TCGREEN-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)') '+TCGREEN: Working...            '
C
C     Reading all the data from file FSEP to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
C
C     Reading the names of the output files:
      CALL RSEP3T('TCGREENE',FOUTE,'tcgreene.out')
      CALL RSEP3T('TCGREENW',FOUTW,'tcgreenw.out')
      CALL RSEP3T('TCGREENQ',FOUTQ,'tcgreenq.out')
      CALL RSEP3T('TCGREENA',FOUTA,'tcgreena.out')
      CALL RSEP3T('TCGREENI',FOUTI,'tcgreeni.out')
      CALL RSEP3T('TCSTRESS',FOUTS,'tcstress.out')
C
C     Recalling the data:
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3R('EPSILON',EPS,999999.)
      CALL RSEP3R('SIN2TH',S2TH,0.)
      CALL RSEP3R('GAMMA',GAMMA,0.)
      CALL RSEP3R('OF',OF,0.)
      CALL RSEP3R('DF',DF,1.)
      CALL RSEP3I('NF',NF,1)
      CALL RSEP3R('TCK',VK,0.)
      CALL RSEP3R('V0',V0,999999.)
      CALL RSEP3R('A44',A44,0.)
      CALL RSEP3R('TCE1',E1,0.)
      CALL RSEP3R('TCE2',E2,0.)
      CALL RSEP3R('TCE3',E3,0.)
      AE2=E1*E1+E2*E2+E3*E3
      IF (AE2.EQ.0.) THEN
        LROT=.FALSE.
      ELSE
        LROT=.TRUE.
        AE=SQRT(AE2)
        COSAE=COS(AE)
        AUX1=(1-COSAE)/AE2
        SINAE=SIN(AE)
        AUX2=SINAE/AE
        REE11=COSAE+E1*E1*AUX1
        REE21=      E2*E1*AUX1                -E3*AUX2
        REE31=      E3*E1*AUX1        +E2*AUX2
        REE12=      E1*E2*AUX1                +E3*AUX2
        REE22=COSAE+E2*E2*AUX1
        REE32=      E3*E2*AUX1-E1*AUX2
        REE13=      E1*E3*AUX1        -E2*AUX2
        REE23=      E2*E3*AUX1+E1*AUX2
        REE33=COSAE+E3*E3*AUX1
        CALL TCGTRA(REE,IME,REET,IMET)
        IF ((FOUTW.NE.' ').OR.(FOUTQ.NE.' ').OR.(FOUTA.NE.' ').OR.
     *      (FOUTI.NE.' ')) THEN
C         TCGREEN-06
          CALL WARN('TCGREEN-06: Only file TCGREENE may be calculated')
C         If the model is to be rotated, only the exact solution may be
C         calculated.  All other output files will not be generated.
          FOUTW=' '
          FOUTQ=' '
          FOUTA=' '
          FOUTI=' '
        ENDIF
      ENDIF
C     Name of the source:
      CALL RSEP3T('SRC',FILSRC,' ')
      IF (FILSRC.NE.' ') THEN
        OPEN(LU,FILE=FILSRC,STATUS='OLD')
        READ(LU,*) (SRCNAM,I=1,20)
        R3=0.
        READ(LU,*) SRCNAM,R1,R2,R3
        CLOSE(LU)
      ENDIF
C     Name of the receiver:
      CALL RSEP3T('REC',FILREC,' ')
      IF (FILREC.NE.' ') THEN
        OPEN(LU,FILE=FILREC,STATUS='OLD')
        READ(LU,*) (RECNAM,I=1,20)
        X3=0.
        READ(LU,*) RECNAM,R1,R2,X3
        CLOSE(LU)
      ENDIF
      R3=X3-R3
      CALL RSEP3R('X3',X3,R3)
      IF (EPS.EQ.999999.) THEN
        AUX=GAMMA*S2TH
        EPS=-AUX/(1.+AUX)
      ENDIF
      EPS2=EPS*EPS
      JEPS2=1.-EPS2
      SJEPS2=SQRT(JEPS2)
      SQJMEP=SQRT(1.-EPS)
      SQJPEP=SQRT(1.+EPS)
      IF (V0.EQ.999999.) THEN
        AUX=GAMMA*S2TH
        V0=SQRT(A44*(1+AUX))
      ENDIF
      IF (V0.EQ.0.) THEN
C       TCGREEN-03
        CALL ERROR('TCGREEN-03: Wrong value of V0')
C       V0 must not equal zero.
      ENDIF
      IF (VK.EQ.0.) THEN
C       TCGREEN-04
        CALL ERROR('TCGREEN-04: Wrong value of TCK')
C       TCK must not equal zero.
      ENDIF
      IF (FOUTS.NE.' ') THEN
        IF (A44.LE.0.) THEN
C         TCGREEN-05
          CALL ERROR('TCGREEN-05: Wrong value of A44')
C         A44 must be positive if TCSTRESS is not blank.
        ENDIF
      ENDIF
      DO 10, I2=1,14
        GREEN(I2)=0.
  10  CONTINUE
      TT=X3/V0
      GREEN(1)=TT
      GREEN(5)=X3
      NGREEN=14+18*NF
      IF (NGREEN.GT.MRAM) THEN
C       TCGREEN-02
        CALL ERROR('TCGREEN-02: Small array RAM.')
C       Parameter MRAM of the file ram.inc should be greater than
C       18 times number of frequencies plus 14.
C       File ram.inc.
      ENDIF
C
C
C     Exact solution:
      IF (FOUTE.NE.' '.OR.FOUTS.NE.' ') THEN
        IF (FOUTS.NE.' ') THEN
          OPEN(LU,FILE=FOUTS)
        END IF
        I2=15
        DO 100, I1=0,NF-1
          F=OF+I1*DF
          K0=2.*PI*F/V0
C         Eq 63:
          VK2=VK*VK
          K02=K0*K0
          K02VK2=K02-VK2
          AUX=1.-EPS2*((VK2/K02VK2)**2)
          REJM=K02VK2*SQRT(JEPS2)
          IF (AUX.GE.0.) THEN
            IMJM=0.
            REJM=REJM*(SJEPS2+SQRT(AUX))
          ELSE
            IMJM=REJM*SQRT(-AUX)
            REJM=REJM*SJEPS2
          ENDIF
          RECIT=EPS*VK*K02
          IF (IMJM.EQ.0.) THEN
            REF1=RECIT/REJM
            IMF1=0.
          ELSE
            AUX=REJM*REJM-IMJM*IMJM
            REF1=RECIT*REJM/AUX
            IMF1=-RECIT*IMJM/AUX
          ENDIF
C         Eq 52:
          REF12=REF1*REF1-IMF1*IMF1
          IMF12=2.*REF1*IMF1
          REPZ=K02/JEPS2+REF12
          IMPZ=IMF12
          REAUX=1.+REF12*JEPS2/VK2
          IMAUX=   IMF12*JEPS2/VK2
          AUX=REAUX*REAUX+IMAUX*IMAUX
          REDZ= REAUX/AUX
          IMDZ=-IMAUX/AUX
          REF02=REPZ*REDZ-IMPZ*IMDZ
          IMF02=REPZ*IMDZ+IMPZ*REDZ
          AF02=SQRT(REF02*REF02+IMF02*IMF02)
          IF (AF02.EQ.0.) THEN
            REF0=0.
            IMF0=0.
          ELSE
            CFIF02=REF02/AF02
            SFIF02=IMF02/AF02
            IF (CFIF02.GE.0.) THEN
              FIF02=ASIN(SFIF02)
            ELSEIF(SFIF02.GE.0.) THEN
              FIF02=ACOS(CFIF02)
            ELSE
              FIF02=-ACOS(CFIF02)
            ENDIF
            AF0=SQRT(AF02)
            FIF0=FIF02/2.
            REF0=AF0*COS(FIF0)
            IMF0=AF0*SIN(FIF0)
            IF (REF0.LT.0.) THEN
              REF0=-REF0
              IMF0=-IMF0
            ENDIF
            IF (IMF0.LT.0.) THEN
              IMF0=-IMF0
              IMF1=-IMF1
            ENDIF
          ENDIF
C         Eq 49:
          REAUX=-EPS+REF1*JEPS2/VK
          IMAUX=     IMF1*JEPS2/VK
          REF3=REF0*REAUX-IMF0*IMAUX
          IMF3=IMF0*REAUX+REF0*IMAUX
C         Eq 47:
          REF2=VK+EPS*REF1
          IMF2=   EPS*IMF1
C
C         Eq 32,33,39:
          CALL EQ39(VK,X3,TT,F,
     *              REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     *              REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22)
C         Remaining components of the 3x3 matrix:
          REU13=0.
          IMU13=0.
          REU23=0.
          IMU23=0.
          REU31=0.
          IMU31=0.
          REU32=0.
          IMU32=0.
          REU33=0.
          IMU33=0.
C
C         Rotation (multiplication E x U x ET):
          IF (LROT) THEN
            CALL TCGMAT(REE,IME,REU,IMU,REA,IMA)
            CALL TCGMAT(REA,IMA,REET,IMET,REU,IMU)
          ENDIF
C
          GREEN(I2   )=REU11*1000000.
          GREEN(I2+ 1)=IMU11*1000000.
          GREEN(I2+ 2)=REU21*1000000.
          GREEN(I2+ 3)=IMU21*1000000.
          GREEN(I2+ 4)=REU31*1000000.
          GREEN(I2+ 5)=IMU31*1000000.
          GREEN(I2+ 6)=REU12*1000000.
          GREEN(I2+ 7)=IMU12*1000000.
          GREEN(I2+ 8)=REU22*1000000.
          GREEN(I2+ 9)=IMU22*1000000.
          GREEN(I2+10)=REU32*1000000.
          GREEN(I2+11)=IMU32*1000000.
          GREEN(I2+12)=REU13*1000000.
          GREEN(I2+13)=IMU13*1000000.
          GREEN(I2+14)=REU23*1000000.
          GREEN(I2+15)=IMU23*1000000.
          GREEN(I2+16)=REU33*1000000.
          GREEN(I2+17)=IMU33*1000000.
          I2=I2+18
C
          IF (FOUTS.NE.' ') THEN
C           Writing the density normalized initial stress, divided by
C           "reference stress" = 2.*PI*F*SQRT(A44) = K0*V0*SQRT(A44):
            IF(F.EQ.0.) THEN
              S11R= 0.
              S11I= SQRT(1.-EPS**2)*V0/SQRT(A44)
              S12R= 0.
              S12I= 0.
              S22R= 0.
              S22I= SQRT(1.-EPS**2)*V0/SQRT(A44)
            ELSE
              S11R=-(1.+EPS)*(IMF0+IMF3)*V0/K0/SQRT(A44)
              S11I= (1.+EPS)*(REF0+REF3)*V0/K0/SQRT(A44)
              S12R=-(1.-EPS**2)*REF1    *V0/K0/SQRT(A44)
              S12I=-(1.-EPS**2)*IMF1    *V0/K0/SQRT(A44)
              S22R=-(1.-EPS)*(IMF0-IMF3)*V0/K0/SQRT(A44)
              S22I= (1.-EPS)*(REF0-REF3)*V0/K0/SQRT(A44)
            END IF
            WRITE(LU,'(7(F9.6,1X))') F,S11R,S11I,S12R,S12I,S22R,S22I
          ENDIF
 100    CONTINUE
C
        IF (FOUTS.NE.' ') THEN
          CLOSE(LU)
        ENDIF
C
        IF (FOUTE.NE.' ') THEN
C         Writing:
          OPEN(LU,FILE=FOUTE)
          WRITE(LU,'(A)') ' /'
          FORMAT(1:4)='(6A,'
          CALL FORM2(14,GREEN,GREEN,FORMAT(5:260))
          WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)),
     *     ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''',
     *                    (' ',GREEN(I),I=1,14)
          DO 110 I2=15,NGREEN-18,18
            FORMAT(1:4)='(1A,'
            CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
            WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17)
  110     CONTINUE
          I2=NGREEN-17
          FORMAT(1:4)='(1A,'
          CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
          WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /'
          WRITE(LU,'(A)') ' /'
          CLOSE(LU)
        ENDIF
      ENDIF
C
C
C     Coupling ray theory:
      IF (FOUTW.NE.' ') THEN
        OPEN(LU,FILE=FOUTW)
        WRITE(LU,'(A)') ' /'
        I2=15
        DO 200, I1=0,NF-1
          F=OF+I1*DF
          K0=2.*PI*F/V0
C         Eq 106:
          REF0=K0*(SQJPEP+SQJMEP)/2./SJEPS2
          IMF0=0.
          REF1=0.
          IMF1=0.
          REF2=VK
          IMF2=0.
          REF3=-K0*(SQJPEP-SQJMEP)/2./SJEPS2
          IMF3=0.
C         Eq 32,33,39:
          CALL EQ39(VK,X3,TT,F,
     *              REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     *              REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22)
C
          GREEN(I2   )=REU11*1000000.
          GREEN(I2+ 1)=IMU11*1000000.
          GREEN(I2+ 2)=REU21*1000000.
          GREEN(I2+ 3)=IMU21*1000000.
          GREEN(I2+ 4)=0.
          GREEN(I2+ 5)=0.
          GREEN(I2+ 6)=REU12*1000000.
          GREEN(I2+ 7)=IMU12*1000000.
          GREEN(I2+ 8)=REU22*1000000.
          GREEN(I2+ 9)=IMU22*1000000.
          GREEN(I2+10)=0.
          GREEN(I2+11)=0.
          GREEN(I2+12)=0.
          GREEN(I2+13)=0.
          GREEN(I2+14)=0.
          GREEN(I2+15)=0.
          GREEN(I2+16)=0.
          GREEN(I2+17)=0.
          I2=I2+18
 200    CONTINUE
C       Writing:
        FORMAT(1:4)='(6A,'
        CALL FORM2(14,GREEN,GREEN,FORMAT(5:260))
        WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)),
     *   ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''',
     *                  (' ',GREEN(I),I=1,14)
        DO 210 I2=15,NGREEN-18,18
          FORMAT(1:4)='(1A,'
          CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
          WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17)
  210   CONTINUE
        I2=NGREEN-17
        FORMAT(1:4)='(1A,'
        CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
        WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /'
        WRITE(LU,'(A)') ' /'
        CLOSE(LU)
      ENDIF
C
C
C     Quasi-isotropic approach:
      IF (FOUTQ.NE.' ') THEN
        OPEN(LU,FILE=FOUTQ)
        WRITE(LU,'(A)') ' /'
        I2=15
        CALL RSEP3R('VREF2',VR2,V0*V0)
        VR=SQRT(VR2)
        DO 300, I1=0,NF-1
          F=OF+I1*DF
          K0=2.*PI*F/V0
C         Eq 111:
          REF0=K0*V0/VR*(3./2.-1./2.*V0*V0/VR/VR)
          IMF0=0.
          REF1=0.
          IMF1=0.
          REF2=VK
          IMF2=0.
          REF3=-K0/2.*EPS*V0*V0*V0/VR/VR/VR
          IMF3=0.
C         Eq 32,33,39:
          CALL EQ39(VK,X3,TT,F,
     *              REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     *              REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22)
C
          GREEN(I2   )=REU11*1000000.
          GREEN(I2+ 1)=IMU11*1000000.
          GREEN(I2+ 2)=REU21*1000000.
          GREEN(I2+ 3)=IMU21*1000000.
          GREEN(I2+ 4)=0.
          GREEN(I2+ 5)=0.
          GREEN(I2+ 6)=REU12*1000000.
          GREEN(I2+ 7)=IMU12*1000000.
          GREEN(I2+ 8)=REU22*1000000.
          GREEN(I2+ 9)=IMU22*1000000.
          GREEN(I2+10)=0.
          GREEN(I2+11)=0.
          GREEN(I2+12)=0.
          GREEN(I2+13)=0.
          GREEN(I2+14)=0.
          GREEN(I2+15)=0.
          GREEN(I2+16)=0.
          GREEN(I2+17)=0.
          I2=I2+18
 300    CONTINUE
C       Writing:
        FORMAT(1:4)='(6A,'
        CALL FORM2(14,GREEN,GREEN,FORMAT(5:260))
        WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)),
     *   ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''',
     *                  (' ',GREEN(I),I=1,14)
        DO 310 I2=15,NGREEN-18,18
          FORMAT(1:4)='(1A,'
          CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
          WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17)
  310   CONTINUE
        I2=NGREEN-17
        FORMAT(1:4)='(1A,'
        CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
        WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /'
        WRITE(LU,'(A)') ' /'
        CLOSE(LU)
      ENDIF
C
C
C     Anisotropic ray theory:
      IF (FOUTA.NE.' ') THEN
        OPEN(LU,FILE=FOUTA)
        WRITE(LU,'(A)') ' /'
        I2=15
        DO 400, I1=0,NF-1
          F=OF+I1*DF
          K0=2.*PI*F/V0
C         Eq 116:
          REF0=K0*(SQJPEP+SQJMEP)/2./SJEPS2
          IMF0=0.
          REF1=0.
          IMF1=0.
          REF2=0.
          IMF2=0.
          REF3=-K0*(SQJPEP-SQJMEP)/2./SJEPS2
          IMF3=0.
C         Eq 32,33,39:
          CALL EQ39(VK,X3,TT,F,
     *              REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     *              REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22)
C
          GREEN(I2   )=REU11*1000000.
          GREEN(I2+ 1)=IMU11*1000000.
          GREEN(I2+ 2)=REU21*1000000.
          GREEN(I2+ 3)=IMU21*1000000.
          GREEN(I2+ 4)=0.
          GREEN(I2+ 5)=0.
          GREEN(I2+ 6)=REU12*1000000.
          GREEN(I2+ 7)=IMU12*1000000.
          GREEN(I2+ 8)=REU22*1000000.
          GREEN(I2+ 9)=IMU22*1000000.
          GREEN(I2+10)=0.
          GREEN(I2+11)=0.
          GREEN(I2+12)=0.
          GREEN(I2+13)=0.
          GREEN(I2+14)=0.
          GREEN(I2+15)=0.
          GREEN(I2+16)=0.
          GREEN(I2+17)=0.
          I2=I2+18
 400    CONTINUE
C       Writing:
        FORMAT(1:4)='(6A,'
        CALL FORM2(14,GREEN,GREEN,FORMAT(5:260))
        WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)),
     *   ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''',
     *                  (' ',GREEN(I),I=1,14)
        DO 410 I2=15,NGREEN-18,18
          FORMAT(1:4)='(1A,'
          CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
          WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17)
  410   CONTINUE
        I2=NGREEN-17
        FORMAT(1:4)='(1A,'
        CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
        WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /'
        WRITE(LU,'(A)') ' /'
        CLOSE(LU)
      ENDIF
C
C
C     Isotropic ray theory:
      IF (FOUTI.NE.' ') THEN
        OPEN(LU,FILE=FOUTI)
        WRITE(LU,'(A)') ' /'
        I2=15
        DO 500, I1=0,NF-1
          F=OF+I1*DF
          K0=2.*PI*F/V0
C         Eq 121:
          REF0=K0*(SQJPEP+SQJMEP)/2./SJEPS2
          IMF0=0.
          REF1=0.
          IMF1=0.
          REF2=VK
          IMF2=0.
          REF3=0.
          IMF3=0.
C         Eq 32,33,39:
          CALL EQ39(VK,X3,TT,F,
     *              REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     *              REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22)
C
          GREEN(I2   )=REU11*1000000.
          GREEN(I2+ 1)=IMU11*1000000.
          GREEN(I2+ 2)=REU21*1000000.
          GREEN(I2+ 3)=IMU21*1000000.
          GREEN(I2+ 4)=0.
          GREEN(I2+ 5)=0.
          GREEN(I2+ 6)=REU12*1000000.
          GREEN(I2+ 7)=IMU12*1000000.
          GREEN(I2+ 8)=REU22*1000000.
          GREEN(I2+ 9)=IMU22*1000000.
          GREEN(I2+10)=0.
          GREEN(I2+11)=0.
          GREEN(I2+12)=0.
          GREEN(I2+13)=0.
          GREEN(I2+14)=0.
          GREEN(I2+15)=0.
          GREEN(I2+16)=0.
          GREEN(I2+17)=0.
          I2=I2+18
 500    CONTINUE
C       Writing:
        FORMAT(1:4)='(6A,'
        CALL FORM2(14,GREEN,GREEN,FORMAT(5:260))
        WRITE(LU,FORMAT) '''',RECNAM(1:LENGTH(RECNAM)),
     *   ''' ''',SRCNAM(1:LENGTH(SRCNAM)),'''',
     *                  (' ',GREEN(I),I=1,14)
        DO 510 I2=15,NGREEN-18,18
          FORMAT(1:4)='(1A,'
          CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
          WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,I2+17)
  510   CONTINUE
        I2=NGREEN-17
        FORMAT(1:4)='(1A,'
        CALL FORM2(18,GREEN(I2),GREEN(I2),FORMAT(5:260))
        WRITE(LU,FORMAT) (' ',GREEN(I),I=I2,NGREEN),' /'
        WRITE(LU,'(A)') ' /'
        CLOSE(LU)
      ENDIF
C
      WRITE(*,'(A)') '+TCGREEN: Done.                  '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE EQ39(VK,X3,TT,F,
     *  REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     *  REU11,IMU11,REU21,IMU21,REU12,IMU12,REU22,IMU22)
C
      REAL VK,X3,TT,F,
     * REF0,IMF0,REF1,IMF1,REF2,IMF2,REF3,IMF3,
     * REU11,IMU11,REU12,IMU12,REU21,IMU21,REU22,IMU22
C
      REAL PI
      PARAMETER (PI=3.1415926)
      REAL AUX,REAUX,IMAUX,REF12,IMF12,REF22,IMF22,REF32,IMF32,
     * REFI,IMFI,AFI,FIFI,
     * REFI2,IMFI2,AFI2,FIFI2,
     * RVF11,IVF11,RVF12,IVF12,RVF21,IVF21,RVF22,IVF22,
     * REVF11,IMVF11,REVF12,IMVF12,REVF21,IMVF21,REVF22,IMVF22,
     * RA11,IA11,RA12,IA12,RA21,IA21,RA22,IA22,
     * RB11,IB11,RB12,IB12,RB21,IB21,RB22,IB22,
     * RC11,IC11,RC12,IC12,RC21,IC21,RC22,IC22,
     * RD11,ID11,RD12,ID12,RD21,ID21,RD22,ID22,
     * CFIFI2,SFIFI2
C
C      Eq 32:
       REF12=REF1*REF1-IMF1*IMF1
       IMF12=2.*REF1*IMF1
       REF22=REF2*REF2-IMF2*IMF2
       IMF22=2.*REF2*IMF2
       REF32=REF3*REF3-IMF3*IMF3
       IMF32=2.*REF3*IMF3
       REFI2=REF32+REF22-REF12
       IMFI2=IMF32+IMF22-IMF12
       AFI2=SQRT(REFI2*REFI2+IMFI2*IMFI2)
       IF (AFI2.EQ.0.) THEN
         REFI=0.
         IMFI=0.
       ELSE
         CFIFI2=REFI2/AFI2
         SFIFI2=IMFI2/AFI2
         IF (CFIFI2.GE.0.) THEN
           FIFI2=ASIN(SFIFI2)
         ELSEIF(SFIFI2.GE.0.) THEN
           FIFI2=ACOS(CFIFI2)
         ELSE
           FIFI2=-ACOS(CFIFI2)
         ENDIF
         AFI=SQRT(AFI2)
         FIFI=FIFI2/2.
         REFI=AFI*COS(FIFI)
         IMFI=AFI*SIN(FIFI)
       ENDIF
C      Eq 33:
       AUX=REFI*REFI+IMFI*IMFI
       IF (AUX.NE.0.) THEN
         REAUX= REFI/AUX
         IMAUX=-IMFI/AUX
         RVF11=REF3
         IVF11=IMF3
         RVF12=-IMF1+IMF2
         IVF12= REF1-REF2
         RVF21=-IMF1-IMF2
         IVF21= REF1+REF2
         RVF22=-REF3
         IVF22=-IMF3
         REVF11=RVF11*REAUX-IVF11*IMAUX
         IMVF11=RVF11*IMAUX+IVF11*REAUX
         REVF12=RVF12*REAUX-IVF12*IMAUX
         IMVF12=RVF12*IMAUX+IVF12*REAUX
         REVF21=RVF21*REAUX-IVF21*IMAUX
         IMVF21=RVF21*IMAUX+IVF21*REAUX
         REVF22=RVF22*REAUX-IVF22*IMAUX
         IMVF22=RVF22*IMAUX+IVF22*REAUX
       ELSE
         REVF11=0.
         IMVF11=0.
         REVF12=0.
         IMVF12=0.
         REVF21=0.
         IMVF21=0.
         REVF22=0.
         IMVF22=0.
       ENDIF
C      Eq 39:
       RA11=COS(VK*X3)
       RA12=-SIN(VK*X3)
       RA21=-RA12
       RA22=RA11
       AUX=SIN(REFI*X3)
       RB11=COS(REFI*X3)-IMVF11*AUX
       IB11=             REVF11*AUX
       RB12=            -IMVF12*AUX
       IB12=             REVF12*AUX
       RB21=            -IMVF21*AUX
       IB21=             REVF21*AUX
       RB22=COS(REFI*X3)-IMVF22*AUX
       IB22=             REVF22*AUX
       RC11=RA11*RB11+RA12*RB21
       IC11=RA11*IB11+RA12*IB21
       RC12=RA11*RB12+RA12*RB22
       IC12=RA11*IB12+RA12*IB22
       RC21=RA21*RB11+RA22*RB21
       IC21=RA21*IB11+RA22*IB21
       RC22=RA21*RB12+RA22*RB22
       IC22=RA21*IB12+RA22*IB22
       AUX=REF0*X3 - 2.*PI*F*TT
       REAUX=COS(AUX)
       IMAUX=SIN(AUX)
       RD11=REAUX*RC11-IMAUX*IC11
       ID11=REAUX*IC11+IMAUX*RC11
       RD12=REAUX*RC12-IMAUX*IC12
       ID12=REAUX*IC12+IMAUX*RC12
       RD21=REAUX*RC21-IMAUX*IC21
       ID21=REAUX*IC21+IMAUX*RC21
       RD22=REAUX*RC22-IMAUX*IC22
       ID22=REAUX*IC22+IMAUX*RC22
       AUX=SINH(IMFI*X3)
       RA11=COSH(IMFI*X3)-REVF11*AUX
       IA11=             -IMVF11*AUX
       RA12=             -REVF12*AUX
       IA12=             -IMVF12*AUX
       RA21=             -REVF21*AUX
       IA21=             -IMVF21*AUX
       RA22=COSH(IMFI*X3)-REVF22*AUX
       IA22=             -IMVF22*AUX
       IF (IMF0.NE.0.) THEN
         AUX=-IMF0*X3
         AUX=EXP(-IMF0*X3)
         RB11=AUX*RA11
         IB11=AUX*IA11
         RB12=AUX*RA12
         IB12=AUX*IA12
         RB21=AUX*RA21
         IB21=AUX*IA21
         RB22=AUX*RA22
         IB22=AUX*IA22
       ELSE
         RB11=RA11
         IB11=IA11
         RB12=RA12
         IB12=IA12
         RB21=RA21
         IB21=IA21
         RB22=RA22
         IB22=IA22
       ENDIF
       REU11=RD11*RB11-ID11*IB11 + RD12*RB21-ID12*IB21
       REU12=RD11*RB12-ID11*IB12 + RD12*RB22-ID12*IB22
       REU21=RD21*RB11-ID21*IB11 + RD22*RB21-ID22*IB21
       REU22=RD21*RB12-ID21*IB12 + RD22*RB22-ID22*IB22
       IMU11=RD11*IB11+ID11*RB11 + RD12*IB21+ID12*RB21
       IMU12=RD11*IB12+ID11*RB12 + RD12*IB22+ID12*RB22
       IMU21=RD21*IB11+ID21*RB11 + RD22*IB21+ID22*RB21
       IMU22=RD21*IB12+ID21*RB12 + RD22*IB22+ID22*RB22
       RETURN
       END
C
C=======================================================================
C
      SUBROUTINE TCGMAT(A,B,C,D,E,F)
C
C----------------------------------------------------------------------
C Subroutine to compute the product of two 3x3 complex matrices
C (E+iF)=(A+iB)*(C+iD)
      REAL A(3,3),B(3,3),C(3,3),D(3,3),E(3,3),F(3,3)
C Input:
C     A,B,C,D ... Real and imaginary parts of the two matrices.
C Output:
C     E,F ... Real and imaginary parts of the resulting matrix.
C
C.......................................................................
      E(1,1)=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1)+
     *       A(1,3)*C(3,1)-B(1,3)*D(3,1)
      E(2,1)=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1)+
     *       A(2,3)*C(3,1)-B(2,3)*D(3,1)
      E(3,1)=A(3,1)*C(1,1)-B(3,1)*D(1,1)+A(3,2)*C(2,1)-B(3,2)*D(2,1)+
     *       A(3,3)*C(3,1)-B(3,3)*D(3,1)
      E(1,2)=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2)+
     *       A(1,3)*C(3,2)-B(1,3)*D(3,2)
      E(2,2)=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2)+
     *       A(2,3)*C(3,2)-B(2,3)*D(3,2)
      E(3,2)=A(3,1)*C(1,2)-B(3,1)*D(1,2)+A(3,2)*C(2,2)-B(3,2)*D(2,2)+
     *       A(3,3)*C(3,2)-B(3,3)*D(3,2)
      E(1,3)=A(1,1)*C(1,3)-B(1,1)*D(1,3)+A(1,2)*C(2,3)-B(1,2)*D(2,3)+
     *       A(1,3)*C(3,3)-B(1,3)*D(3,3)
      E(2,3)=A(2,1)*C(1,3)-B(2,1)*D(1,3)+A(2,2)*C(2,3)-B(2,2)*D(2,3)+
     *       A(2,3)*C(3,3)-B(2,3)*D(3,3)
      E(3,3)=A(3,1)*C(1,3)-B(3,1)*D(1,3)+A(3,2)*C(2,3)-B(3,2)*D(2,3)+
     *       A(3,3)*C(3,3)-B(3,3)*D(3,3)
C
      F(1,1)=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1)+
     *       A(1,3)*D(3,1)+B(1,3)*C(3,1)
      F(2,1)=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1)+
     *       A(2,3)*D(3,1)+B(2,3)*C(3,1)
      F(3,1)=A(3,1)*D(1,1)+B(3,1)*C(1,1)+A(3,2)*D(2,1)+B(3,2)*C(2,1)+
     *       A(3,3)*D(3,1)+B(3,3)*C(3,1)
      F(1,2)=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2)+
     *       A(1,3)*D(3,2)+B(1,3)*C(3,2)
      F(2,2)=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2)+
     *       A(2,3)*D(3,2)+B(2,3)*C(3,2)
      F(3,2)=A(3,1)*D(1,2)+B(3,1)*C(1,2)+A(3,2)*D(2,2)+B(3,2)*C(2,2)+
     *       A(3,3)*D(3,2)+B(3,3)*C(3,2)
      F(1,3)=A(1,1)*D(1,3)+B(1,1)*C(1,3)+A(1,2)*D(2,3)+B(1,2)*C(2,3)+
     *       A(1,3)*D(3,3)+B(1,3)*C(3,3)
      F(2,3)=A(2,1)*D(1,3)+B(2,1)*C(1,3)+A(2,2)*D(2,3)+B(2,2)*C(2,3)+
     *       A(2,3)*D(3,3)+B(2,3)*C(3,3)
      F(3,3)=A(3,1)*D(1,3)+B(3,1)*C(1,3)+A(3,2)*D(2,3)+B(3,2)*C(2,3)+
     *       A(3,3)*D(3,3)+B(3,3)*C(3,3)
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TCGTRA(A,B,C,D)
C
C----------------------------------------------------------------------
C Subroutine to compute the transpose matrix to the 3x3 complex matrix
C (C+iD)=(A+iB)^T
      REAL A(3,3),B(3,3),C(3,3),D(3,3)
C Input:
C     A,B ... Real and imaginary parts of the input matrix.
C Output:
C     C,D ... Real and imaginary parts of the transposed matrix.
C
C.......................................................................
      C(1,1)=A(1,1)
      C(2,1)=A(1,2)
      C(3,1)=A(1,3)
      C(1,2)=A(2,1)
      C(2,2)=A(2,2)
      C(3,2)=A(2,3)
      C(1,3)=A(3,1)
      C(2,3)=A(3,2)
      C(3,3)=A(3,3)
C
      D(1,1)=B(1,1)
      D(2,1)=B(1,2)
      D(3,1)=B(1,3)
      D(1,2)=B(2,1)
      D(2,2)=B(2,2)
      D(3,2)=B(2,3)
      D(1,3)=B(3,1)
      D(2,3)=B(3,2)
      D(3,3)=B(3,3)
C
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C=======================================================================
C