# History file 'tc.h' to analytically compute propagator matrices in
# model 'twisted crystal', and to compare the results with Green
# functions computed by CRT.
# Synthetic seismograms are generated and a figure of comparison of
# different ray-tracing methods is plotted in the second part of the
# history file.
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Executables required
# ~~~~~~~~~~~~~~~~~~~~
# This history file requires programs 'greentc.for', 'greenmul.for',
# 'tcgreen.for', 'tccomp.for' and 'tcfd.for' to be compiled.
# The programs are located in directory 'data/tc' and may be compiled
# by perl script 'ftc.pl'.
  chk.pl: "data/tc/" "greentc.for"
  chk.pl: "data/tc/" "greenmul.for"
  chk.pl: "data/tc/" "tcgreen.for"
  chk.pl: "data/tc/" "tccomp.for"
  chk.pl: "data/tc/" "tcfd.for"
  chk.pl: "data/tc/" "ftc.pl"
  chk.pl: "forms/"   "f.pl"
  ftc.pl:

# Input files required
  chk.pl: "data/tc/" "tc-cod.dat"
  chk.pl: "data/tc/" "tc-dcr.dat"
  chk.pl: "data/tc/" "force1.dat"
  chk.pl: "data/tc/" "tc-hil.dat"
  chk.pl: "data/tc/" "tc-hilc.dat"
  chk.pl: "data/tc/" "tc-mod.dat"
  chk.pl: "data/tc/" "tc-rec.dat"
  chk.pl: "data/tc/" "tc-rpa.dat"
  chk.pl: "data/tc/" "tc-src.dat"
  chk.pl: "data/tc/" "unit.dat"
  chk.pl: "crt/"     "writall.dat"
  chk.pl: "data/tc/" "tc-fdx.dat"
  chk.pl: "data/tc/" "tc-fdy.dat"
  chk.pl: "data/tc/" "tc-err1.pro"
  chk.pl: "data/tc/" "tc-err2.pro"
  chk.pl: "forms/"   "append.pl"
  chk.pl: "forms/"   "copy.pl"
  chk.pl: "forms/"   "echo.pl"

# Computation of propagator matrices and their comparison
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Model
  MODEL='tc-mod.dat'

# Source and receivers
  SRC='tc-src.dat'  REC='tc-rec.dat'

# Other input data files for crt.for
  DCRT='tc-dcr.dat'
  INIPAR=3
  CODE='tc-cod.dat'
  RPAR='tc-rpa.dat'
  WRIT='writall.dat'
  CRTLOG='log.out'

# Frequency step and interval
  OF=0. DF=0.0125 NF=209

# Calculating Green function by programs of package CRT
  GREENAMP=8171.6852       # conversion of point source to plane wave
  crt:                     # two-point ray tracing
# Coupling ray theory Green function
  GREEN='green.out' GREENMUL='greenmul.out' GREENTC='greencrt.out'
  green:                   # coupling ray theory Green function
  greenmul:                # multiplication by GREENAMP
  greentc:                 # conversion to propagator matrix U
# Result - propagator matrix written to the file 'greencrt.out'
# Isotropic ray theory Green function
  MODEL=                   # disabling anisotropy in 'green.for'
  GREEN='green.out' GREENMUL='greenmul.out' GREENTC='greeniso.out'
  green:                   # isotropic ray theory Green function
  greenmul:                # multiplication by GREENAMP
  greentc:                 # conversion to propagator matrix U
# Result - propagator matrix written to the file 'greeniso.out'

# Analytical solutions
  SIN2TH=.75  GAMMA=0.15  TCK=0.032  A44=6.  VREF2=6.9
  TCGREENE=  TCGREENW=  TCGREENQ=  TCGREENA=  TCGREENI=  TCSTRESS=' '
  tcgreen:                 # analytical Green functions
# Output in the form of Green functions written to files 'tcgreene.out',
#           'tcgreenw.out','tcgreenq.out','tcgreena.out','tcgreeni.out'.
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:                 # format conversion to matrix U exact
  GREEN='tcgreenw.out' GREENTC='greentcw.out'
  greentc:                 # format conversion to matrix U weak
  GREEN='tcgreenq.out' GREENTC='greentcq.out'
  greentc:                 # format conversion to matrix U qi
  GREEN='tcgreena.out' GREENTC='greentca.out'
  greentc:                 # format conversion to matrix U aniso
  GREEN='tcgreeni.out' GREENTC='greentci.out'
  greentc:                 # format conversion to matrix U iso
# Result - propagator matrices written to the files 'greentce.out',
#           'greentcw.out','greentcq.out','greentca.out','greentci.out'.

# Comparison inv(U1)*(U2-U3) of plane-wave propagator matrices
# Analytic isotropic ray method - isotropic ray method computed by CRT
  GREENTC1='unit.dat'     GREENTC2='greentci.out' GREENTC='greeniso.out'
  TCCOMP='erriso.out'      # output file
  tccomp:
# Analytic coupling ray theory - coupling ray theory computed by CRT
  GREENTC1='unit.dat'     GREENTC2='greentcw.out' GREENTC='greencrt.out'
  TCCOMP='errcrt.out'      # output file
  tccomp:
# Exact analytic solution - analytic ray methods
  GREENTC1='unit.dat'     GREENTC2='greentce.out'
  GREENTC='greentcw.out'  TCCOMP='errw.out'
  tccomp:
  GREENTC='greentcq.out'  TCCOMP='errq.out'
  tccomp:
  GREENTC='greentca.out'  TCCOMP='erra.out'
  tccomp:
  GREENTC='greentci.out'  TCCOMP='erri.out'
  tccomp:


# Synthetic seismograms
# ~~~~~~~~~~~~~~~~~~~~~
# Note  For X1 component, change BoundingBox from B1 B2   B3 B4
#                                              to B1 B2+4 B3 B4-5

# Anisotropic travel times for plotting
  SPHILI='tc-hil.dat'  SPHIWI=7.0

# Data describing the source time function
  SOURCE='force1.dat'
  KSIG=1 # (Gabor signal)
  SIGT=0.  SIGF=1.30  SIGW=4.  SIGPH=0.  SIGA=1.

# Data describing the filtration of the source time function
  FMIN=0.00  FLOW=0.13  FHIGH=2.47  FMAX=2.60  # (cosine filter)

# Time step and time interval for the Fast Fourier Transform
  DT=0.01953125  NFFT=4096

# Optional plotting of the source time function in program SS
  SIGPLOT='tc-sig.ps'  SIGGSE='tc-sig.gse'  MPTS=200
  ss:                        # source time function
  SIGPLOT=' '  SIGGSE=' '    # avoiding repeated plots of the s.t.f

# Computation of synthetic seismograms
  SIGPLOT=' '  SIGGSE=' '
# Exact Green function
  GREEN='tcgreene.out' SS='tc-sse.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
# Analytical coupling ray theory Green function
  GREEN='tcgreenw.out' SS='tc-ssw.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
# Analytical quasi-isotropic Green function
  GREEN='tcgreenq.out' SS='tc-ssq.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
# Analytical anisotropic Green function
  GREEN='tcgreena.out' SS='tc-ssa.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
# Analytical isotropic Green function
  GREEN='tcgreeni.out' SS='tc-ssi.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms

# Plotting seismograms
  SS='tc-ssi.gse'   KOLOR=2    # isotropic       red
  SS1='tc-ssa.gse'  KOLOR1=4   # anisotropic     blue
  SS2='tc-ssq.gse'  KOLOR2=19  # quasi-isotropic deep yellow
  SS3='tc-ssw.gse'  KOLOR3=3   # coupling        green
  SS4='tc-sse.gse'  KOLOR4=1   # exact           black
  SPTMIN=35.  SPTMAX=42.  SPTLEN=14.0  SPTDIV=7   SPTSUB=10
  KODESP=0    SPXLEN=7.0  NORMSP=1     SPAMP=3.5  SPCHRH=0.3
  SP1='tc-ss1.ps'  SP2='         '  SP3=' '
  sp:                        # seismogram plotting
  SPTDIV=-7
  SP1='         '  SP2='tc-ss2.ps'  SP3=' '
  sp:                        # seismogram plotting


# Synthetic seismograms at coupling frequency
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Anisotropic travel times for plotting
  SPHILI='tc-hilc.dat'  SPHIWI=7.0

# Data describing the coupling-frequency source time function
  SIGT=0.  SIGF=0.26  SIGW=4.  SIGPH=0.  SIGA=1.
  FMIN=0.00  FLOW=0.026  FHIGH=2.47  FMAX=2.60  # (cosine filter)

# Time step and time interval for the Fast Fourier Transform
  DT=0.0390625  NFFT=2048

# Computation of synthetic seismograms
  GREEN='tcgreene.out' SS='tc-ssec.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
  GREEN='tcgreenw.out' SS='tc-sswc.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
  GREEN='tcgreenq.out' SS='tc-ssqc.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
  GREEN='tcgreena.out' SS='tc-ssac.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms
  GREEN='tcgreeni.out' SS='tc-ssic.gse'
  greenss:                   # response function
  ss:                        # synthetic seismograms

# Plotting seismograms
  SS='tc-ssic.gse'   KOLOR=2
  SS1='tc-ssac.gse'  KOLOR1=4
  SS2='tc-ssqc.gse'  KOLOR2=19
  SS3='tc-sswc.gse'  KOLOR3=3
  SS4='tc-ssec.gse'  KOLOR4=1
  SPTMIN=32.  SPTMAX=46.  SPTLEN=14.0  SPTDIV=14  SPTSUB=10
  KODESP=0    SPXLEN=7.0  NORMSP=1     SPAMP=3.5  SPCHRH=0.3
  SP1='tc-ss1c.ps'  SP2='          '  SP3=' '
  sp:                        # seismogram plotting
  SPTDIV=-14
  SP1='          '  SP2='tc-ss2c.ps'  SP3=' '
  sp:                        # seismogram plotting


# Figure comparing the errors of four ray methods in frequency domain
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Comparisons inv(U1)*(U2-U3) are plotted, where U1 is unit matrix,
# U2 is the exact propagator matrix, and U3 is the analytic propagator
# matrix corresponding to the selected method

  GREENTC1='unit.dat'  GREENTC2='greentce.out'

  OF=0.001  DF=0.00025   NF=397
  TCSTRESS='stress1.out'
  tcgreen:                 # analytic Green functions
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:
# Comparisons
  GREEN='tcgreenw.out' GREENTC='greentcw.out'  TCCOMP='errw1.out'
  greentc:
  tccomp:
  GREEN='tcgreenq.out' GREENTC='greentcq.out'  TCCOMP='errq1.out'
  greentc:
  tccomp:
  GREEN='tcgreena.out' GREENTC='greentca.out'  TCCOMP='erra1.out'
  greentc:
  tccomp:
  GREEN='tcgreeni.out' GREENTC='greentci.out'  TCCOMP='erri1.out'
  greentc:
  tccomp:
# Storing the exact analytical propagator matrix for comparison with FD
  copy.pl:   "greentce.out" "tc-anal.out"

  OF=0.1  DF=0.0025  NF=361
  TCSTRESS='stress2.out'
  tcgreen:                 # analytic Green functions
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:
# Comparisons
  GREEN='tcgreenw.out' GREENTC='greentcw.out'  TCCOMP='errw2.out'
  greentc:
  tccomp:
  GREEN='tcgreenq.out' GREENTC='greentcq.out'  TCCOMP='errq2.out'
  greentc:
  tccomp:
  GREEN='tcgreena.out' GREENTC='greentca.out'  TCCOMP='erra2.out'
  greentc:
  tccomp:
  GREEN='tcgreeni.out' GREENTC='greentci.out'  TCCOMP='erri2.out'
  greentc:
  tccomp:
# Storing the exact analytical propagator matrix for comparison with FD
  append.pl: "tc-anal.out" "greentce.out"

  OF=1.0  DF=0.025  NF=361
  TCSTRESS='stress3.out'
  tcgreen:                 # analytic Green functions
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:
# Comparisons
  GREEN='tcgreenw.out' GREENTC='greentcw.out'  TCCOMP='errw3.out'
  greentc:
  tccomp:
  GREEN='tcgreenq.out' GREENTC='greentcq.out'  TCCOMP='errq3.out'
  greentc:
  tccomp:
  GREEN='tcgreena.out' GREENTC='greentca.out'  TCCOMP='erra3.out'
  greentc:
  tccomp:
  GREEN='tcgreeni.out' GREENTC='greentci.out'  TCCOMP='erri3.out'
  greentc:
  tccomp:

# Changing reference velocity
  VREF2=6.675
  TCSTRESS=' '
  OF=0.001  DF=0.00025   NF=397  TCCOMP='err01.out'
  tcgreen:
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:
  GREEN='tcgreenq.out' GREENTC='greentc0.out'
  greentc:
  tccomp:
  OF=0.1    DF=0.0025    NF=361  TCCOMP='err02.out'
  tcgreen:
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:
  GREEN='tcgreenq.out' GREENTC='greentc0.out'
  greentc:
  tccomp:
  OF=1.0    DF=0.025     NF=361  TCCOMP='err03.out'
  tcgreen:
  GREEN='tcgreene.out' GREENTC='greentce.out'
  greentc:
  GREEN='tcgreenq.out' GREENTC='greentc0.out'
  greentc:
  tccomp:

# Writing the PostScript figure
  copy.pl: tc-err1.pro                     tc-err.ps
  echo.pl: "stroke 1 0 0 setrgbcolor"   ">>tc-err.ps"
  append.pl: tc-err.ps erri1.out
  append.pl: tc-err.ps erri2.out
  append.pl: tc-err.ps erri3.out
  echo.pl: "stroke 0 0 1 setrgbcolor"   ">>tc-err.ps"
  append.pl: tc-err.ps erra1.out
  append.pl: tc-err.ps erra2.out
  append.pl: tc-err.ps erra3.out
  echo.pl: "stroke 1 0.8 0 setrgbcolor" ">>tc-err.ps"
  append.pl: tc-err.ps errq1.out
  append.pl: tc-err.ps errq2.out
  append.pl: tc-err.ps errq3.out
  echo.pl: "stroke 1 0.8 0 setrgbcolor" ">>tc-err.ps"
  append.pl: tc-err.ps err01.out
  append.pl: tc-err.ps err02.out
  append.pl: tc-err.ps err03.out
  echo.pl: "stroke 0 1 0 setrgbcolor"   ">>tc-err.ps"
  append.pl: tc-err.ps errw1.out
  append.pl: tc-err.ps errw2.out
  append.pl: tc-err.ps errw3.out
  append.pl: tc-err.ps tc-err2.pro

# Comparison of the analytical propagator matrix with FD by V.Vavrycuk
  copy.pl:   "stress1.out" "stress.out"
  append.pl: "stress.out" "stress2.out"
  append.pl: "stress.out" "stress3.out"
  TCFDX='tc-fdx.dat'  TCFDY='tc-fdy.dat'  TCSTRESS='stress.out'
  TCFD='tc-fd.out'  GREENTC='tc-anal.out'  TCCOMP='tc-fderr.out'
  tcfd: