C Auxiliary main program RTCOEF (for testing C of routine COEF52) C ******************************* C For testing purposes, a brief main program RTCOEF is included C to control the computation of R/T coefficients using the routine C COEF52. Main program first calls COEF52 to read the data for the C model, and, after this, it reads two records specifying the C required computations. The computations are performed either in an C angle loop, or in a frequency loop. For any required angle of C incidence and frequency, routine COEF52 is called to compute the C relevant R/T coefficient. C C The complete set of input data is as follows: C ******************************************** C 1) One record: NZ C NZ ... number of layers, including both halfspaces C 2) NZ records: VP,VS,RHO,QP,QS,D C Parameters of the layers; one record for one layer. The C first record corresponds to the halfspace with the incident C wave (first halfspace), the last record to the NZ-th layer, C that is to the second halfspace. C 3) One record: FREF C FREF ... reference frequency (in Hz). C 4) One record: NC,NH,NQ,NF,NA C NC ... type of R/T coefficient C NH ... NH=0...the second halfspace is solid or liquid C NH=1...the second halfspace is a vacuum C NQ ... NQ=0... model is non-dissipative C NQ=1... model is dissipative C NF ... number of frequencies in the frequency loop. C For NF=1, the frequency loop is not considered C NA ... number of angles of incidence in the angle loop. C For NA=1, the angle loop is not considered C 5) One record: FMIN,DF,AMIN,DA,GAMMA C FMIN,DF ... specify frequency loop, in Hz C AMIN,DA ... specify angle loop, in degrees C GAMMA ... attenuation angle, in degrees C The records 4 and 5 can be repeated any number of times. The C computation finishes if NC=0 is used in record 4. C If we wish to compute the R/T coefficients in an angle loop C (for one given frequency FMIN), we use NF=1. If we wish to C compute the R/T coefficients in a frequency loop (for one C given angle of incidence AMIN), we use NA=1. Use always NF=1 C or NA=1. C The input data are stored in the file input.dat. C The results of computations are stored in the file C output.dat. For convenience, all input data are first stored C in output.dat (even those for model). Then, after an empty C line, the results of computations are stored. Each line diplays: C ANGLE,FREQ,RMOD,RPHASE,RMOD0,RPH0. C C Test examples C ************* C Four test examples are included. See the description in the file C 'rtcoef.htm'. C C The test data are located in subdirectory C rtcoef C of package DATA. C The test examples may be executed by command C perl go.pl rtcoef.h C running the demonstration history file C rtcoef.h. C C References C ********** C Brokesova,J. (2000). Reflection/transmission coefficients at a C plane interface in dissipative and non-dissipative media: C A comparison. J.Comput.Acoustics, 9,623 -641. C Brokesova,J., and Cerveny,V. (1998). Inhomogeneous plane waves C in dissipative, isotropic and anisotropic media. Reflection/ C transmission coefficients. In Seismic waves in complex 3-D C structures, Report No. 7, pp. 57 - 146. Department of C Geophysics, Charles University, Prague. C Cerveny,V. (1989). Synthetic body wave seismograms for laterally C varying media containing thin transition layers. Geophys. J. C Int., 99, 331-349, C Cerveny,V. (2001). Seismic ray theory. Cambridge Univ. Press, C Cambridge. C Muller,G. (1985). The reflectivity method. A tutorial. J.Geophys., C 58, 153-174. C C C Consortium Project "Seismic Waves in Complex 3-D Structures" C Praha, April 2003 C J.Brokesova, V.Cerveny C ************************************************************************ c c program rtcoef c ************** c c Auxiliary program rtcoef is designed for computations c of R/T coefficients of inhomogeneous P, SV and SH plane wave c at stack of layers between two isotropic anelastic halfspaces. c It uses the routine coef52. c open(7,file='input.dat') open(9,file='output.dat') c c reading the model call coef52(0,0,0,0,7,9,0.,0.,0.,0.,0.,0.,0.) c c reading the data for coef52 1 read(7,100)nc,nh,nq,nf,na write(9,100)nc,nh,nq,nf,na if(nc.eq.0.or.nc.eq.11.or.nc.eq.12.or.nc.gt.14)go to 10 read(7,101)fmin,df,amin,da,gamma write(9,101)fmin,df,amin,da,gamma write(9,102) c c angle loop if(na.eq.1)go to 5 freq=fmin do 2 i=1,na angle=amin+(i-1)*da call coef52(nc,nh,nq,1,7,9,angle,gamma,freq,rmod,rphase, 1rmod0,rph0) write(9,101)angle,freq,rmod,rphase,rmod0,rph0 2 continue go to 10 c c frequency loop 5 angle=amin do 6 i=1,nf freq=fmin+(i-1)*df call coef52(nc,nh,nq,1,7,9,angle,gamma,freq,rmod,rphase, 1rmod0,rph0) write(9,101)angle,freq,rmod,rphase,rmod0,rph0 6 continue c 7 go to 1 100 format(6i5) 101 format(8f10.3) 102 format (/) 10 continue stop end c ************************************************************************ c INCLUDE 'coef52.for' c coef52.for c ************************************************************************ c