Description of the program F R E S A N Program FRESAN is designed for the computation of the frequ- ency response at a receiver or a system of receivers distributed regularly or irregularly along the Earth's surface, an internal interface or a vertical line simulating borehole. Short description of the program FRESAN Program FRESAN is a modification of program GBSYN from prog- ram package BEAM87 written by V.Cerveny. As input data, data ge- nerated in the program ANRAY and stored in the file LU2 are used. The file LU2 contains travel times, amplitudes and phase shifts for all specified receivers and all elementary waves. The file also contains the slowness vectors specifying corresponding rays at a source. If standard ray computations are required the data from LU2 are sufficient. If quasi-isotropic computations are required, data generated in the program WEAKAN and stored in the file LU6 are necessary. Program FRESAN generates file LU7 containing tables of frequency response for all considered receivers. The file LU7 can be used in program SYNFAN for computing synthetic seismograms for various high-frequency source time functions. To evaluate the frequency response, numerically very efficient fast frequency response (FFR) algorithm is used. The frequency response is computed completely only for one frequency for each elementary wave at a considered receiver. The successive evaluation of the frequency response for other frequencies requi- res only one complex-valued multiplication. In program FRESAN, several types of point sources can be con- sidered, namely explosive (implosive) source, single force sour- ce, moment type source including double couple source. The force should be specified in Newtons, the moments in Newtonmeters. If length and velocity units used in the ANRAY program are m and m/s then the resulting displacements due to the force or moment source are in 10**(-3)m. If the length and velocity units used in ANRAY program are km and km/s, the resulting displacements are in 10**(-12)m. The program FRESAN is organized so that it allows to include any frequency-dependent effects at the endpoints of rays through the routine FDEP to be supplied by the user. In the program FRESAN, the frequency-dependent amplitude- distance curves may be optionally plotted. Description of input and output files Main input data are read from the standard input by list-directed input (free format) and consist of a single line containing following data: 'LIN' 'LOU' 'LU2' 'LU7' 'LU6'/ Here: 'LIN' is the name of the input data file LIN. 'LOU' is the name of the output log file LOU. 'LU2' is the name of the input data file LU2, generated by program ANRAY. 'LU7' is the name of the output data file LU7. 'LU6' is the name of the input data file LU6, generated by program WEAKAN (read in only in the QI mode). / is a slash recomended in batch and script files to enable future extensions. Defaults: 'LIN'='fresan.dat' 'LOU'='fresan.out' 'LU2'='lu2.dat' 'LU7'='lu7.dat' 'LU6'='lu6.dat' Example of the main input data: 'fresan.ani' / 'fresan.qi' / Input data consist partially of the data generated by the program ANRAY, stored in the formatted form in the file LU2, of the data generated by the program WEAKAN (if QI approximation is used) stored in the formatted form in the file LU6, and partially of the additional input data read by list-directed input (free format), prepared by the user and stored in the file LIN. Output data describing the computations are stored in the file LOU. Output data containing computed frequency responses are stored in the file LU7. The program can generate a postscript file with a plot of frequency-dependent amplitude-distance curve. The data stored in LU2 The data are stored in LU2 in a formatted form. For details see the description of the content of the file LU2 in program ANRAY. 1) MTEXT FORMAT(A) 2) NDST,KSH,ILOC FORMAT(26I3) 3) XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS FORMAT(8F10.5) 4) IF(NDST.NE.1) DST(I),I=1,NDST FORMAT(8F10.5) IF(NDST.EQ.1) XREC,YREC,XPRF,XATAN FORMAT(8F10.5) 5) NCODE,II,T,CAX,CAY,CAZ,TAST FORMAT(2I3,F12.7,6E12.6,F10.6) 6) IF(NC.LT.0) CAX1,CAY1,CAZ1 FORMAT(6E12.6) 7) (P(K),K=1,3) FORMAT(3F10.5) 8) (POL(K),K=1,3) FORMAT(3F10.5) 9) IF(NC.LT.0) (POL1(K),K=1,3) FORMAT(3F10.5) The data 5-9 can be repeated in LU2 many times for different receivers and different elementary waves. In the program FRESAN, this number may not exceed 2000. Otherwise it would be necessary to change the dimensions in the program. Note that the maximum number of the receiver positions is 100. The data stored in LU6 The data are stored in LU6 in a formatted form. For details see the description of the content of the file LU2 in program WEAKAN. 1) FL,FD,NF FORMAT(2F10.5,I5) 2) (E(1,K,NTOT),K=1,3),(E(2,L,NTOT),L=1,3) FORMAT(16E15.5) 3) FF,W FORMAT(16E15.5) The data 3 repeat NF times for each considered ray. In the program FRESAN, the number NF may not exceed 1000. Otherwise it would be necessary to change the dimensions in the program. The data 2 and 3 repeat for every ray. The additional input data in the file LIN These data are specified by the user. The data control the computation of synthetic frequency responses. The places where the data from LU2 and LU6 are read in are denoted by **LU2/1, **LU2/2, **LU6/1, **LU6/2 etc. 1) Arbitrary alphanumeric text specifiyng the following data set. ITEXT Default 'FRESAN'. 2) Switches controlling the output of results into the output fi- le LOU and specifying the input and output files LU2 and LU7. IPRINT,JPRINT,PSTEXT IPRINT... controls the storage of results during the computation of the frequency response and the type of the name of the postscript file with a plot. IPRINT.GE.0... automatic generation of postscript file names: plot00.ps, plot01.ps, etc. IPRINT.LT.0... postscript file name specified by the user, see PSTEXT. IPRINT=0... input data only are stored. Default value. IABS(IPRINT)=1... only data sub 5 from LU2 are stored. IABS(IPRINT)=2... print of complex amplitude for each elementary wave and receiver. IABS(IPRINT)=3... auxiliary parameters in the frequency response computation are stored. JPRINT... controls the storage of frequency-dependent amplitude-distance curves. JPRINT=0... no amplitudes are stored. JPRINT=1... amplitudes are stored, see the additio- nal input data No.10 for details. PSTEXT... the name of the postscript file specified by the user; it is used only if IPR.LT.0. 3) Various switches controlling the computations. MCOMP,MABS,MSOUR,MSELEC,MEPIC,MFR,MPLOT MCOMP... controls the choice of the component of the displa- cement vector. MCOMP=0... vertical component (positive down). Default value. MCOMP=1... component along the x-axis. MCOMP=2... component along the y-axis. MABS... specifies the dissipation model. MABS=0... no dissipation, perfectly elastic medium. Default value. MABS=1... non-causal absorption. MABS=2... causal absorption, Muller's model, strict frequency-dependent approach. MABS=3... causal absorption, Muller's model, linea- rization approach. Note: In this version, only perfectly elastic medium with is considered, MABS=0. MSOUR... controls the type of the considered point source. MSOUR=0... unit isotropic radiation pattern. Default value. MSOUR=1... single force. MSOUR=2... double couple. MSOUR=3... explosive (implosive) source. MSOUR=4... moment tensor specified by its components. MSELEC... controls selection of elementary waves. MSELEC=0... all elementary waves stored at LU2 are used to construct synthetic frequency response. Default value. MSELEC=1... only selected waves are used to construct synthetic frequency response, see the additional input data No.7. MEPIC... controls selection of receivers. MEPIC=0... the frequency responses are constructed for all receiver positions. Default value. MEPIC=1... the frequency responses are constructed at selected receiver positions only, see additional data No.4. MFR... controls the application of user-supplied routine for frequency filtering FDEP. This filtering may differ from receiver to receiver. It may include effects of local geology, instrumental effects, etc. MFR=0... no filtering. Default value. MFR.GT.0.... frequency-dependdent filtering, the routine FDEP is called. Note: In this version, the routine FDEP is empty, MFR=0 by default. MPLOT... controls plotting of frequency-dependent amplitude- distance curves, see additional input data No.11. MPLOT=0... no plotting. Default value. MPLOT=1... plotting of non-reduced amplitude- distance curves. MPLOT=2... plotting of amplitude-distance curves with maximum amplitude normalized to unity . 4) Auxiliary quantities for absorption computations and the angle of rotation of the recording coordinate system. The auxiliary quantities for absorption have no meaning in this version. RFR,GM,QRED,FREQM,AROT RFR... reference frequency for absorption computations.The values of velocity and quality factor used in the model are supposed to be specified for the frequency RFR. For MABS=0,1, RFR has no meaning. Default value, RFR=1. GM... exponent in the power law of the Muller's dissipation model. In this version GM=0. This choice corresponds to the Kjartansson's model, and approximately also to the Futterman's model. QRED... factor by which the quantity TAST (global absorption factor) stored in LU2 is to be multiplied. QRED changes the quality factor distribution in the whole model. The actual quality factor used in the computations in FRESAN is thus Q/QRED where Q is the quality factor considered in ANRAY. Default value, QRED=1. FREQM... frequency at which Muller's dissipation model is li- nearized if MABS=3. The linearization approach gives good results if the signals used in the program SYNFAN (in which the file LU7 generated in FRESAN is processed) have prevailing frequencies close to FREQM, and their amplitude spectrum is narrow. The linearization approach increases considerably the speed of computations. AROT... an angle, in radians, by which the horizontal coordinate axes of the recording coordinate system are to be rotated. For AROT=0, the recording coordinate system coincides with the "model" coordinate system. For AROT.NE.0, the axes of the recording system are rotated anticlockwise in the horizontal plane by the angle AROT. Default AROT=0. 5) Specification of the frequency range considered in the compu- tations. In case of QI computations, i.e. when LU6.NE.0, the following data are read in only formally. The information about the considered frequency range is read in from the file LU6, line 1, see above. NFREQ,FL,FR,FD NFREQ... specifies the form of the sampling. NFREQ=0... sampling in frequency. Default value. NFREQ.GT.0... sampling in time. FL... lower limit of the considered frequency range. FR... upper limit of the considered frequency range. FD... sampling step. For NFREQ=0, frequency step. For NFREQ.GT.0, time step. The frequency step is then determined as: FD=1./(FD*FLOAT(NFREQ)). 6) Selection of receivers. Included only if MEPIC.NE.0. NEPIC,(IEP(I),I=1,NEPIC) NEPIC... number of excluded receiver positions. IEP(1),IEP(2),..,IEP(NEPIC)... sequential numbers of excluded receiver positions. 7) Selection of elementary waves. Included only if MSELEC.NE.0. NSELEC,(ISEL(I),I=1,NSELEC) NSELEC... number of excluded elementary waves. ISEL(1),ISEL(2),..,ISEL(NSELEC)... successive numbers of exc- luded elementary waves. **LU2/1 **LU2/2 **LU2/3 **LU2/4 8) Specification of source parameters. Included only if MSOUR.NE.0. The sources are specified with respect to the model coordinate system. IPAR(1),...,IPAR(4),PAR(1),...,PAR(6) MSOUR=1: Single force. IPAR(1)-IPAR(4)... free parameters. PAR(1)... force azimuth - the angle, in radians, made by the projection of the force vector into a horizontal plane, and the positive x-axis. Positive clockwise. PAR(2)... magnitude of the force (in N, i.e. kg*m/s**2); recommended value is of the order of 1000. PAR(3)... force declination - the angle, in radians, made by the force vector and a horizontal plane. PAR(4)-PAR(6)... free parameters. MSOUR=2: Double couple. IPAR(1)-IPAR(4)... free parameters. PAR(1)... dip angle - the angle, in radians, between the fault and a horizontal plane, in radians. PAR(2)... source moment (in Nm, i.e. kg*m**2/s**2). ; recommended value is of the order of 1000. PAR(3)... strike angle - the angle, in radians, between posi- tive x-axis and the intersection of the fault with a horizontal plane. PAR(4)... rake angle - the angle, in radians, between the di- rection of the displacement on the fault and the intersection of the fault with a horizontal plane. Positive anticlockwise. PAR(5)-PAR(6)... free parameters. Note: The formulae for the double couple follow the definition in Quantitative Seismology by Aki and Richards. MSOUR=3: Explosive (implosive) source. IPAR(1)=1 (-1)... for explosive (implosive) source. IPAR(2)-IPAR(4)... free parameters. PAR(1)... magnitude of the source (in Nm, i.e. kg*m**2/s**2). PAR(2)-PAR(6)... free parameters. MSOUR=4: Moment tensor (MIJ) specified by its components. IPAR(1)-IPAR(4)... free parameters. PAR(1)=M11, PAR(2)=M12, PAR(3)=M13, PAR(4)=M22, PAR(5)=M23, PAR(6)=M33. 9) Frequency filtering: routine FDEP. Included only if MFR.NE.0. In the present version, the routine FDEP is empty, use MFR=0. **LU6/1 **LU2/5 **LU2/6 **LU2/7 **LU2/8 **LU2/9 **LU6/2 **LU6/3 Reading from LU2 and LU6 is repeated unti the end of the files. 10) Specification of frequencies for which the tables of amplitude-distance curves are to be printed. Included only if JPRINT.NE.0. MTAB, (IFREQ(I),I=1,MTAB) MTAB... number of frequencies for which the tables of amplitude-distance curves are to be printed. IFREQ(I)... the number of the I-th frequency for which the table is to be printed. The frequency is given by the relation: FREQ=FL+FD*(IFREQ(I)-1). 11) Plotting of amplitude-distance curves for selected frequencies. Included only if MPLOT.NE.0. 11a) Arbitrary alphanumeric text to be shown in the plot TEXT 11b) Description of the x-axis of the plot. XMIN,XMAX,XLEN,DTICX,SC XMIN,XMAX... the minimum and maximum values along the range axis (in the user length units, e.g. km). Remember, the epicenter is automatically considered to be situated at x=0. XLEN... length of the x axis, in cm. DTICX... the distance between two neighbouring tics on the x-axis which are denoted by the corresponding coordinate values (in the user length units). DTICX.GT.0.0: Tic marks starting from XMIN and ap- pearing at the subsequent points XMIN+DTICX, XMIN+2.0*DTICX,... DTICX.LT.0.0: Tic marks start and continue to be plotted from the first integer multiple of ABS(DTICX) greater than XMIN. SC... scaling factor, controls the scales of tics and alphanumeric texts. For SC=1., the tics are 0.15cm long and coordinates and text describing the plot are 0.4 and 0.45 cm high, respectively. Default SC=1. 11c) Description of the amplitude axis of the plot. AMIN,AMAX,ALEN,DTICA AMIN,AMAX,ALEN,DTICA... the same meaning as in the case of the x-axis, but for the amplitude axis (decadic logarithms of amplitudes). 11d) Additional data for the plot. NLOG,NTICX,NTICA,NDX,NDA NLOG... specifies the variable plotted along the amplitude axis: NLOG=0... linear amplitude scale. Default value. NLOG.NE.0... logarithmic amplitude scale. NTICX... NTICX.NE.0: The number of marked intervals on the x-axis between two neighbouring tics with corresponding coordinate values. Default NTICX=1. NTICA... NTICA.NE.0: The same as NTICX, but for amplitude axis. Default NTICA=1 NDX,NDA...control the precision of numbers describing the coordinate axes in the plots. NDX corresponds to range axis, NDA to the amplitude axis. ND.GT.0... the number of digits to the right of the decimal point. ND=0... only integer portions of the numbers with decimal points. Default value. ND.LT.0... integers. 11e) Specification of frequencies, for which the amplitude-distance curves are to be plotted. MTAB, (IFREQ(I),I=1,MTAB) MTAB... number of frequencies for which the plots of amplitude-distance curves are to be plotted within one frame. MTAB=0 indicates the end of plotting. IFREQ(I)... the number of the I-th frequency for which the curve is to be plotted. The frequency is given by the relation: FREQ=FL+FD*(IFREQ(I)-1). The input data sub 11e may repeat several times. MTAB=0 indicates the end of plotting. Example of data LIN for anisotropic model Example of data LIN for isotropic model for QI computations Termination of computations. Input data sub 11 can be repeated several times to produce several plots of frequency-dependent amplitude-distance curves. To terminate the plotting and the whole computation in the prog- ram FRESAN insert two blank lines after the last line 11e. OUTPUT TABLES Output to the file LOU All the additional input data are automatically reproduced to the file LOU. Besides them the data LU2/1, LU2/2, LU2/3 and LU2/4 are automatically reproduced. The storage of additional data is controlled by parameters IPRINT and JPRINT, see the additional input data sub 2. For IPRINT=1, the following data are stored: the successive number of the considered elementary wave, the successive number of the receiver hit by the ray of the elementary wave, the arri- val time, the ray amplitude and phase of the considered component of the elementary wave in the "model" coordinate system, directivity of the source (ray amplitude at the source) and the global absorption factor (not considered in this version). This table is followed by two-column table consisting of coordinates of the receiver coordinates and corresponding maximum amplitudes. For IPRINT=2, the following data are stored: the real and imaginary parts of the ray amplitude of the considered component of the displacement vector including the source effects, for each considered receiver and elementary wave in the "model" coordinate system. This table is followed by the table of normalized values of frequency response for each receiver. For IPRINT=3, the following data are stored: the real and imaginary parts of the complex travel time CT (real part represents travel time, the imaginary part global absorption factor), two constant factors for recurent fast frequency response computation (cos(2*PI*FD*CT) and sin(2*PI*FD*CT)), real and imaginary parts of the amplitude corresponding to the lowest considered frequency FL. These six values are stored for each considered receiver and elementary wave. For JPRINT=1, the tables of frequency-dependent amplitude-distance curves are stored for selected frequencies, see additional input data sub 10. Each table is stored under the heading: FREQUENCY= (HZ), MAXIMUM AMPLITUDE= with value of frequency and corresponding maximum amplitude on the profile. Under the heading, reduced amplitudes for all consi- dered receivers are stored. The reduction is such that the maxi- mum amplitude has value 999. Output to the file LU7 File LU7 contains the frequency responses for specified sys- tem of receivers in the form required by the program SYNFAN, whe- re they are used for the computation of ray synthetic seismo- grams. The data are stored in the formatted form in the following order: 1) MTEXT FORMAT(A) Arbitrary alphanumeric text describing the computations; specified in the program ANRAY. 2) ITEXT FORMAT(A) Arbitrary alphanumeric text describing the computations; specified in the program FRESAN. 3) XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,FL,FD FORMAT(5F10.5,2E15.7) XSOUR,YSOUR,ZSOUR... coordinates of the source. TSOUR... initial time. RSTEP... the distance between receivers for a regularly dist- ributed receivers; the average distance between irregularly distributed receivers (in used length units). FL... lower limit of the considered frequency range. FD... frequency step. 4) NRDST,NF,NF1,MCOMP,ILOC FORMAT(26I3) NRDST.. number of selected receivers. NF... number of frequency samples within the frequency range [FL,FL+(NF-1)*FD]. NF1... number of frequency samples between zero frequency and FL. MCOMP... specifies the considered component of the displace- ment vector. MCOMP=0... vertical component. MCOMP=1... component along the x-axis. MCOMP=2... component along the y-axis. ILOC... controls the orientation and position of the profile, on which receivers are situated. ILOC=0: Receivers along a surface profile. ILOC=1: Receivers along a vertical profile. ILOC.GT.1: Receivers along the ILOC-th interface. 5) The tables of frequency response, successively for all selec- ted receivers. 5a) DST,A FORMAT(F10.3,E12.5) DST... profile coordinate of the receiver (range coordinate in the surface profile mode, depth coordinate in the VSP mode). A... maximum value of the real and imaginary parts of the frequency response. 5b) The table of normalized real and imaginary parts of the fre- quency response for the receiver DST, in a reduced form. For NF frequencies, starting from the frequency FD*FLOAT(NF1+1), with the step FD. Each line contains 12 integers III representing 6 normalized complex-valued values of frequency response (maximum value of III is 99999). III(J) FORMAT(12I6) Graphical output The frequency-dependent amplitude-distance curves for selec- ted frequencies may be plotted, see additional input data sub 11. The number of plotted figures is unlimited. One or several amplitude-distance curves corresponding to different frequencies may be plotted into one frame.