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.ant' /
    'fresan.anv' /
    'fresan.qit' /
    'fresan.qiv' /

    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,
see also the input SEP parameters
controlling the form of the output PostScript file.

                                                      
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 of the "model"
              coordinate system.
              MCOMP=2... component along the y-axis of the "model"
              coordinate system.
    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  horizontal  axes of  the  recording
              system are  rotated  clockwise around the z-axis 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 anisotropic model
Example of data LIN for isotropic model for QI computations
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.