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.

***
    The components of the displacement vector into the x and y
coordinates of the model coordinate system are formally called
"radial" and "transverse" components, respectively. They are
truly radial and transverse if the profile of receivers lies
along the x-axis. By use of the angle AROT, see below, horizontal
components can be transformed into a rotated coordinate system.
***

                                                     
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' 'LU8'/
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.
    'LU8' is the name of the output data file LU8.
    / is a slash recomended in batch and script files to enable future
        extensions.
Defaults:
    'LIN'='fresan.dat'
    'LOU'='fresan.out'
    'LU2'='lu2.in'
    'LU8'='lu8.out'
Example of the main input data:
    'fresan.sch' /

    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  LU8.
The program  generates a postscript file with the desired plot.

                                                      
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 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 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                            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.

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.
              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.
              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   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.

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.

   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.

**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.

Example of data LIN

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 LU8 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.