*** Modified in Nov., 1995 **** "grflt6.f" is modified from "grflt4.f" to include not only the "ramp function" but also the "smoothed ramp function" whose velocity is triangle. Recently, the smoothed ramp function is more popular than the classical ramp function in Seismology. You should use the new input file "grflt.in" in this directory to use "grflt6.f". (Y. Hisada) ************************************************ *** README File for grflt4 (September, 1995) *** ************************************************ ** The files in this directory are Fortran codes and m-files for Matlab to efficiently compute and plot seismograms and normal mode solusions in layered half-spaces. Those are improved versions of "grfault.f" files. They are much faster than "grfault.f". ** The all codes are open to academic users, although we feel that the codes are far from being complete and needs some modifications in near future. Please use them by your own risks. ** For comments or questions, please e-mail to "hisada@cc.kagakuin.ac.jp" or mail to: Yoshiaki Hisada Dept. of Architecture, Kogakuin Univ., Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan (FAX: 3-3340-0140) ** Please refer the following papers for the parts of the theories: 1) X. CHEN, Geophys. J. Intern., Vol.115, 391-409, 1993, for the normal mode solutions using the R/T matrix method. 2) Y.HISADA, Bull. Seismo. Soc. Am., V.84, Oct. 1994, (Part 1), for the formulation of the R/T matrix and an asymptotic method. 3) Y.HISADA, Bull. Seismo. Soc. Am., V.85, Aug. 1995, (Part 2), for the improved asymptotic method, which are used in the code here. 4) Y.HISADA, Bull. Seismo. Soc. Am., V.85, Oct. 1995, for comments and reply. ** In order to grasp an basic idea what you can get from those files, please follow the procedures listed below: 1) Compile and link "normal.f", and change the name of the exective file. If you are using Unix system, those procedures can be > f77 normal.f -o normal.exe 2) Run "normal.exe" to get the phase velocities of the surface waves. The input data for frequencies and the layered structutes are in "grflt.in", which are implicitly linked to the code. In addition, you have to enter some parameters to controll the accuracies and speed. As an example, we put those parameters in "a.dat". Therefore, the procedures will be > normal.exe < a.dat Then the code will start to compute the phase plus group velocities. For details, please check "normal.manual". 3) If you can use "Matlab" on the Unix system, you can plot the dispersion curves using the results obtained above and "m-files" in this directory. On the Matlab, we can just type >> phsl for the dispersion curves of the Love waves, and >> phsr for those of the Rayleigh waves. The graphs are stored in the postscript files "phsl.ps" and "phsr.ps", respectively. Those who cannot use Matlab, please check the files in this directory. 4) Compile and link "grflt4.f" and change the name of the exective file. > f77 grflt4.f -o grflt4.exe 5) Run the "grflt4.f" to get the waveforms from seismic sources. The input data for frequencies, structures, source parameters, the number of the wavenumber integration points, observation points and so on are in "grflt.in". In addition, this code uses the phase velocities obtained from "normal.f". The procedure is > grflt4.exe Please check "grflt4.manual" for the details. 6) Because the results obtained above are in the frequency domain, you have to FFT them to the time domain. For this purpose, you can use "grfftm.f". > f77 grfftm.f -o grfftm.exe Then run it, > grfft.exe < b.dat "b.dat" contains some parameters for the change of the signs of the results, differentiations, band-pass filtering, units and so on. Please check "grfftm.manual" for the details. The results obtained above are the NS, EW, and UD components of the high-cut filtered velocity waveforms with unit "cm". 7) If you use Matlab, you can plot the waveforms as follows; >> wavex for the X (NS) components, >> wavey for the Y (EW) components, and >> wavez for the Z (UD) components. The graphs are stored in the postscript files "wavex.ps", "wavey.ps" and "wavez.ps". If you can not use Matlab, Please check those files in this directory. ** If you want to check the details for the normal mode solutions and the integrands of the Green's functions at each frequency, seismic source and observation point, please follow the procedures below; 1) Change the parts of data in "grflt.in" as to be the frequency, the source, and the observation point which you want to check. As an example, we made an input file named "grflt.in.check". Therefore, please change its name as "grflt.in". The procedures can be > cp grflt.in grflt.in.original > cp grflt.in.check grflt.in In this example, we consider the circular frequency (omega) = 2.5, the first of the two seismic sources, and the last of the eleven observation points in the original "grflt.in". In addition, we changed the final velocity for the wavenumber integration from 500 (m/sec) to 2500 (m/sec) to enlarge the ranges including the poles. Please check the manuals for the details. 2) Run "normal.exe" to get the normal mode solutions. > normal.exe < a.dat 3) If you use Matlab, you can plot the secular functions (the characteristic functions) and the eigen vectors of the surface waves. To see the secular functions, we can type >> secl for the Love waves, and >> secr for the Rayleigh waves. In the figures, the circles are the roots of the secular functions, which represent the phase velocities at the frequency we chose. To see the eigen vectors, you can type >> vld >> vls for the displacement and stress vectors of the Love modes, and >> vrd >> vrs for those of the Rayleigh modes. In the figues, the horizontal dotted lines represent the locations of the layer boundaries. The all figures are stored in the postscript files "vld.ps", "vls.ps", "vrd.ps" and "vrs.ps". Those who cannot use Matlab, please check them. 3) Run "grflt4.exe" to output the integrands. > grflt4.exe Then, the code will ask you if you want to output them. "Output Integrands?: Enter No=0 or Yes=1:" Please enter "1" to output. 4) If you use Matlab, please type >> integ for the integrands without differentiations, >> integr for the integrands differentiated with respect to "r", and >> integz for the integrands differentiated with respect to "z". All figures are stored in the postscript files "integ.ps", "integr.ps", and "integz.ps", which you can find in this directory. ** The followings are brief explanation for the files. Plotting results using the following "m-files" of Matlab can be seen in postscript files in this directory. For more details, please check the manuals: "normal.manual", "grflt4.manual", and "grfftm.manual". 1) "normal.f" (Fortran code for the normal mode solutions) This code computes the phase and group velocities, secular functions, and eigen vectors of the Love and Rayleigh waves using the R/T matrix method. Regarding input data for the frequencies and structure, parts of "grflt.in" for "grflt4.f" are used. Please check "normal.manual" for more details. 2) "m-files" of Matlab for the secular functions and the eigen vectors of the normal mode solutions In order to use the following "m-files", number of frequancy (omega) must be one. Therefore, NOM and IOM in "grflt.in" must be same value. For example, when NOM=IOM=1 and DOM=1.0, the solusions are at omega=1. a. secl.m: Plot the secular functions of Love wave at one omega (Input data: sclc_dat, sclr_dat and scli_dat from "normal.f") secr.m: Plot the secular functions of Rayleigh wave at one omega (Input data: scrc_dat, scrr_dat and scri_dat from "normal.f") Note: Circles in the figure indicate the locations of phase velocities. b. vld.m: Plot the eigen displacement vectors of Love wave at one omega (Input data: vzn_dat and vl1_dat from "normal.f") vls.m: Plot the eigen stress vectors of Love wave at one omega (Input data: vzn_dat and vl2_dat from "normal.f") vrd.m: Plot the eigen displacement vectors of Rayleigh wave (Input data: vzn_dat vr1_dat and vr2_dat from "normal.f") vrs.m: Plot the eigen stress vectors of Rayleigh wave at one omega (Input data: vzn_dat vr3_dat and vr4_dat from "normal.f") 3) "m-files" for the phase and group velocities (dispersion curves) In order to use the following "m-files", number of frequancy (omega) must be more than one. Threrefore, NOM and IOM in "grflt.in" must not be same value. a. phsl.m: Plot the phase and group velocities of Love wave (Input data: om_dat, phsl_dat, and grpl_dat from "normal.f") phsr.m: Plot the phase and group velocities of Rayleigh wave (Input data: om_dat, phsr_dat, and grpr_dat from "normal.f") 4) "grflt4.f" (Fortran code for Green's function of layered half-space) This code computes displacements at multi observation points from seismic sources. The input data is "grflt.in". Before running this code, you must run "normal.f" to get the phase velocities, which correspont to "poles" in the integrands of the wavenumber integration. This code has the following features; a. We adopted the R/T matrix method as the propagator matrix. Therefore, it is numerically stable up to the infinite frequency. b. We used the asymptotic method by Hisada (1995). Therefore, the integrands of wavenumber integrations behave very well, even for cases that the depth of a source point is equal or nearly equal to that of observation points. c. Before starting the wavenumber integrations, we compute and store the R/T coefficients. This greatly reduces the computation time, when a lot of source points are needed. d. For large k*r (k: wavenumber, r: horizontal distance), Filon's integration formula is used instead of Simpson's. This also reduces the computation efforts greatly, because the integration points are needed to cover only the envelopes of oscillating integrands, but not the integrands themselves. e. Finer increments for the numerical integrations are used around the poles, which correspond to the phase velocities of the surface waves. You must run "normal.f" to get those data before running this code. Please check "grflt4.manual" for more details. 5) "m-files" of Matlab to plot the integrands In order to use the following "m-files", number of frequancy (omega) and observation point must be one. a. integ.m : Plot wavnumber vs integrands (Qx Wx Vx Qz Wz) (Input Data: qxr_dat qxi_dat wxr_dat wxi_dat vxr_dat vxi_dat qzr_dat qzi_dat wzr_dat wzi_dat from "grflt4.f") integr.m: Plot wavnumber vs integrands (Qx,r Wx,r Vx,r Qz,r Wz,r; the derivatives of Qx and so on with respect to r) (Input Data: qxrr_dat qxri_dat wxrr_dat wxri_dat vxrr_dat vxri_dat qzrr_dat qzri_dat wzrr_dat wzri_dat from "grflt4.f") integz.m: Plot wavnumber vs integrands (Qx,z Wx,z Vx,z Qz,z Wz,z; the derivatives of Qx and so on with respect to z) (Input Data: qxzr_dat qxzi_dat wxzr_dat wxzi_dat vxzr_dat vxzi_dat qzzr_dat qzzi_dat wzzr_dat wzzi_dat from "grflt4.f") b. qx.m : wavnumber vs integrand Qx (Input Data: qxr_dat qxi_dat) qxr.m: wavnumber vs integrand Qx,r (Input Data: qxrr_dat qxri_dat) qxz.m: wavnumber vs integrand Qx,z (Input Data: qxzr_dat qxzi_dat) qz.m : wavnumber vs integrand Qz (Input Data: qzr_dat qzi_dat) qzr.m: wavnumber vs integrand Qz,r (Input Data: qzrr_dat qzri_dat) qzz.m: wavnumber vs integrand Qz,z (Input Data: qzzr_dat qzzi_dat) Note: Please modify them to plot other components by yourself. c. qx1.m : wavnumber vs integrand Qx with locations of the poles (Input Data: qxr_dat qxi_dat, and phsl_dat and phsr_dat) qxr1.m: wavnumber vs integrand Qx,r with locations of poles (Input Data: qxrr_dat qxri_dat, and phsl_dat and phsr_dat) qxz1.m: wavnumber vs integrand Qx,z with locations of poles (Input Data: qxzr_dat qxzi_dat, and phsl_dat and phsr_dat) Note: Please modify them to plot other components by yourself. 6) Fortran codes for FFT All computations in "grflt4.f" are carried out in the frequency domain. Therefore, we need to Fourier transfer the results in order to obtain those in the time domain. a. grfft.f : FFT and filter three components at one observation point grfftm.f: FFT and filter three components at multi observation points Please check "grfftm.manual" for more details. 7) "m-files" of Matlab to plot waveforms a. wave.m : Plot three components at one observation point (Input-data: wave_dat from "grfft.f") b. orbit.m: Plot orbits and 3D trajectory of selected parts of waveforms (Input-data: wave_dat from "grfft.f") c. wavex.m : Plot X components at multi observation points (Input-data: x_dat, time_dat, and xyz_data from "grfftm.f") wavey.m : Plot Y components at multi observation points (Input-data: y_dat, time_dat, and xyz_data from "grfftm.f") wavez.m : Plot Z components at multi observation points (Input-data: z_dat, time_dat, and xyz_data from "grfftm.f")