% ============= defs.tex ==================
\gdef\cjklm{{c_{jklm}}}
\gdef\tjk{{\tau_{jk}}}
\gdef\tjkk{{\tau_{jk,k}}}
\gdef\xb {{ \sv x }}
\gdef\sb {{ \sv s }}
\gdef\ab {{ \sv a }}
\gdef\shb {{\hat { \sv s }}}
\gdef\ahb {{\hat { \sv a }}}
\gdef\gb {{ \sv g }}
\gdef\xp {{ ( \xb ) }}
\gdef\xtp {{ ( \xb , t ) }}
\gdef\xfp {{ ( \xb , f ) }}
\gdef\ut {{\tilde u }}
\gdef\Ft {{\tilde F }}
\gdef\at{{\tilde a }}
\gdef\phit{{\tilde \phi }}
\gdef\exp{ \hbox{\rm  exp} }
% ============= elas.tex ==================
%\tracingcommands=1
\def\repno {Elastic notes: W.S. Harlan}
%\myeqad \let\myeqmodel=\myeqnow
\gdef\marginal#1{}

\input report

\startreport
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\input title

%\baselineskip=25pt
\input defs

%\input goals

%\input summary

%\input conclusions

%\input recommendations

%\input introduction

\input prel

\input flow

\input mode

\input modeflow

\input impu

%\input acknowledgements

%\input references

%\input figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

\endreport 

% ============= flow.tex ==================
\majorhead {Energy Flow}
A raypath should indicate the flow and direction of energy in
a wavefield.  Only in isotropic media can we assume that this
flow of energy is perpendicular to wavefronts.   To describe
energy flow, we do not even need to know how to decompose
the wavefield into modes and wavefronts.

Take the dot product of each side of the 
elastic wave equation \myeqwave\ with $\dot u_j$.
\myeqad \let\myeqcons=\myeqnow
$$
\eqalign{
F_j \dot u_j &= 
\rho \ddot u_j \dot u_j - (c_{jklm} u_{l,m} )_{,k} \dot u_j 
\cr
&= {d \over dt} ( {1\over2} \rho \dot u_j \dot u_j )
- ( c_{jklm} u_{l,m} \dot u_j )_{,k}
+ c_{jklm} u_{l,m} \dot u_{j,k} 
\cr
&= {d \over dt} ( {1\over2} \rho \dot u_j \dot u_j
+ {1\over2} c_{jklm} u_{j,k} u_{l,m} )
- (c_{jklm} u_{l,m} \dot u_{k} )_{,j} 
~ . \cr
}
\eqno\myeqcons
$$
\marginal{cons}(The symmetries of equation \myeqsym\ were used.)
This equation is in the form of a conservation law:
\myeqad \let\myeqconsa=\myeqnow
$$
{ {d \sigma \xtp \over dt } + \nabla \cdot {\sv J} \xtp } = 
{ {d \sigma \xtp \over dt } +  J_{j,j} \xtp } = 0,
\eqno\myeqconsa
$$
\marginal{consa}where 
\myeqad \let\myeqdens=\myeqnow
$$
\sigma \xtp \equiv {1\over2} \rho \xp | \dot {\sv u} \xtp |^2  + {1\over2} 
c_{jklm} \xp
u_{j,k} \xtp u_{l,m} \xtp
\eqno\myeqdens
$$
\marginal{dens}is the total energy density, and  
\myeqad \let\myeqpoyn=\myeqnow
$$
J_j \xtp  \equiv - c_{jklm} \xp u_{l,m} \xtp \dot u_j \xtp 
= - \tau_{jk} \xtp \dot u_k \xtp 
\eqno\myeqpoyn
$$
\marginal{poyn}is the energy flow, or Poynting vector.

Define the group velocity (or energy velocity) $g_j$ by
\myeqad \let\myeqgvel=\myeqnow
$$
J_j \xtp \equiv \sigma \xtp g_j \xtp .
\eqno\myeqgvel
$$
\marginal{gvel} The energy flow equals the energy density times
the energy (group) velocity.  Raypaths are generally defined
to follow energy.  Let us define a raypath as a curve which is
tangent to the group velocity vector at every point.
In the anisotropic case, raypaths will not
necessarily be perpendicular to wavefronts.
% ============= impu.tex ==================
\majorhead{Impulsive source}

Assume a body force with an arbitrary wavelet at an impulsive location:
${\sv x}^0$:
\myeqad \let\myeqimp=\myeqnow
$$
F_j \xtp = w(t) b_j \delta ( \xb - \xb^0 ) , \hbox { and } 
\Ft_j \xfp  =  \tilde w (f) \delta ( \xb - \xb^0 ) .
\eqno\myeqimp$$
\marginal{imp}The only term of equation \myeqdifb\ that could cancel
this impulsive term would involve the second spatial derivative
of $a_j \xp$.  Thus,
\myeqad \let\myeqimpu=\myeqnow
$$
c_{jklm} \xp a_{l,km} \xp \phi (t) = - F_j \xtp =
- w(t) b_j \delta ( \xb - \xb^0 ) .
\eqno\myeqimpu$$
\marginal{impu}Since the approximation \myeqein\ allows for an arbitrary
scaling of the amplitude and phase terms, we can assume
that our solution sets
\myeqad \let\myeqphi=\myeqnow
$$
\eqalign {
\phi (t) &= w(t) , \hbox { and } \cr
c_{jklm} \xp a_{l,km} \xp
&= b_j \delta ( \xb - \xb^0 ) .
\cr}
\eqno\myeqphi$$
\marginal{phi}~
% ============= mode.tex ==================
\majorhead {Wave Mode Decomposition}

Assume that a solution for the displacement vector $u_j \xtp $
is well approximated locally by the approximation:
\myeqad \let\myeqein=\myeqnow
$$ u_j \xtp  \approx a_j \xp \phi [ t - T \xp ]
; \hbox { and }
\ut_j \xfp \approx \at_j \xp \phit (f) \exp [ - i 2 \pi f T \xp ] .
\eqno\myeqein$$
\marginal{ein}where $a_j$ gives the local amplitude and
$\phi$ modulates the phase according to the time delay $T$.

Substituting equation \myeqein\ into \myeqwave\ gives the following
(suppressing arguments):
\myeqad \let\myeqdifa=\myeqnow
$$ \eqalign{
0
= &- \rho a_j \ddot \phi + [ c_{jklm} (a_{l,m} \phi + a_l \dot \phi s_m )
]_{,k} + F_j\cr
=  &- \rho a_j \ddot \phi 
+ c_{jklm,k} ( a_{l,m} \phi + a_l \dot \phi s_m ) \cr
 &+ c_{jklm} ( a_{l,km} \phi + a_{l,m} \dot \phi s_k + a_{l,k} \dot \phi
s_m  + a_l \ddot \phi s_k s_m +a_l \dot \phi s_{m,k} ) + F_j\cr
= &- ( \rho a_j - c_{jklm} a_l s_k s_m ) \ddot\phi 
\cr
&+ [c_{jklm,k} a_l s_m + c_{jklm} 
( a_{l,m} s_k + a_{l,k} s_m  + a_l s_{m,k} ) ] \dot\phi
\cr &+ ( c_{jklm,k} a_{l,m} + c_{jklm} a_{l,km} ) \phi + F_j
, \cr}
\eqno\myeqdifa$$
\marginal{difa}where $s_j \xp = T_{,j} \xp $ 
is the phase slowness vector, pointing in the
direction of a local plane wave, with a magnitude equal to the reciprocal
of the wave velocity.

In the frequency domain, we find a similar expression, grouped by
powers of $f$.
\myeqad \let\myeqdifb=\myeqnow
$$ 
\eqalign {
0 = \{ &4 \pi^2 f^2 ( \rho a_j - c_{jklm} a_l s_k s_m ) \cr
+ &i 2 \pi f
[ c_{jklm,k} a_l s_m + c_{jklm} 
( a_{l,m} s_k + a_{l,k} s_m  + a_l s_{m,k} ) ] \cr
+ &(c_{jklm,k} a_{l,m} + c_{jklm} a_{l,km} )  
\} \tilde \phi + \tilde F_j .
\cr
}
\eqno\myeqdifb$$
\marginal{difb}To allow non-vanishing $\phi$, 
each of these three scaled terms should
vanish in areas with a vanishing body force $F_j$.

The first term produces the equation appropriate for a high-frequency
limit---a vector version of the Eikonal equation:
\myeqad \let\myeqeig=\myeqnow
$$ a_j \xp = [ c_{jklm} \xp  s_k \xp s_m \xp / \rho \xp ] 
\cdot a_l \xp ,~~ \forall j .
\eqno\myeqeig$$
\marginal{eig}Notice that this equation is independent of time or frequency.
If we scale equation \myeqeig\ with a phase velocity 
$v \xp  = | {\sv s} \xp |^{-1} $, 
the reciprocal of the magnitude of the phase slowness, then
\myeqad \let\myeqeigen=\myeqnow
$$
v \xp^2 a_j \xp 
= [ c_{jklm} \xp  \hat s_k \xp \hat s_m \xp / \rho \xp ]
\cdot a_l \xp  ,  ~~ \forall j.
\eqno\myeqeigen$$
\marginal{eigen}Choose a unit normal $\hat s_j = v s_j$ 
to a particular plane wave.  The eigenvalues of 
equation \myeqnow\ 
equal the squared velocities, and the
eigenvectors give the polarity of different wave modes (i.e. compressional
and shear waves).

Equations \myeqeig\ and \myeqeigen\ resemble the Christoffel equation,
which Fourier transforms the spatial dimensions and assumes
a homogeneous material.  Rather than assume global homogeneity, I
prefer a high-frequency approximation that assumes stiffness
to change slowly over the spatial wavelengths of the propagating
wavefront.

A plot of $| \sb |$ as a function of $\shb$ is called a slowness surface.
A plot of $v = | \sb |^{-1} $ 
as a function of $\hat s$ is called a normal surface.

The remaining terms of \myeqdifb\ give ``Transport equations'' in regions of vanishing
body sources:
\myeqad \let\myeqtrans=\myeqnow
$$
\eqalign {
[c_{jklm} \xp s_m \xp a_l \xp ]_{,k} + 
c_{jklm} \xp s_k \xp a_{l,m} \xp & =0, \hbox { and }  \cr
[c_{jklm} \xp a_{l,m} \xp]_{,k} &= 0 ,~~ \forall j.
\cr}
\eqno\myeqtrans$$
\marginal{trans}These functions are also independent of time and frequency.

% ============= modeflow.tex ==================
\majorhead {Mode energy flow}
Substituting the approximation \myeqein\ into the total energy
density \myeqdens, we find
\myeqad\let\myeqdensb=\myeqnow
$$
\sigma\xtp = {1\over2} \rho \xp | {\sv a} \xp |^2 
\dot \phi [t - T\xp ]^2
+ {1\over2} c_{jklm} \xp a_j \xp s_k \xp a_l \xp s_m \xp \dot 
\phi [t - T\xp ]^2 ~.
\eqno\myeqdensb
$$
\marginal{densb}Because of equation \myeqeig, these two terms
are equal---i.e., the kinetic equals the potential energy density.
\myeqad\let\myeqdensc=\myeqnow
$$
\sigma\xtp = 
c_{jklm} \xp a_j \xp a_l \xp s_k \xp s_m \xp \dot 
\phi [t - T\xp ]^2 ~.
\eqno\myeqdensc
$$
\marginal{densc}Similarly, the energy flow \myeqpoyn\ becomes
\myeqad\let\myeqpoynb=\myeqnow
$$
J_i \xtp = 
c_{jklm} \xp a_k \xp a_m \xp s_l \xp \dot \phi [t - T\xp ]^2 ~.
\eqno\myeqpoynb
$$
\marginal{poynb}Substituting energy density \myeqdensc\ and 
energy flow \myeqpoynb\
into the definition of group velocity \myeqgvel, we find the
following equality:
\myeqad\let\myeqgroup=\myeqnow
$$
\eqalign {
&c_{jklm} \xp a_k \xp a_m \xp s_l \xp
\equiv
\left[
c_{nklm} \xp a_n \xp a_l \xp s_k \xp s_m \xp 
\right] 
\cdot g_j \xp 
; \hbox {~ and } \cr
&g_j \xp = {
c_{jklm} \xp a_k \xp a_m \xp s_l \xp
\over 
c_{nklm} \xp a_n \xp a_l \xp s_k \xp s_m \xp 
}
= [ c_{jklm} \xp \hat a_k \xp \hat a_m \xp / \rho \xp ] s_l \xp  .
\cr}
\eqno\myeqgroup
$$
\marginal{group}Equation \myeqeig\ has allowed the last equality.
Notice that, under this approximation, the group
velocity is function only of position, not time or frequency.
A plot of $| \sv g |$ as a function of direction $ \hat {\sv g}$
is called a ray surface.

If both sides of equation \myeqgroup\ are dotted with the slowness
vector $ s_j \xp $, then we find that
\myeqad\let\myeqcos=\myeqnow
$$
g_j \xp s_j \xp = 1.
\eqno\myeqcos
$$
\marginal{cos}This important result states that the dot product
of the group velocity vector with the slowness vector is unity.
The slowness vector is perpendicular to the wavefront, by construction,  
but the group velocity vector is not, unless the group and
phase velocities are equal.

To propagate a wavefront one small step $\Delta t$ in time,
we first calculate the normal $\shb$ to each point on
the wavefront.  The phase velocity $v=|\shb|^{-1}$ along
the wavefront is given by the appropriate eigenvalues of 
equation \myeqeig.  With this value of $\sb = \shb / v$, 
we can calculate the group velocity $\gb$ by substitution 
into equation \myeqgroup.  Each point on the wavefront 
can be extrapolated in the direction of $\gb$ by 
a perturbation $\gb \cdot \Delta t$.  Finally,
we connect the revised points on the wavefront 
and recalculate the normal vectors $\shb$.

It is possible to calculate a group
velocity surface from a slowness surface numerically, and
vice versa.  Assuming the elastic material desribed by
the stiffness tensor $\st C$ and the density $\rho$ to be
fixed, let us perturb the slowness vector $\sb + \delta \sb$ 
and particle motion $\ahb + \delta \ahb$ in a way
that continues to satisfy equation \myeqeig.  Expanded to first order,
\myeqad\let\myeqpert=\myeqnow
$$
\delta a_j =
c_{jklm}/\rho \delta s_k s_m a_l
+ c_{jklm}/\rho s_k \delta s_m a_l
+ c_{jklm}/\rho s_k s_m \delta a_l
\eqno\myeqpert
$$
\marginal{pert}Take the dot product of both sides with $a_j$
and regroup,
\myeqad\let\myeqpertu=\myeqnow
$$
a_j \delta a_j  =
2 [ c_{jklm}/\rho s_m a_j a_l ] \delta s_k 
+ [ c_{jklm}/\rho s_k s_m a_j ] \delta a_l
=
2 g_k \delta s_k + a_l \delta a_l .
\eqno\myeqpertu
$$
\marginal{pertu}We have used the definition of group velocity \myeqgroup,
the Eikonal equation \myeqeig, and the symmetries of stiffness \myeqsym\
to simplify the terms in brackets.
Subtracting the identical terms we find the simple result
\myeqad\let\myeqperp=\myeqnow
$$
g_j \delta s_j = 0.
\eqno\myeqperp
$$
\marginal{perp}And because of \myeqcos, we also have 
\myeqad\let\myeqperpe=\myeqnow
$$
s_j \delta g_j = - g_j \delta s_j = 0.
\eqno\myeqperpe
$$
The perturbation $\sb + \delta \sb$ must lie along the slowness
surface to be a valid perturbation.  Thus, $\sb$ is tangent
to the slowness surface, and $\gb$ is {\sl perpendicular} to the slowness
surface.  Similarly, the phase slowness is perpendicular to the ray surface.

The normal relations implied by \myeqperpe\ allow us to calculate
the angle between $\sb$ and $\gb$.  If we know the magnitude of
one vector, then equation \myeqcos\ allows us to calculate the magnitude
of the other vector.
If we have already constructed a slowness surface, then group velocity
directions are calculated as normals to this surface.  The magnitudes
of the group velocities derive from the cosine relation \myeqcos,
providing all information necessary to draw the ray surface.
Similarly, the slowness surface can be constructed from a ray surface.
% ============= prel.tex ==================
\majorhead {Preliminaries}

\myeqad \let\myeqnewton=\myeqnow
Begin with conservation of momentum for a continuous medium:
$$
\rho \xp \cdot \ddot u_j \xtp - \tjkk \xtp = F_j \xtp ,
\eqno\myeqnewton
$$
\marginal{newton}where $t$ is time, and $\xb$  (with elements $x_j$) 
the vector of spatial coordinates.
$\rho$ is the density,  $u_j$ the displacement vector, 
$F_j$ the body force vector,
and $\tau_{jk}$ the stress tensor.  Each dot indicates a time
derivative.  

Tensor notation will be used.  Ordinary subscripts 
index the spatial components of vectors
and tensors.   Subscripts appearing after a comma indicate a spatial
derivative in the indexed direction.  The summation convention
requires that repeated subscripts, like $k$ in equation \myeqnewton,
be implicitly summed and eliminated.

Next, Hooke's Law assumes a linear relationship between the stress
tensor and a spatially differentiated displacement vector, called the strain
tensor:
\myeqad \let\myeqhook=\myeqnow
$$
\tjk \xtp = \cjklm \xp \cdot u_{l,m} \xtp.
\eqno\myeqhook
$$
\marginal{hook}The tensor $\cjklm$ is the elastic stiffness, with a maximum of 21
independent components in a three-dimensional coordinate system.  
The following symmetries always hold:
\myeqad \let\myeqsym=\myeqnow
$$
\eqalign{
c_{jklm} &= c_{kjlm} = c_{lmjk} ; ~ {\rm so} \cr
c_{jklm} &=
c_{kjlm} =
c_{jkml} =
c_{kjml} =
c_{lmjk} =
c_{mljk} =
c_{lmkj} =
c_{mlkj} .\cr}
\eqno\myeqsym
$$
\marginal{sym}Combining the equations \myeqnewton\ and \myeqhook\ to eliminate
the stress tensor produces the elastic wave equation:
\myeqad \let\myeqwave=\myeqnow
$$
\rho \xp \cdot \ddot u_j \xtp 
- [\cjklm \xp u_{l,m} \xtp ]_{,k}
= F_j \xtp  
.
\eqno\myeqwave
$$
\marginal{wave}We use the Fourier kernal $\exp ( i  2 \pi f t )$
to transform from frequency $f$ to time $t$.  The inverse transformation
uses the complex conjugate.
\myeqad \let\myeqwavef=\myeqnow
$$
4 \pi^2 f^2 \rho \xp \ut_j \xfp +
[\cjklm \xp \ut_{l,m} \xfp ]_{,k} + \Ft_j \xfp = 0.
\eqno\myeqwavef
$$
\marginal{wavef}Tildes will mark Fourier transformed functions.
% ============= report.tex ==================
% report.tex
% Old format for Research reports, used by elas.tex
\input twelvept
\input utils

\hsize=6truein
\hoffset=0.3in
\voffset=24pt
\advance\vsize by-\voffset
\nopagenumbers
\parindent=0pt
\parskip=1.5ex plus0.1ex         


%\raggedbottom
%\interlinepenalty=1000

\headline={\ifnum\pageno>1
  {\repno \hfil Page \folio}
  \else
  {\repno \hfil}
  \fi}
% ============= title.tex ==================
\title {Raypath extrapolation of vector wavefields}
\title {in arbitrary elastic media}
\firstauthor {William S. Harlan}
% ============= twelvept.tex ==================
% twelvept.tex
%
\font\twelverm=cmr12 \font\tenrm=cmr10 \font\ninerm=cmr9
\font\eightrm=cmr8 \font\sevenrm=cmr7 \font\sixrm=cmr6
\font\fiverm=cmr5
\font\twelvei=cmmi12 \font\teni=cmmi10 \font\ninei=cmmi9
\font\eighti=cmmi8 \font\seveni=cmmi7 \font\sixi=cmmi6
\font\fivei=cmmi5
\font\twelvesy=cmsy10 \font\tensy=cmsy10 \font\ninesy=cmsy9 % cmsy12 not avail.
\font\eightsy=cmsy8 \font\sevensy=cmsy7 \font\sixsy=cmsy6
\font\fivesy=cmsy5
\font\twelvebf=cmbx12 \font\tenbf=cmbx10 \font\ninebf=cmbx9
\font\eightbf=cmbx8 \font\sevenbf=cmbx7 \font\sixbf=cmbx6
\font\fivebf=cmbx5
\font\twelvett=cmtt12 \font\tentt=cmtt10 \font\ninett=cmtt9
\font\eighttt=cmtt8 \font\seventt=cmtt8 \font\sixtt=cmtt8
\font\fivett=cmtt8 % 7,6,5 not available
\font\twelveit=cmti12 \font\tenit=cmti10 \font\nineit=cmti9
\font\eightit=cmti8 \font\sevenit=cmti7 \font\sixit=cmti7
\font\fiveit=cmti7  % 6 5 not available
\font\twelvesl=cmsl12 \font\tensl=cmsl10 \font\ninesl=cmsl9
\font\eightsl=cmsl8 \font\sevensl=cmsl8 \font\sixsl=cmsl8
\font\fivesl=cmsl8 % 7,6,5 not available
\def\twelvepoint{
  \def\rm{\fam0\twelverm}
  \textfont0=\twelverm \scriptfont0=\eightrm \scriptscriptfont0=\sixrm
  \textfont1=\twelvei  \scriptfont1=\eighti  \scriptscriptfont1=\sixi
  \textfont2=\twelvesy \scriptfont2=\eightsy \scriptscriptfont2=\sixsy
  \textfont3=\tenex    \scriptfont3=\tenex   \scriptscriptfont3=\tenex
  \textfont\itfam=\twelveit  \def\it{\fam\itfam\twelveit}%
  \textfont\slfam=\twelvesl  \def\sl{\fam\slfam\twelvesl}%
  \textfont\ttfam=\twelvett  \def\tt{\fam\ttfam\twelvett}%
  \textfont\bffam=\twelvebf  \scriptfont\bffam\eightbf
  \scriptscriptfont\bffam=\sixbf  \def\bf{\fam\bffam\twelvebf}%
  \normalbaselineskip=13pt  \normalbaselines\rm
  \spaceskip 4pt plus 2pt minus 1.33333pt }
\twelvepoint
\let\rf=\rm  \let\sf=\sl  
\font\bigfont cmssbx10 scaled\magstep2
%\font\smallfont cmss10
\font\smallfont=cmr10 at 10truept
\font\smallerheadfont=cmr10 at 10truept
\font\smallheadfont=cmr12 at 12truept
%\font\largeheadfont=cmdunh10 at 17.28truept
%\font\largeheadfont=cmssbx10 at 17.28truept
\font\largeheadfont=cmr10 at 24.866667truept
\font\logofont=manfnt at 14.4truept
\def\prog#1{{\tenrm #1}}
\def\fivepoint{
  \def\rm{\fam0\fiverm}
  \textfont0=\fiverm \scriptfont0=\fiverm \scriptscriptfont0=\fiverm
  \textfont1=\fivei  \scriptfont1=\fivei  \scriptscriptfont1=\fivei
  \textfont2=\fivesy \scriptfont2=\fivesy \scriptscriptfont2=\fivesy
  \textfont3=\fiverm    \scriptfont3=\fiverm   \scriptscriptfont3=\fiverm
  \textfont\itfam=\fiveit  \def\it{\fam\itfam\fiveit}%
  \textfont\slfam=\fivesl  \def\sl{\fam\slfam\fivesl}%
  \textfont\ttfam=\fivett  \def\tt{\fam\ttfam\fivett}%
  \textfont\bffam=\fivebf  \scriptfont\bffam\fivebf
  \scriptscriptfont\bffam=\fivebf  \def\bf{\fam\bffam\fivebf}%
  \normalbaselineskip=6pt  \normalbaselines\rm
  \spaceskip 2pt plus 1pt minus .666666pt }
% ============= utils.tex ==================
% utils.tex
\def\ifundefined#1{\expandafter\ifx\csname#1\endcsname\relax}
\def\startreport{\bigskip}
\def\title#1{\line {\hfill {\rm #1} \hfill} }
\def\author#1{\bigskip \line {\hfill {\sl #1} \hfill} \bigskip}
\def\firstauthor#1{\bigskip \line {\hfill {\sl #1} \hfill}         }
\def\nextauthor#1{          \line {\hfill {\sl #1} \hfill}         }
\def\lastauthor#1{          \line {\hfill {\sl #1} \hfill} \bigskip}
\def\majorhead#1{\bigskip {\hfill {\bf #1} \hfill} \medskip }
\def\minorhead#1{\bigskip { \bf #1 } \medskip }
\def\endreport{\vfil\eject\end}
% figure captions
\def\figure#1. #2
{\vfil \eject ~  \vfill FIG. #1. #2 }
\def\figurex#1. #2
{\par FIG. #1. #2 }
% reference
\outer\def\reference#1\par{ 
  \hangafter=1 
  \hangindent2em
  #1 
}
% special strings
\def\Geophysics{{\sl Geophysics}}
% equation numbering ; use \myeqad \let\thisequation=\myeqnow
\newcount\myeqno
\myeqno=0
\def\myeqad{\advance\myeqno by1 \xdef\myeqnow{(\number\myeqno)} }
% subheads for memos
\def\iocitem#1{{ \parindent=0pt \vskip 1ex {\tenbf #1} } \rm \nobreak}
\def\hyphen{\discretionary{-}{}{-}}
\def\squeeze{\leftskip=1cm \rightskip=1cm}
\def\unsqueeze{\leftskip=0pt \rightskip=0pt}
  % SEG vector: roman boldface, tilde underneath
\gdef\sv#1{\oalign{{$\bf#1$}\crcr\hidewidth 
\vbox to.2ex{\hbox{\bf\char126}\vss}\hidewidth}} 
 % SEG tensor: roman boldface, two tildes underneath
\gdef\st#1{{\sv{\sv#1}}} 
\def\fillpage{\voffset=-.5in \vsize=9.9in \hoffset=-.5in 
\hsize=7.5in \parindent=0pt}
%\def\marginal#1{{\llap{\fivett #1 ~~~~~}}}
\def\marginal#1{}
\def\capsule{ {\special{psfile=/u/pogrsa/harlaws/docs/tex/input/capsule.ps 
	hscale=20 vscale=20 voffset=-60 hoffset=-20} \hskip 3cm } }
\def\ttywidth{ \hoffset=-.5in \hsize=5.in }
\def\Cerveny{\v{C}erven\'{y}}