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 Reference:
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 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 Data describing the twisted-crystal model:
C     SIN2TH=real ... Second power of sinus theta
C              (equation 15 of the paper).
C              Default: SIN2TH=0.
C     GAMMA=real ... Gamma of equation 15.
C              Default: GAMMA=0.
C     TCK=real ... Uppercase K of equation 9.
C              Default: TCK=0.
C     EPSILON=real ... Epsilon of equation 9.
C              Default: EPSILON=-(GAMMA*SIN2TH)/(1.+GAMMA*SIN2TH)
C     A44=real ...  a44 of equation 152.
C              Default: A44=0.
C     V0=real ... v0 of equation 8.
C              Default: V0=SQRT(A44*(1.+GAMMA*SIN2TH))
C     VREF2=real ... second power of reference velocity
C              (see equation 111).
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     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
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,EQ39
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C             File sep.for.
C     FORM2 ... File forms.for.
C     EQ39 ... 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,FOUT,FILREC,FILSRC,RECNAM,SRCNAM
      CHARACTER*260 FORMAT
      INTEGER LU
      PARAMETER (LU=1)
      REAL PI
      PARAMETER (PI=3.1415926)

C     Auxiliary storage locations:
      INTEGER I1,I2,I,NF,NGREEN
      REAL EPS,S2TH,GAMMA,OF,DF,F,VK,K0,A44,V0,X3,K02,VK2,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,REU12,IMU12,REU21,IMU21,REU22,IMU22,
     * CFIF02,SFIF02,TT,VR,VR2,R1,R2,R3
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     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.)
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
      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:
      CALL RSEP3T('TCGREENE',FOUT,'tcgreene.out')
      OPEN(LU,FILE=FOUT)
      WRITE(LU,'(A)') ' /'
      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
        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
 100  CONTINUE
C
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 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)
C
C
C     Coupling ray theory:
      CALL RSEP3T('TCGREENW',FOUT,'tcgreenw.out')
      OPEN(LU,FILE=FOUT)
      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)
C
C
C     Quasi-isotropic approach:
      CALL RSEP3T('TCGREENQ',FOUT,'tcgreenq.out')
      OPEN(LU,FILE=FOUT)
      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)
C
C
C     Anisotropic ray theory:
      CALL RSEP3T('TCGREENA',FOUT,'tcgreena.out')
      OPEN(LU,FILE=FOUT)
      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)
C
C
C     Isotropic ray theory:
      CALL RSEP3T('TCGREENI',FOUT,'tcgreeni.out')
      OPEN(LU,FILE=FOUT)
      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)
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
      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