Computation of amplitudes

To compute the amplitudes, we need to solve equations (14) and (16). In polar coordinates, these two equations become:  

$$\tau_r\gamma_r+{1 \over r^2}\tau_\theta\gamma_\theta=0$$ (24)

and  

$$\tau_r\eta_r+{1 \over r^2}\tau_\theta\eta_\theta=m(x_0,z_0),$$ (25)

respectively. Because the traveltimes are known now, these two equations are first-order linear partial differential equations in $\gamma$ and $\eta$. The algorithms for solving these two equations are very similar to the one for traveltime calculation. To solve equation (24), we replace equations (18) to (20) by

$$f_{i,j} = {\gamma_{i+1,j}-\gamma_{i,j} \over r_j\Delta\theta},$$ (26)

$$g_{i,j} = -{u_{i,j} \over w_{i,j}}f_{i,j}$$ (27)

and  

$$\hat{g}_{i,j} = {\gamma_{i,j+1}-\gamma_{i,j} \over \Delta r},$$ (28)

respectively. As in the calculation of $\hat{w}_{i,j}$, the calculation of $\hat{g}_{i,j}$ also involves four different cases:

1.

$u_{i-1,j} \ge 0$ and $u_{i,j} \le 0$

$$\hat{g}_{i,j}=\left\{ \begin{array} {ll} g_{i-1,j} & \hbox{i... ...\ \hat{w}_{i,j}=s_{i-1,j}\ \hbox{or}\ s_{i,j}\end{array}\right.$$

2.

$u_{i-1,j} \ge 0$ and u_i,j > 0

$$\hat{g}_{i,j}=\left\{ \begin{array} {ll} g_{i-1,j} & \hbox{i... ...-1,j} \\ 0 & \hbox{if}\ \hat{w}_{i,j}=s_{i,j}\end{array}\right.$$

3.

u_i-1,j < 0 and $u_{i,j} \le 0$

$$\hat{g}_{i,j}=\left\{ \begin{array} {ll} g_{i,j} & \hbox{if}... ...,j} \\ 0 & \hbox{if}\ \hat{w}_{i,j}=s_{i-1,j}\end{array}\right.$$

4.

u_i-1,j < 0 and u_i,j > 0

$$\hat{g}_{i,j}=0$$

If $\hat{w}_{i,j}$ and $\hat{g}_{i,j}$ are calculated in one subroutine, then the implementation is simple and the program is fully vectorizable, which is shown as follows:

         do i = 1,n+1
            sign1 = 0.5*(1.+sign(1.,u(i-1)))
            sign2 = 0.5*(1.+sign(1.,-u(i)))
            hatw1 = sign1*min(w(i-1),s(i))+(1.-sign1)*s(i-1)
            hatw2 = sign2*min(w(i),s(i-1))+(1.-sign2)*s(i)
            hatw(i) = min(hatw1,hatw2)
            sign3 = 0.5*(1.+sign(1.,w(i-1)-s(i)))
            sign4 = 0.5*(1.+sign(1.,w(i)-s(i-1)))
            sign5 = 0.5*(1.+sign(1.,hatw2-hatw1))
            hatg1 = sign1*(1.-sign3)*g(i-1)
            hatg2 = sign2*(1.-sign4)*g(i)
            hatg(i) = sign5*hatg1+(1.-sign5)*hatg2
         enddo

From equation (28), we have

$$\gamma_{i,j+1}=\gamma_{i,j}+\hat{g}_{i,j}\Delta r.$$ (29)

I can apply Podvin and Lecomte's method to the computation of function $\gamma(r,\theta)$ so that backward extrapolation can be handled.

To solve equation (25), I use

$$\eta_r = {1 \over \tau_r}[m(x_0,z_0)-{1 \over r^2}\tau_\theta\eta_\theta]$$

to extrapolate the function $\eta(r,\theta)$. Referring to the algorithm for calculating $\tau(r,\theta)$ and $\gamma(r,\theta)$, one can easily construct the algorithm for calculating function $\eta(r,\theta)$.

The Jacobian function in polar coordinates can be determined from (13) as follows:

$$\vert J(r,\theta)\vert = {1 \over \sqrt{\left({\partial \gam... ...t)^2+\left({\partial \gamma \over r\partial \theta}\right)^2}}.$$ (30)

We can either map this Jacobian function to Cartesian coordinates and use equation (11) to calculate the amplitudes, or calculate the amplitudes in polar coordinates and then map them to Cartesian coordinates.