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 initial angles of corresponding rays. Program FRESAN generates file LU8 containing tables of frequency response for all considered receivers. The file LU8 can be used in the 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, 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 the 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 data Input data consist partially of the data generated by the program ANRAY, stored in the formatted form in the file LU2, and partially of the additional input data prepared by the user and stored in the file LIN. Output data describing the computations are stored in the file LOU. Output data containing the computed frequency responses are stored in the file LU9. Specification of the files LIN, LOU, LU8 and LU9 can be made through the routine SERV, which is part of this program package. If the routine SERV is not used (put 'C' in front of 'CALL SERV' in the beginning of the main program), the files LIN and LOU are automatically specified LIN=5, LOU=6, the file with plot has number 7 and the number LU9 is read from the file LIN, see below. Program FRESAN can generate plots of frequency-dependent amplitude-distance curves. For this purpose the CALCOMP routines PLOTS, PLOT, SYMBOL and NUMBER are used. These routines are not included in the package, but must be linked with FRESAN. 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,VPS,VSS 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,CAX1, CAY1,CAZ1,DEC,AZIM,TAST FORMAT(2I3,F10.5,12E12.6,3F10.6) The data sub 5 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 additional input data These data are specified by the user. The data control the computation of synthetic frequency responses. The places where the data from LU2 are read in are denoted by **LU2/1, **LU2/2, etc. 1) Arbitrary alphanumeric text specifiyng the following data set. ITEXT FORMAT(A) 2) Switches controlling the output of results into the output fi- le LOU and specifying the input and output files LU2 and LU8. IPRINT,JPRINT,LU2,LU8 FORMAT(16I5) IPRINT... controls the storage of results during the computa- tion of the frequency response. IPRINT=0... input data only are stored. IPRINT=1... only data sub 5 from LU2 are stored. IPRINT=2... print of complex amplitude for each elementary wave and receiver. 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. LU2... the number of the file generated in ANRAY, in which the data for constructing synthetic frequency responses are stored. LU8... the number of the file in which tables of frequency response computed in FRESAN are stored. 3) Various switches controlling the computations. MCOMP,MABS,MSOUR,MSELEC,MEPIC,MFR,MPLOT FORMAT(16I5) MCOMP... controls the choice of the component of the displa- cement vector. MCOMP=0... vertical component. MCOMP=1... radial component. MCOMP=2... transverse component. MABS... specifies the dissipation model. MABS=0... no dissipation, perfectly elastic medium. 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. MSOUR=1... single force. MSOUR=2... double couple. MSOUR=3... explosive (implosive) source. Note: MSOUR.NE.0 can be used only for the source situated in an isotropic layer. MSELEC... controls selection of elementary waves. MSELEC=0... all elementary waves stored at LU2 are used to construct synthetic frequency response. 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. 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. 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.12. MPLOT=0... no plotting. 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 FORMAT(8F10.5) 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 LU8 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 by which the horizontal coordinate axes of the recording coordinate system are to be rotated. For AROT=0, the radial axis is parallel to the horizontal projection of the line connecting source with receivers. Transverse axis is perpendicular to it. For AROT.NE.0, the axes are rotated anticlockwise in the horizontal plane by the angle AROT. 5) Specification of the frequency range considered in the compu- tations. NFREQ,FL,FR,FD FORMAT(I5,6F10.5) NFREQ... specifies the form of the sampling. NFREQ=0... sampling in frequency. 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) FORMAT(16I5) 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) FORMAT(16I5) 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. IPAR(1),...,IPAR(4),PAR(1),...,PAR(6) FORMAT(4I5,6F10.5) 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). 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). 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. 9) Frequency filtering: routine FDEP. Included only if MFR.NE.0. In the present version, the routine FDEP is empty, use MFR=0. **LU2/5 Reading from LU2 is repeated upto the end of the file LU2. 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) FORMAT(16I5) 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 FORMAT(A) 11b) Description of the x-axis of the plot. XMIN,XMAX,XLEN,DTICX,SC FORMAT(8F10.5) 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.0, the tics are 0.15cm long and coordinates and text describing the plot are 0.4 and 0.45 cm high, respectively. 11c) Description of the amplitude axis of the plot. AMIN,AMAX,ALEN,DTICA FORMAT(8F10.5) 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 FORMAT(16I5) NLOG... specifies the variable plotted along the amplitude axis: NLOG=0... linear amplitude scale. 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. NTICA... NTICA.NE.0: The same as NTICX, but for amplitude axis. 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. ND.LT.0... integers. 11e) Specification of frequencies for which the tables of amplitude-distance curves are to be plotted. MTAB, (IFREQ(I),I=1,MTAB) FORMAT(16I5) MTAB... number of frequencies for which the tables 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 table 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. 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. In addition the heading of each table with frequency response stored in LU8/6 is stored. The storage of additional data is controlled by parameters IPRINT and JPRINT, see the additional input data sub 1. 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 amplitude and phase of the considered component of the elementary wave, declination and azimuth of the ray at the source and the global absorption factor (not considered in this version). This table is followed by a two-column table consisting of coordinates of the receivers and corresponding maximum ampli- tudes. 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 for each considered receiver and ele- mentary wave. This tablle 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 LU8 File LU3 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. 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 a (FL,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... radial component. MCOMP=2... transverse component. 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.