Sven Ivansson Feb. 24, 1998
National Defence Research Establishment (FOA)
S-17290 Stockholm, Sweden
email: sveni@sto.foa.se
General introduction to the program rmodfgr
-------------------------------------------
The program rmodfgr belongs to the rpress package. It is intended for
tracking the horizontal slowness of a particular mode as a function of
frequency. To find the modal horizontal slowness at the initial frequency,
one should first run the program rmodfaut. With the modal slowness obtained
at one frequency, the derivative of the slowness with respect to frequency
(which is closely related to the group-slowness) can be used as an initial
prediction for the modal slowness at the next frequency. The required
derivative is computed from derivatives of the dispersion function, see
sections 5.2 and 3.1 of Ref. 7. An accurate slowness value can subsequently
be computed using the secant rule (which is also used for final estimates in
the program rmodfaut). A range-independent (laterally homogeneous)
fluid-solid medium is assumed, with a fluid region interspaced between two
(possibly vanishing) solid regions. The lower and upper boundaries may
independently be either free, rigid, or of "homogeneous half-space" type
(either fluid or solid). In addition, an entirely solid medium (without
fluid) can be handled. Note, for example, that "solid-fluid-solid" media
with surrounding half-spaces of either fluid or solid type can be handled.
The program rmodfgr was a useful tool for verifying the theoretical results
in Ref. 9 numerically.
One should be aware that the mode tracking can fail. Automatic shifts of
Riemann sheets as a branch-cut is approached have not been implemented. One
should look manually for branch-cut crossings and restart the tracking as
desired. (An undesired jump across a branch-cut cannot be excluded in the
automatic procedure.) Furthermore, it is well known that dispersion curves
of different modes may pass very close to each other and that they may even
cross. Even if the frequency steps are chosen adaptively, there is a risk of
jumping to the dispersion curve of another mode during the tracking process.
Hence, it may be useful to check the 'modal slowness' pattern at isolated
frequencies using the program rmodfaut. (No argument-variation checks are
made within rmodfgr.)
Units are s for time, km for length and "10**12 kg" for mass. In particular,
the unit for slowness is s/km.
An underlying harmonic time dependence "exp(-i*omeg*t)" is understood,
omeg>0 being the angular frequency and t being time.
As specified in the information file read.me in this directory rpress, there
are other programs in the rpress package. Some important general remarks on
the implementation are also given there.
Input to the program rmodfgr
----------------------------
FILM
MHQEFF
ICUTTYP
if (icuttyp=3) then
AIMCUT
end if
ISHEETTYP
EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX
EPSFSUB HSTEPFSUBPRE
FREQMIN FREQINT NFR
if ( freqmin<0.0 and freqint=0.0 ) then
OFREQ(1:nfr)
end if
ALADAUER PREALA
IUROT0
UPSTART
FILFU
DSFMIN DSFMAX DSFFAKTMAX
DSUMAX DRELGOALFAKT
RSUG2
REDUCFAKT DDELT XYTOL
Explanation of the parameters
-----------------------------
FILM, character*32, is the name of the infile that specifies the
range-independent medium. The medium is built up by a number of layers,
each of which is divided into a number of sublayers. If so is indicated in
the infile FILM, the number of sublayers for a given layer will be
multiplied by MHQEFF.
Note that frequency-dependent layering and/or characteristic medium
slownesses, as made possible in the other programs by the parameters
"mhq0,freqmhq0,iazimi", is not allowed here.
ICUTTYP specifies the type of branch-cuts to be used if the lower or upper
boundary of the medium is a homogeneous half-space. The values 1 and 2 will
imply hyperbolic and vertical branch-cuts, respectively. (The value 3
will imply a kind of mixture, with a change of type where the imaginary
part of horizontal slowness is +-AIMCUT s/km.)
ISHEETTYP specifies the type of Riemann sheet to be used if the lower or
upper boundary of the medium is a homogeneous half-space. The value 1,2,3,4
will imply the ++,-+,+-,-- sheet, respectively. (Here, the P-wave sheet is
given before the S-wave sheet.) In the presence of two homogeneous
half-spaces, the possibility of different sheets for the lower and upper
half-spaces has not been implemented.
The parameters
EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX
EPSFSUB HSTEPFSUBPRE
are relevant only when numerical integration, as specified in the infile
FILM, is to take place through certain layers. This option has not yet
been fully tested, however, and we omit the detailed presentation of these
parameters.
NFR frequencies OFREQ(1:nfr) Hz are considered in sequence. They can be
given explicitly, or they can be obtained as
OFREQ(jf)=FREQMIN+jf*(FREQINT/NFR) for jf=1,..,nfr.
PREALA is a "temporary reduction factor during the time ALADAUER". For
PREALA less than one, complex frequencies will be used during the
computations.
IUROT0 is used to indicate whether or not given values for horizontal
slowness are to be rotated in the complex plane because of nonreal
frequencies. IUROT0 is irrelevant when PREALA=1.0, which is the normal case
when running rmodfgr. Rotation is indicated by IUROT0=1 as opposed to
IUROT0=0. After rotation, a given real horizontal slowness will correspond
to a real horizontal wavenumber.
One should note that internally, the computations are always done with
"iurot=0" (no rotation). When IUROT0=1, rotation is performed to
get the initial horizontal slowness "ustart" from UPSTART (see below), and
rotation is also performed in connection with the output.
UPSTART is the horizontal slowness (possibly rotated as specified by the
parameter IUROT0) for the desired mode at the initial frequency as given
by FREQMIN or OFREQ(1).
FILFU is the name of the outfile in which three data sets will be provided:
modal horizontal slownesses (possibly rotated as specified by the parameter
IUROT0), derivatives of (unrotated) horizontal slowness with respect to
frequency, and (unrotated) complex group-slownesses (i.e. derivatives of the
modal wavenumbers with respect to angular frequency). For each data set,
results for the different frequences will be provided in sequence.
A good way to study the variation of mode shape with frequency is to run the
program rmodpfw with "nfr=-1" and "mainmodu" as the FILFU outfile from
rmodfgr.
DSFMIN, DSFMAX and DSFFAKTMAX are control parameters for the frequency steps
that are taken during the computations. The frequency step is adjusted
adaptively to achieve a reasonable trade-off between efficiency and safety
in the mode-tracking process.
DSFMIN (typically 0.0 Hz) is the minimal frequency step that is handled.
DSFMAX (typically 1.0 Hz) is the maximal frequency step. If a larger step is
suggested, the next step-size will nevertheless be DSFMAX.
DSFFAKTMAX (typically 5.0) is the maximal step-change factor for two
subsequent steps. If a larger step is suggested, the next step-size will
nevertheless be DSFFAKTMAX times the previous step-size.
DSUMAX and DRELGOALFAKT are further control parameters for the frequency
steps. The derivative of the modal slowness with respect to frequency is
used as an initial prediction for the (unrotated) slowness at the next
frequency. If the implied change of slowness is absolutely greater than
DSUMAX (typically 0.01 s/km), a smaller frequency step is chosen.
DRELGOALFAKT (typically 0.75, always less than 1.0) is a factor for choosing
the next step-size. The largest step-size that is expected to fulfil the
step-size requirements is reduced by multiplication with DRELGOALFAKT.
RSUG2, REDUCFAKT, DDELT and XYTOL are tolerance parameters for the final
determination of the modal horizontal slowness at the next frequency using
the secant rule.
RSUG2 (typically 0.01 s/km) is the half-width for a tolerance square that is
put out centered at the prediction of (unrotated) horizontal slowness
obtained from the derivative of the modal slowness with respect to frequency.
The size of the square is reduced if it is crossed by a branch-cut, but it is
not checked automatically if an undesired jump across a branch-cut has
appeared. If the required zero (modal slowness) is not found within the
square by the secant rule, one goes back and tries with a smaller frequency
step.
REDUCFAKT (less than one) is the required reduction factor in each step for
the dispersion function when the secant rule is applied.
DDELT is the 'horizontal slowness' distance-difference limit for choosing a
second point to initiate the secant rule. (The first point is of course
chosen as the prediction from the derivative of the modal slowness with
respect to frequency.) Since the required derivatives are actually available
at little extra computational work, Newton's method could be used as an
alternative.
XYTOL is the desired accuracy for the modal horizontal slownesses.
Format for the medium infile FILM
---------------------------------
ZSOURC0 ZSOURC1 NS NSHC
NSOURCFW
if (nsourcfw>0) then
NSOURCFW
do i = 1, nsourcfw
NSOFW0(i) AMSOURCFW(i)
end do
end if
AA BB RRHO DBLA DBLB
a number of data blocks of the following type
ZZ IIFSMTYP IIMHQANI
AA BB RRHO DBLA DBLB NSH NH
Explanation of the format for the medium infile FILM
----------------------------------------------------
When NSOURCFW is positive, the file contains data for a vertical source
array (see the information file rpressfw.txt in this directory rpress).
NSOURCFW and NSOFW0() are integers, while AMSOURCFW() are complex numbers.
However, the source information is not relevant for running rmodfgr. The
receiver depths are not relevant either.
In other respects, things are as for a medium infile FILM to rpress, see
the information file rpress.txt in this directory rpress. There is an
additional option, however, to handle an entirely solid medium (without
interior fluid). This option is specified by NS=NSHC=0. In this case, one
must take ZSOURC0=ZSOURC1="the initial ZZ". The specified matching depth
is not relevant for the rmodfgr computations, and the "actual" matching
depth will be at "the initial ZZ".
As mentioned in the information file rmodfaut.txt in this directory rpress,
the number of layers should be kept as small as possible. This will enhance
computational speed.
Format for the outfiles
-----------------------
The outfile FILFU will consist of three data sets. There is an initial line
for each data set, specifying the number of lines with corresponding data
that will follow. Each data line is for a particular frequency that is given
together with a complex number for the desired quantity at that frequency.
Before the three data sets, however, there appear several lines (their
number given as the first integer in the initial line of the file) with
repeated input data and an additional line where the number of data sets
(three) is written.
Subroutine directories
----------------------
The source code for the program and its subroutines can be found in the
following directories:
rpress/cpm (program and main routines)
rpress/cpm/single (routines that have different "multi" versions)
rpress/cpm/dhcomp (routines that have different "full-wave" versions)
rpress/rut (general routines)
rpress/slatec (slatec routines for complex Airy functions)
Example
-------
The following input to rmodfgr was used for tracking the Scholte mode in
Fig. 1c of Ref. 7 from the original frequency 100 Hz down to 1 Hz:
ref7fig1.mx ! film
1 ! mhqeff
2 ! icuttyp
1 ! isheettyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsubpre
100.0 -100.0 100 ! freqmin freqint nfr
1.0 1.0 ! aladauer,preala
0 ! iurot0
(2.39439017180364,4.917060644640056E-03) ! upstart
slowfigfu ! filfu
0.0 1.0 5.0 ! dsfmin,dsfmax,dsffaktmax
0.01 0.75 ! dsumax,drelgoalfakt
0.01 ! rsug2
0.7 0.00001 0.000000001 ! reducfakt,ddelt,xytol
The medium infile ref7fig1.mx has already been presented in the information
file rmodpfw.txt (in this directory rpress), and we do not repeat it here.
To enhance speed, however, it would be advantageous to remove the layer
splitting when rmodfgr is run.
To study the variation of mode shape with frequency, we may then run the
program rmodpfw with the following input, note that "nfr=-1":
ref7fig1.mx ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
2 ! icuttyp
1 ! isheettyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsubpre
0.0 0.0 -1 ! freqmin freqint nfr
1.0 1.0 ! aladauer,preala
0 ! iurot
slowfigfu ! mainmodu
1 ! iuttyp
slowfigp ! mainmodp
1 ! inolx
1 ! iystab
The results are actually very interesting, and they are included in the two
figures that can be printed from the two postscript files slowfig.ps and
slowfigp.ps, respectively, in this directory rpress. In both figures, the
horizontal scale is for frequency and it covers -10 Hz to 110 Hz, with a tick
mark every 10 Hz.
The first figure, the one in slowfig.ps, shows the dispersion curve for our
mode. The vertical slowness scale covers 0 to 4 s/km. Actually, three curves
are shown. The curve for modal horizontal slowness is the one with a minimum
of about 1 s/km at about 11 Hz. The corresponding group-slowness curve is
also shown, with a minimum at about 10 Hz and a maximum at about 18 Hz. These
two curves were obtained by running the program rmodfgr in the way shown
above, tracking our mode from 100 Hz down to 1 Hz.
For a solid-fluid-solid medium of this type, case (iv) of Theorem 6.1 in
Ref. 9 predicts the existence of a slow wave such that the modal slowness
tends to infinity as the inverse 2/3:th power of frequency when the frequency
tends to zero. The third curve in our slowfig.ps figure shows the
corresponding theoretical low-frequency slowness asymptote, with the
constant factor obtained from equation (6.7) of Ref. 9. (As stated in Ref. 9,
this equation is for the case with a rigid lower boundary for our fluid, but
it applies to our actual case with a lower homogeneous solid half-space as
well.) At low frequencies, there is a close agreement to the "rmodfgr" modal
slowness curve, suggesting that our 100 Hz Scholte mode is transformed to
this kind of slow mode when the frequency is lowered.
To study the transformation process, we consider the second figure, the one
in slowfigp.ps. It gives a "waterfall plot" for the modal vertical
displacement, as obtained by running the program rmodpfw in the way shown
above. The vertical scale is for depth, from 0.000 km at the surface down to
0.250 km at the bottom. The Scholte character of the mode persists down to
11 Hz. At lower frequencies, however, the main characteristic of the mode is
an almost constant nonzero vertical displacement within the upper solid (the
upper 0.030 km). This is consistent with the bending-wave appearance of the
slow mode, as predicted by Theorem 6.6 of Ref. 9. Hence, we conclude that the
Scholte wave for the lower fluid-solid interface of our solid-fluid-solid
medium is transformed to a slow bending wave involving the upper solid as the
frequency tends to zero.
Both our figures show the real parts of slowness and vertical displacement,
respectively. However, the corresponding imaginary parts are negligible. For
a physically realistic visco-elastic medium model, absorptions in
dB/wavelength should actually tend to zero as the frequency tends to zero.
In Ref. 9, this was called the "truly" visco-elastic case. With vanishing
anelastic attenuation in the medium, our imaginary parts would vanish
identically.
Using our programs, it turns out to be difficult to follow the mode in this
example to lower frequencies than those of the order of a Hz. Numerical
cancellations emerge among the terms of the propagator matrices, and a
different numerical implementation would be needed. (The Runge-Kutta
procedure implied by nonvanishing "iifsmtyp" in the medium infile "film"
might be useful in this context.)
References
----------
See the information file read.me in this directory rpress.