NAME
     hspec91p - evaluate F(f,p) for plane  wave  (tau-p)  elastic
     modeling

VERSION
     (2.0) - New Formulation
     (3.0) - Added Kjartansson (1979) constant Q operator

AUTHOR
     R. B. Herrmann, Saint Louis University, 1991

LANGUAGE
     FORTRAN 77

LIBRARY REQUIRED
     None

     Computations are performed in complex double precision.  The
     following non-standard double precision complex routines are
     required: zsqrt(), zabs(), dcmplx(), and zlog().

SYNOPSIS
     hspec91p
      [-O  name4 ]
      [-tmax  tmax ]
      [-dt  delt ]
      [-pmax  pmax ]
      [-pmin  pmin ]
      [-pnum  pnum ]
      [-F  fref ]
      [-K]
      [-SD]
      [-SU]
      [-SPUP]
      [-SPDN]
      [-SSUP]
      [-SSDN]
      [-RD]
      [-RU]
      [-RPUP]
      [-RPDN]
      [-RSUP]
      [-RSDN]
      [-?]

DESCRIPTION
     The purpose of this program is to evaluate  the  medium  p-t
     response  for  sixteen basic Green's functions as a function
     of distance and frequency. The Green's functions are  for  a
     point  vertical  downward force, a point horizontal force, a
     point explosion, and double couple earthquake sources.

     The meaning the the command line arguments are as follow:

     -F fref
          Set the reference frequency for causal Q. (fref  =  1.0
          is the default).

      [-K ]
          Use Kjartansson Constant  Q  operator  instead  of  the
          Futterman  operator  (default  =  Futterman).  Causal Q
          means that the layer velocity is complex.


      [-SD]
          Perform wavefield separation at source to consider only
          downgoing waves from the source.  Default no separation
          done at the source.

      [-SU]
          Perform wavefield separation at source to consider only
          upgoing  waves  from the source.  Default no separation
          done at the source.

      [-SPUP]

      [-SPDN]

      [-SSUP]

      [-SSDN]
          These permit wavefield separation  at  the  source  for
          upgoing  and  downgoing  P and upgoing and downgoing S,
          respectively. The values are logically  additive.   The
          -SUP  is equivalent to -SPUP -SSUP together. Default no
          separation done at the source. Due to the nature of the
          wavefield  at  low  frequencies  mathematically correct
          linear trends may be seen if the P-wave contribution is
          desired independently of the corresponding S-wave.

      [-RD]
          Perform wavefield separation at  receiver  to  consider
          only  downgoing  waves  from  the  source.   Default no
          separation done at the receiver.

      [-RU]
          Perform wavefield separation at  receiver  to  consider
          only upgoing waves from the source.  Default no separa-
          tion done at the receiver.

      [-RPUP]

      [-RPDN]

      [-RSUP]

      [-RSDN]
          These permit wavefield separation at the  receiver  for
          upgoing  and  downgoing  P and upgoing and downgoing S,
          respectively. The values are logically  additive.   The
          -SUP  is equivalent to -RPUP -RSUP together. Default no
          separation done at the source. Due to the nature of the
          wavefield  at  low  frequencies  mathematically correct
          linear trends may be seen if the P-wave contribution is
          desired independently of the corresponding S-wave.

          Note that this may not be too useful for  receivers  at
          the  free  surface  since  these  options do not permit
          separation of the wavefield for free surface motion for
          incident  P  or incident S because of the difficulty of
          separating the composite reflected P and S  wavefields.
          The use of these flags is a filtering operation applied
          to the true motion. These receiver flags may have value
          for receivers at depth, as for VSP.

      [-?]
          Do not run the program. Provide online help for command
          line arguments.

MODEL INPUT (UNIT 05)
     The control input for the model and the source and  receiver
     depths is on UNIT LIN:

     -Oname4
          Enter the name of output  file  containing  the  F(f,r)
          spectra. No default.

     -alphaalpha
          Enter the time domain damping.

     -dtdelt
          Enter the sampling interval in ms. No default.

     -tmaxtmax
          Enter the maximum trace time (ms). No default.

     -pminpmin
          Enter the minimum ray parameter. Default = .0000001

     -pmaxpmax
          Enter the maximum ray parameter. No default.

     -pnumpnum
          Enter the number of ray parameters desired. No default.


     -typeieqex
          Enter the control for Green's functions  computed.  The
          choices are as follows:

                    = 0 : earthquake and explosion
                    = 1 : point force and explosion
                    = 2 : earthquake, explosion and point force
                    = 3 : explosion only
                    = 4 : earthquake only
                    = 5 : force only


     -surfisurf
          Enter the surface flag (see explanation below). Default
          = 1

     -halfihalf
          Enter the half  space  flag  (see  explanation  below).
          Default = 0
                                surface   = 0 - elastic
                                            1 - free
                                            2 - rigid
                                halfspace = 0 - elastic
                                            1 - free
                                            2 - rigid
          e.g., surf = 1, half = 1 consists of a model with  both
          surfaces  free,  and   surf = 0, half = 0 consists of a
          model inserted into a wholespace, with the property  of
          the upper halfspace equal to that of the first layer in
          the model, and that of the  lower  halfspace  equal  to
          that of the last layer of the model.

     LINE 1

           read(LIN,*)mmax

     mmax - number of layers in the model

     LINE 2 (mmax lines)

          read(LIN,*)
     (d(i),a(i),b(i),rho(i),qa(i),qb(i),i=1,mmax)

     d    - layer thickness - if  d(mmax)   is  zero,  the  layer
          thickness  will be internal increased to permit sources
          at deeper depths, however, this  may  not  be  what  is
          desired  if  the base boundary condition is not that of
          an elastic halfspace.

     a    - compressional wave velocity

     b    - shear wave velocity
     rho  - layer density

     qa   - compressional wave quality factor,  or  its  inverse,
          for  the  layer.   The  program ASSUMES Q $ >=$ 1. If a
          smaller value is input, it is ASSUMED to be inverse Q.

     qb   - shear wave quality factor, or its  inverse,  for  the
          layer.   The  program  ASSUMES  Q $ >=$ 1. If a smaller
          value is input, it is ASSUMED to be inverse Q.

     LINE 3

           read(LIN,*)mdpths

     mdpths
          - number of source  depths  for  which  synthetics  are
          desired

     LINE 4 (mdpths lines)

           do 2008 i=1,mdpths
                 read(LIN,*)depths(i)
      2008 continue

     depths(i)
          - depth of source

     LINE 5

           read(LIN,*)mdpthr

     mdpthr
          - number of receiver depths for  which  systhetics  are
          desired

     LINE 6 (mdpthr lines)

           do 2009 i=1,mdpthr
                 read(LIN,*)depthr(i)
      2009 continue

     depthr(i)
          - depth of receiver

     LINE 7

THEORY
     (see hardcopy available on request):


DIAGNOSTIC OUTPUT
     Diagnostic output is placed on UNIT 06. The desired  spectra
     are written onto UNIT 04 in unformatted, sequential access.

SEE ALSO
     hspec91(1), rhfoc(1)

SAMPLE DATA SET
Plane Wave Response to Simple Layer Over 1/2-Space
     Below is a script that computes the plane wave resopnse  for
     a  simple  layer  ober  a  1/2-space.  The thickness is 2500
     (units), the p-wave speeds are 8000 and 8500  in  the  layer
     and  1/2-space respectively. The corresponding s-wave speeds
     are 2300 and 3063 and the densities are  2.3  and  2.4.  The
     inverse  p  and s Qs are .001 and The ray parameters go from
     normal incidence to .0005 in 20 steps.


     # SITE RESPONSE

     #!/bin/sh


     cat > hspec91.dat << EOF
         2
      0.2500E+04  0.7000E+04  0.2300E+04  0.2300E+01   0.1000E-02
     0.1000E-02
      0.0000E+00  0.8500E+04  0.3063E+04  0.2400E+01   0.1000E-02
     0.1000E-02
         1
      100
         1
       50
      EOF


     hspec91p  -Oswae  -tmax2000  -dt4   -pmax.0005   -pnum20   <
     hspec91.dat > outssup.sed

     rhfoc -fswae-t -l8 -P \|

     gasp -Ojunk11 -rec













Man(1) output converted with man2html