The Reflectivity Method for Seismogram Synthesis at Surface and Buried Receivers
                      Version 1.0 (Jan 03, 2026)           

          Copyright 2026  K. Koketsu.  All rights reserved.            
                                                                        
        Before using this package, please read the following.      
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^       

                                                Aug 13, 2026  K. Koketsu
------------------------------------------------------------------------
<Notes>

+ When you publish results produced by this, you have to refer to:

     Koketsu. K. (2026). The reflectivity method for synthetic seismograms
	at a buried receiver, Seismol. Res. Lett., doi:10.1785/0220260001.

+ The author does not accept any responsibility for troubles resulting 
  from the use of this software.

+ You may distribute this to any person without modification.           
  If you would like to distribute this with modifications               
  (including partial revisions and file deletions), you have to         
  contact the author.

             author: Kazuki Koketsu (koketsu@eri.u-tokyo.ac.jp)
                     Earthquake Research Institute, University of Tokyo
                     Yayoi 1-1-1, Bunkyo-ku, Tokyo 113-0032, Japan
------------------------------------------------------------------------
<files>

 README  : this file
 rx.f    : Fortran 77 code for this version of the reflectivity method
 pplot.f : to plot obtained seismograms using a Calcomp-like library   
 rxdemo1 : sample input No.1 to `rx.f' [Fig.1B of Kohketsu (1985)]
 rxdemo2 : sample input No.2 to `rx.f' [Fig.2B of Kohketsu (1985)]
 rxdemo3 : sample input No.3 to `rx.f' [Fig.2b of Koketsu (2026); DU is inverted]
 rxdemo4 : sample input No.4 to `rx.f' [Fig.3b of Koketsu (2026)]
 ppdemo1.in: pplot input for rxdemo1 
 ppdemo2.in: pplot input for rxdemo2
 ppdemo3.in: pplot input for rxdemo3
 ppdemo4.in: pplot input for rxdemo4
------------------------------------------------------------------------
<installation>

1. Create an empty directory and unzip the package there.
2. gfortran -o rx rx.f
   gfortran -o pplot  pplot.f -lcalpak
   (Please replace 'calpak' to suit your environment.)
3. rx < rxdemo1; pplot < ppdemo1.in; and view fort.8 using an application for EPS files.
4. rx < rxdemo2; pplot < ppdemo2.in; and view fort.8 using an application for EPS files.
5. rx < rxdemo3; pplot < ppdemo3.in; and view fort.8 using an application for EPS files.
6. rx < rzdemo4; pplot < ppdemo4.in; and view fort.8 using an application for EPS files.
------------------------------------------------------------------------
<hints>

+ Computational overflows may happen, when a seismic wave of high
  frequency passes through a thick, slow layer.  If so and you cannot
  change the frequency and phase velocity ranges, please divide the
  layer into several having the same physical parameters.
 
+ Waveforms are sensitive to frequency window (FU - FO), phase velocity
  window (ALIM - BLIM), step size of FFT (1/(DT*2**INDEX)) and step size
  of wavenumber integration (DK).
  FU - FO and ALIM - BLIM should be wide enough to cover a seismogram
  spectrum.  Unless the FFT step is fine enough (INDEX is large enough),
  or DK is small enough, aliasing or ripples may happen.  However, too
  wide ranges lead to a long computation time.  Since their optimum
  values depend on medium and source parameters, and there is no way
  to find out them, please determine them by referring to the sample
  inputs.

+ WI avoids surface wave poles along the wavenumber integration path,
  and can depress aliasing, but too large WI may cause other numerical
  noises.  Please refer to the sample inputs again.

+ For checking whether resultant waveforms include numerical noises,
  please slightly modify the parameters, recalculate waveforms and
  compare them with the previous ones.
------------------------------------------------------------------------
<outputs>

'rx.f' writes outputs on 'fort.10'.
 
Line 1       WRITE(10,'(2G12.5)') DT,VRED
DT,VRED same as the input parameters
 
Line 2       WRITE(10,'(5I5)') NNC
NNC     number of waveforms
 
Line 3       WRITE(10,'(6F12.7)) R(KR),T(KR)
R       epicentral distance
T       absolute time of the initial point (s, TI+R/VRED)
 
Line 4       WRITE(10,'(5I5)') NDR
NDR     number of sampling points ((TL-TI)/DT + 1)
 
Line 5       WRITE(10,'(E15.7)') AMAX
 AMAX   maximum amplitude (cm, cm/s, cm/s**2)
 
Line 6-      WRITE(10,'(6F12.7)') (AMW(I),I=1,NDR)
AMW()   waveforms normalized by AMAX
 
(The lines 3 to 6- repeat NNC times.)
------------------------------------------------------------------------
<pplot.f>

This reads the above outputs and plots waveforms.  This was written using
a CALCOMP-like library called 'calpak', which was developed at RSES ANU. 
------------------------------------------------------------------------
