# Frequency domain Acoustic FWI 

## Acoustic wave equation
We start with the first order acoustic wave-equations expressed in particle velocity $v_z$, $v_x$ and pressure $p$ as
$$
\begin{align}
\frac{\partial v_z}{\partial t} &= b \frac{\partial p}{\partial z} \\ \frac{\partial v_x}{\partial t} &= b \frac{\partial p}{\partial x} \\
\frac{\partial p}{\partial t} &= k ( \frac{\partial v_z}{\partial z} + \frac{\partial v_x}{\partial x}) + s(t)
\end{align}    (1)
$$

## Split PML absorbing boundary condition
To accomdate the split PML boundary conditions, we split the pressure term $p$ into two unphysical terms $p_z$ and $p_x$ with $p=p_z + p_x$. Equation (1) can be reformulated as
$$
\begin{align}
\frac{\partial v_z}{\partial t} + \gamma_z v_z &= b \frac{\partial (p_z+P_x)}{\partial z} \\
\frac{\partial v_x}{\partial t} + \gamma_x v_x &= b \frac{\partial (p_z+p_x)}{\partial x} \\
\frac{\partial p_z}{\partial t} + \gamma_z p_z &= k  \frac{\partial v_z}{\partial z} \\
\frac{\partial p_x}{\partial t} + \gamma_x p_x &= k  \frac{\partial v_x}{\partial x} + s(t) \\
\end{align}    (2)
$$


where $b$ is the bouyancy which equal to $b = \frac{1}{\rho}$, $k$ is the bulk modulus, $\gamma_z$ and $\gamma_x$ are the PML coefficients. The source term is additive, so it can be add to either one of the last two equations of Equation (2).

$\gamma(x) = c_{pml} cos(\frac{\pi}{2} \frac{x}{L})$ where $x$ is the local coordinate in the PML layer, the origin is located at the outer edges of the model, $L$ is the width of PML layer. the scalar $c_{pml}$ is a scalar defined by trial and error. After disrectization of the PML coefficients, both $\gamma_z$ and $\gamma_x$ are a 1D vector.

## Transform to frequency domain
Perform Fourier transform to both side of Equation (2), we can get the wave-equations in frequency domain
$$
\begin{align}
iw v_z + \gamma_z v_z &= b \frac{\partial (p_z+P_x)}{\partial z}  \\
iw v_x + \gamma_x v_x &= b \frac{\partial (p_z+p_x)}{\partial x}  \\
iw p_z + \gamma_z p_z &= k  \frac{\partial v_z}{\partial z}       \\
iw p_x + \gamma_x p_x &= k  \frac{\partial v_x}{\partial x} + s(w)\\
\end{align}    (3)
$$
where $i$ is the imagary unit and $w$ is the radian frequency $w = 2 \pi f$, where $f$ is frequency. Combine the two terms on the left side, we can get
$$
\begin{align}
(iw + \gamma_z) v_z &= b \frac{\partial (p_z+P_x)}{\partial z}  \\
(iw + \gamma_x) v_x &= b \frac{\partial (p_z+p_x)}{\partial x}  \\
(iw + \gamma_z) p_z &= k  \frac{\partial v_z}{\partial z}       \\
(iw + \gamma_x) p_x &= k  \frac{\partial v_x}{\partial x} + s(w)\\
\end{align}    (4)
$$
We define new variable $\epsilon_z = 1+ \frac{\gamma_z}{iw} = 1 - \frac{i \gamma_z}{w}$ and $\epsilon_x = 1+ \frac{\gamma_x}{iw} = 1 - \frac{i \gamma_x}{w}$. As we have discussed for $\gamma_z$ and $\gamma_x$, $\epsilon_z$ and $\epsilon_x$ are also two 1D vectors.


So Equation (4) can be transformed to 
$$
\begin{align}
iw \epsilon_z v_z &= b \frac{\partial (p_z+P_x)}{\partial z}  \\
iw \epsilon_x v_x &= b \frac{\partial (p_z+p_x)}{\partial x}  \\
iw \epsilon_z p_z &= k  \frac{\partial v_z}{\partial z}       \\
iw \epsilon_x p_x &= k  \frac{\partial v_x}{\partial x} + s(w)\\
\end{align} (5)
$$

We can furture transform the Equations (5) into follow form which facilitate to derive the finite-difference stencil.
$$
\begin{align}
iw v_z &= \frac{b}{\epsilon_z} \frac{\partial p}{\partial z}  \\
iw v_x &= \frac{b}{\epsilon_x} \frac{\partial p}{\partial x}  \\
\frac{iw}{k} p_z &= \frac{1}{\epsilon_z} \frac{\partial v_z}{\partial z}       \\
\frac{iw}{k} p_x &= \frac{1}{\epsilon_x} \frac{\partial v_x}{\partial x} + \frac{s(w)}{k \epsilon_x}\\
\end{align} (6) 
$$
Note: the main purpose of moving $\epsilon_z$ and $\epsilon_x$ to the right side is that we can summ the last two equations to recover the pressure componnent $p$.

## Discretization
We use stagger grid to discretize the wave-equation. we place the pressure component $p_z$ and $p_x$, bulk modulus $k$, buoyancy $b$ on integer grid, partical velocity $v_z$ and $v_x$ on half grid along vertical and horizontal direction, respectively. So the discretized field is given as 
$$
v_z[i+\frac{1}{2},j] \\
v_x[i,j+\frac{1}{2}] \\
p_z[i,j],\, p_x[i,j],\, b[i,j],\, k[i,j]
$$ 

With this discretizatino, assume the grid size is $h$ in both vertical and horizontal directions. The partial derivative is approximated by second-order finite-difference method. so we can get
$$
\begin{align}
\frac{\partial v_z}{\partial z}[i,j] &= \frac{v_z[i+\frac{1}{2},j] - v_z[i-\frac{1}{2},j]}{h} \\
\frac{\partial v_x}{\partial x}[i,j] &= \frac{v_x[i,j+\frac{1}{2}] - v_z[i,j-\frac{1}{2}]}{h} \\ 
\end{align}  (7)\\
$$

From the first two equations of Equation (6) and Moving the $iw$ term from left side of Equation (6) to the right hand side, we can get
$$
\begin{align}
v_z[i+\frac{1}{2},j] &= \frac{b[i+\frac{1}{2},j]}{iw\, \epsilon_z[i+\frac{1}{2}]} \frac{1}{h}(p[i+1,j] - p[i,j]) \\
v_z[i-\frac{1}{2},j] &= \frac{b[i-\frac{1}{2},j]}{iw\, \epsilon_z[i-\frac{1}{2}]} \frac{1}{h}(p[i,j] - p[i-1,j]) \\
v_x[i,j+\frac{1}{2}] &= \frac{b[i,j+\frac{1}{2}]}{iw\, \epsilon_x[j+\frac{1}{2}]} \frac{1}{h}(p[i,j+1] - p[i,j]) \\
v_x[i,j-\frac{1}{2}] &= \frac{b[i,j-\frac{1}{2}]}{iw\, \epsilon_x[j-\frac{1}{2}]} \frac{1}{h}(p[i,j] - p[i,j-1])
\end{align}  (8)
$$


Using the expressions in Equation (8) to replace the terms in equation (7), we can get a new system as 
$$
\begin{align}
\frac{\partial v_z}{\partial z}[i,j] &=   \frac{b[i+\frac{1}{2},j]}{iw\, \epsilon_z[i+\frac{1}{2}]} \frac{1}{h^2}(p[i+1,j] - p[i,j]) - \frac{b[i-\frac{1}{2},j]}{iw\, \epsilon_z[i-\frac{1}{2}]} \frac{1}{h^2}(p[i,j] - p[i-1,j])\\
\frac{\partial v_x}{\partial x}[i,j] &= \frac{b[i,j+\frac{1}{2}]}{iw\, \epsilon_x[j+\frac{1}{2}]} \frac{1}{h^2}(p[i,j+1] - p[i,j])   -   \frac{b[i,j-\frac{1}{2}]}{iw\, \epsilon_x[j-\frac{1}{2}]} \frac{1}{h^2}(p[i,j] - p[i,j-1]) \\ 
\end{align}  (9)\\
$$

The last two equations of Equation (6) can be discretized as
$$
\begin{align}
\frac{iw}{k[i,j]} p_z[i,j] &= \frac{1}{\epsilon_z[i]} \frac{\partial v_z}{\partial z}[i,j]       \\
\frac{iw}{k[i,j]} p_x[i,j] &= \frac{1}{\epsilon_x[j]} \frac{\partial v_x}{\partial x}[i,j] + \frac{s(w)}{k[i,j] \epsilon_x[j]}\\
\end{align}  (10)
$$

We know the source is always inject inside the computational area, so $\gamma_x = 0$ and $\epsilon_x = 1$, so equations (10) can be simplified as 
$$
\begin{align}
\frac{iw}{k[i,j]} p_z[i,j] &= \frac{1}{\epsilon_z[i]} \frac{\partial v_z}{\partial z}[i,j]       \\
\frac{iw}{k[i,j]} p_x[i,j] &= \frac{1}{\epsilon_x[j]} \frac{\partial v_x}{\partial x}[i,j] + \frac{s(w)}{k[i,j]}\\
\end{align}  (11)
$$

Using the expressions in equations (9), equations (11) can be represented as
$$
\begin{align}
\frac{iw}{k[i,j]} p_z[i,j] &= \frac{1}{\epsilon_z[i] h^2} \left(  \frac{b[i+\frac{1}{2},j]}{iw\, \epsilon_z[i+\frac{1}{2}]} (p[i+1,j] - p[i,j]) - \frac{b[i-\frac{1}{2},j]}{iw\, \epsilon_z[i-\frac{1}{2}]} (p[i,j] - p[i-1,j])  \right)       \\
\frac{iw}{k[i,j]} p_x[i,j] &= \frac{1}{\epsilon_x[j] h^2} \left( \frac{b[i,j+\frac{1}{2}]}{iw\, \epsilon_x[j+\frac{1}{2}]} (p[i,j+1] - p[i,j])   -   \frac{b[i,j-\frac{1}{2}]}{iw\, \epsilon_x[j-\frac{1}{2}]} (p[i,j] - p[i,j-1]) \right) + \frac{s(w)}{k[i,j]}\\
\end{align}  (11)
$$

Summing these two equations and multiply $iw$ on both side of equations (11), we get
$$
\begin{align}
\frac{-w^2}{k[i,j]} p[i,j] &= \frac{1}{\epsilon_z[i] h^2} \left(  \frac{b[i+\frac{1}{2},j]}{\epsilon_z[i+\frac{1}{2}]} (p[i+1,j] - p[i,j]) - \frac{b[i-\frac{1}{2},j]}{\epsilon_z[i-\frac{1}{2}]} (p[i,j] - p[i-1,j])  \right) \\
&+ \frac{1}{\epsilon_x[j] h^2} \left( \frac{b[i,j+\frac{1}{2}]}{\epsilon_x[j+\frac{1}{2}]} (p[i,j+1] - p[i,j])   -   \frac{b[i,j-\frac{1}{2}]}{\epsilon_x[j-\frac{1}{2}]} (p[i,j] - p[i,j-1]) \right) \\
&+ \frac{iw \, s(w)}{k[i,j]}
\end{align}  (11)
$$

The bouyancy at intermidiate grid point is computed via arithmatic average, for example
$$b[i+\frac{1}{2},j] = \frac{b[i,j] + b[i+1,j]}{2}$$

To facilitate building the Helmhotz operator as a sparse matrix, we can move the pressure component to the right size of the equation and the source term to the left side of the equation, we get 
$$
\begin{align}
&\frac{1}{\epsilon_z[i] h^2} \left(  \frac{b[i+\frac{1}{2},j]}{\epsilon_z[i+\frac{1}{2}]} (p[i+1,j] - p[i,j]) - \frac{b[i-\frac{1}{2},j]}{\epsilon_z[i-\frac{1}{2}]} (p[i,j] - p[i-1,j])  \right) \\
&+ \frac{1}{\epsilon_x[j] h^2} \left( \frac{b[i,j+\frac{1}{2}]}{\epsilon_x[j+\frac{1}{2}]} (p[i,j+1] - p[i,j])   -   \frac{b[i,j-\frac{1}{2}]}{\epsilon_x[j-\frac{1}{2}]} (p[i,j] - p[i,j-1]) \right) \\
&+ \frac{w^2}{k[i,j]} p[i,j] = -\frac{iw \, s(w)}{k[i,j]} 
\end{align}  (12)
$$

So the coefficients in front of
$$
\begin{align}
p[i,j] & \rightarrow \frac{w^2}{k[i,j]} - \frac{1}{\epsilon_z[i] h^2} \left( \frac{b[i+\frac{1}{2},j]}{\epsilon_z[i+\frac{1}{2}]} +  \frac{b[i-\frac{1}{2},j]}{\epsilon_z[i-\frac{1}{2}]} \right) - \frac{1}{\epsilon_x[j] h^2} \left( \frac{b[i,j+\frac{1}{2}]}{\epsilon_x[j+\frac{1}{2}]} + \frac{b[i,j-\frac{1}{2}]}{\epsilon_x[j-\frac{1}{2}]} \right) \\
p[i+1,j] & \rightarrow \frac{1}{\epsilon_z[i] h^2} \frac{b[i+\frac{1}{2},j]}{\epsilon_z[i+\frac{1}{2}]} \\
p[i-1,j] & \rightarrow \frac{1}{\epsilon_z[i] h^2} \frac{b[i-\frac{1}{2},j]}{\epsilon_z[i-\frac{1}{2}]} \\
p[i,j+1] & \rightarrow \frac{1}{\epsilon_x[j] h^2} \frac{b[i,j+\frac{1}{2}]}{\epsilon_x[j+\frac{1}{2}]} \\
p[i,j-1] & \rightarrow \frac{1}{\epsilon_x[j] h^2} \frac{b[i,j-\frac{1}{2}]}{\epsilon_x[j-\frac{1}{2}]} 
\end{align}
$$

## Rotated stagger grid

This method still put the pressure component, the new coordinate system is rotated $45$ degrees.
$$
\begin{align}
iw v_z &= \frac{\sqrt{2}}{2} \frac{b}{\epsilon_z} \left(\frac{\partial p}{\partial z^{\prime}} -\frac{\partial p}{\partial x^{\prime}} \right)  \\
iw v_x &= \frac{\sqrt{2}}{2} \frac{b}{\epsilon_x} \left( \frac{\partial p}{\partial z^{\prime}} + \frac{\partial p}{\partial x^{\prime}}  \right)   \\
\frac{iw}{k} p_z &= \frac{\sqrt{2}}{2\epsilon_z} \left( \frac{\partial v_z}{\partial z^{\prime}} - \frac{\partial v_z}{\partial x^{\prime}} \right)     \\
\frac{iw}{k} p_x &= \frac{\sqrt{2}}{2\epsilon_x} \left( \frac{\partial v_x}{\partial z^{\prime}} + \frac{\partial v_x}{\partial x^{\prime}} \right) + \frac{s(w)}{k}
\end{align}    (13)
$$

We first approximate the spatial partial derivative of partical velocity by finite difference, so we have
$$
\begin{align}
\frac{\partial v_z}{\partial z^{\prime}}[i,j] &= \frac{v_z[i+\frac{1}{2},j+\frac{1}{2}] - v_z[i-\frac{1}{2},j-\frac{1}{2}] }{\sqrt{2}h} \\
\frac{\partial v_z}{\partial x^{\prime}}[i,j] &= \frac{v_z[i-\frac{1}{2},j+\frac{1}{2}] - v_z[i+\frac{1}{2},j-\frac{1}{2}] }{\sqrt{2}h} \\
\frac{\partial v_x}{\partial z^{\prime}}[i,j] &= \frac{v_x[i+\frac{1}{2},j+\frac{1}{2}] - v_x[i-\frac{1}{2},j-\frac{1}{2}] }{\sqrt{2}h} \\
\frac{\partial v_x}{\partial x^{\prime}}[i,j] &= \frac{v_x[i-\frac{1}{2},j+\frac{1}{2}] - v_x[i+\frac{1}{2},j-\frac{1}{2}] }{\sqrt{2}h} \\
\end{align} (14)
$$

From the first equation of the system (13):
$$
iw v_z = \frac{\sqrt{2}}{2} \frac{b}{\epsilon_z} \left(\frac{\partial p}{\partial z^{\prime}} -\frac{\partial p}{\partial x^{\prime}} \right)
$$

we represent particle velocity component $v_z$ by pressure component as
$$
\begin{align}
v_z[i-\frac{1}{2}, j-\frac{1}{2}] &= \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{2h\, iw \epsilon_z[i-\frac{1}{2}]} \left(p[i,j] - p[i-1,j-1] + p[i,j-1] - p[i-1,j] \right) \\
v_z[i+\frac{1}{2}, j-\frac{1}{2}] &= \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{2h\, iw \epsilon_z[i+\frac{1}{2}]} \left(p[i+1,j] - p[i,j-1] + p[i+1,j-1] - p[i,j] \right) \\
v_z[i-\frac{1}{2}, j+\frac{1}{2}] &= \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{2h\, iw \epsilon_z[i-\frac{1}{2}]} \left(p[i,j+1] - p[i-1,j] + p[i,j] - p[i-1,j+1] \right) \\
v_z[i+\frac{1}{2}, j+\frac{1}{2}] &= \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{2h\, iw \epsilon_z[i+\frac{1}{2}]} \left(p[i+1,j+1] - p[i,j] + p[i+1,j] - p[i,j+1] \right) \\
\end{align} (15)
$$

From the second equation of system (13):
$$
iw v_x = \frac{\sqrt{2}}{2} \frac{b}{\epsilon_x} \left( \frac{\partial p}{\partial z^{\prime}} + \frac{\partial p}{\partial x^{\prime}}  \right)   
$$

By taking advantage of this equation, we can also represent $v_x$ via pressure component as:
$$
\begin{align}
v_x[i-\frac{1}{2}, j-\frac{1}{2}] &= \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{2h\, iw \epsilon_x[j-\frac{1}{2}]} \left(p[i,j] - p[i-1,j-1] + p[i-1,j] - p[i,j-1] \right) \\
v_x[i+\frac{1}{2}, j-\frac{1}{2}] &= \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{2h\, iw \epsilon_x[j-\frac{1}{2}]} \left(p[i+1,j] - p[i,j-1] + p[i,j] - p[i+1,j-1] \right) \\
v_x[i-\frac{1}{2}, j+\frac{1}{2}] &= \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{2h\, iw \epsilon_x[j+\frac{1}{2}]} \left(p[i,j+1] - p[i-1,j] + p[i-1,j+1] - p[i,j] \right) \\
v_x[i+\frac{1}{2}, j+\frac{1}{2}] &= \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{2h\, iw \epsilon_x[j+\frac{1}{2}]} \left(p[i+1,j+1] - p[i,j] + p[i,j+1] - p[i+1,j] \right) \\
\end{align} (16)
$$

Summing the last two equations of the system (13), we can get
$$
\frac{iw}{k} p = \frac{\sqrt{2}}{2\epsilon_z} \left( \frac{\partial v_z}{\partial z^{\prime}} - \frac{\partial v_z}{\partial x^{\prime}} \right) + \frac{\sqrt{2}}{2\epsilon_x} \left( \frac{\partial v_x}{\partial z^{\prime}} + \frac{\partial v_x}{\partial x^{\prime}} \right) + \frac{s(w)}{k} \,\,\, (17)
$$

Using the expressions in equation (14) to replace the terms in equation (16), we can get
$$
\begin{align}
\frac{iw}{k[i,j]}p[i,j] & = \frac{1}{2h \epsilon_z[i]} \left( v_z[i+\frac{1}{2}, j+\frac{1}{2}] - v_z[i-\frac{1}{2},j-\frac{1}{2}] + v_z[i+\frac{1}{2}, j-\frac{1}{2}] - v_z[i-\frac{1}{2}, j+\frac{1}{2}] \right) \\
&+ \frac{1}{2h \epsilon_x[j]} \left(v_x[i+\frac{1}{2},j+\frac{1}{2}]- v_x[i-\frac{1}{2},j-\frac{1}{2}] + v_x[i-\frac{1}{2},j+\frac{1}{2}] - v_x[i+\frac{1}{2},j-\frac{1}{2}]\right) + \frac{s}{k[i,j]} \\
\end{align} (18)
$$

The next step is to use the expression in equations (15, 16) to replace the particle velocity terms in above equation, take out the terms $2h$ and $iw$,  then we can get a equation only with pressure component.

$$
\begin{align}
\frac{-w^2}{k[i,j]}p[i,j] &= \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} \left(p[i+1,j+1] - p[i,j] + p[i+1,j] - p[i,j+1] \right) \\
 &- \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} \left(p[i,j] - p[i-1,j-1] + p[i,j-1] - p[i-1,j] \right) \\
 &+ \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} \left(p[i+1,j] - p[i,j-1] + p[i+1,j-1] - p[i,j] \right) \\
 &- \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} \left(p[i,j+1] - p[i-1,j] + p[i,j] - p[i-1,j+1] \right) \\
 &+ \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} \left(p[i+1,j+1] - p[i,j] + p[i,j+1] - p[i+1,j] \right) \\
 &- \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} \left(p[i,j] - p[i-1,j-1] + p[i-1,j] - p[i,j-1] \right) \\
 &+ \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} \left(p[i,j+1] - p[i-1,j] + p[i-1,j+1] - p[i,j] \right) \\
 &- \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} \left(p[i+1,j] - p[i,j-1] + p[i,j] - p[i+1,j-1] \right) \\
 &+ \frac{iw \,S}{K[i,j]}
\end{align} (19)
$$

To facilitate coding, we group the common terms together, we get a new equations
$$
\begin{align}
(\frac{-w^2}{k[i,j]} &+ \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} \\
                          &+ \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} ) p[i,j] \\
                          &- (\frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]}) p[i+1,j+1] \\
                          &- (\frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]})p[i-1,j-1]  \\
                          &- (\frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]})p[i+1,j-1] \\
                          &- (\frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]})p[i-1,j+1] \\
                          &-(\frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} +\frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} -\frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} -\frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} )p[i+1,j]  \\
                          &-(\frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} - \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} )p[i,j+1] \\
                          &-(\frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} - \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} )p[i,j-1] \\
                          &-(\frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} )p[i-1,j] \\
                          &= \frac{iw\,S}{k[i,j]}
\end{align} (20)
$$

## Mixed grid 
By combine the finte-difference grid from cartesian and rotated system, we can get the mixed-grid finite-difference stencial. This process includes two steps:

### Mass acceleration averaging
$$
\begin{align}
\frac{w^2}{k[i,j]}p[i,j] = w^2( &c \cdot \frac{p[i,j]}{k[i,j]} \\ 
+ &d (\frac{p[i+1,j]}{k[i+1,j]}+\frac{p[i-1,j]}{k[i-1,j]}+\frac{p[i,j+1]}{k[i,j+1]}+\frac{p[i,j-1]}{k[i,j-1]}) \\
+ & \frac{1-c-4d}{4} (\frac{p[i+1,j+1]}{k[i+1,j+1]}+\frac{p[i-1,j+1]}{k[i-1,j+1]}+\frac{p[i+1,j-1]}{k[i+1,j-1]}+ \frac{p[i-1,j-1]}{k[i-1,j-1]}) ) \\ 
\end{align}
$$

### Weighted averaging of the two stencil
The weight between the two grid is determined by parameter $\alpha$
$$
\frac{w^2}{k}p + a \cdot \left( \frac{\partial}{\partial z}(b \frac{\partial }{\partial z}) + \frac{\partial}{\partial x}(b \frac{\partial }{\partial x}) \right)p + (1-a) \cdot \left( \frac{\partial}{\partial z^{\prime}}(b\frac{\partial }{\partial z^{\prime}}) + \frac{\partial}{\partial x^{\prime}}(b \frac{\partial }{\partial x^{\prime}}) \right)p  = -\frac{iw\, S}{k}
$$

where $c=0.6287326$, $d=0.092816675$, $\alpha=0.5617365$. so the coefficient $\frac{1-c-4d}{4} = 1.75e-7$.

At here, we give the coefficients for the 9 grid points
$$
\begin{align}
p[i,j] \rightarrow &c \frac{w^2}{k[i,j]} \\
                   -\alpha \left\{  \frac{1}{\epsilon_z[i] h^2}\left(\frac{b[i+\frac{1}{2},j]}{\epsilon_z[i+\frac{1}{2}]} +  \frac{b[i-\frac{1}{2},j]}{\epsilon_z[i-\frac{1}{2}]} \right) + \frac{1}{\epsilon_x[j] h^2} \left( \frac{b[i,j+\frac{1}{2}]}{\epsilon_x[j+\frac{1}{2}]} + \frac{b[i,j-\frac{1}{2}]}{\epsilon_x[j-\frac{1}{2}]}\right) \right\} \\
                   &- (1-\alpha) \left\{ \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} ) \right\} \\
&p[i-1,j-1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i-1,j-1]} + (1-\alpha) \left\{ \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]}  \right\} \\
&p[i,j-1] \rightarrow d \frac{w^2}{k[i,j-1]} + \alpha \frac{1}{\epsilon_x[j] h^2} \frac{b[i,j-\frac{1}{2}]}{\epsilon_x[j-\frac{1}{2}]} + (1-\alpha)\left\{ \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} - \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} \right\} \\
&p[i+1,j-1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i+1,j-1]} + (1-\alpha)\left\{\frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} \right\} \\
&p[i-1,j] \rightarrow d \frac{w^2}{k[i-1,j]} + \alpha \left\{ \frac{1}{\epsilon_z[i] h^2} \frac{b[i-\frac{1}{2},j]}{\epsilon_z[i-\frac{1}{2}]}\right\} + (1-\alpha) \left\{ \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} \right\} \\
&p[i+1,j] \rightarrow d\frac{w^2}{k[i+1,j]} + \alpha \left\{\frac{1}{\epsilon_z[i] h^2} \frac{b[i+\frac{1}{2},j]}{\epsilon_z[i+\frac{1}{2}]}\right\} + (1-\alpha)\left\{\frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} +\frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} -\frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} -\frac{b[i+\frac{1}{2}, j-\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j-\frac{1}{2}]}\right\} \\
&p[i-1,j+1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i-1,j+1]} + (1-\alpha)\left\{\frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]} + \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]}  \right\} \\
&p[i,j+1] \rightarrow d\frac{w^2}{k[i,j+1]} + \alpha \left\{\frac{1}{\epsilon_x[j] h^2} \frac{b[i,j+\frac{1}{2}]}{\epsilon_x[j+\frac{1}{2}]}\right\} + (1-\alpha) \left\{\frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]} - \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} - \frac{b[i-\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i-\frac{1}{2}]}  \right\} \\
&p[i+1,j+1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i+1,j+1]} + (1-\alpha)\left\{\frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_z[i] \epsilon_z[i+\frac{1}{2}]} + \frac{b[i+\frac{1}{2}, j+\frac{1}{2}]}{4h^2 \epsilon_x[j] \epsilon_x[j+\frac{1}{2}]}  \right\}
\end{align}
$$

Outside of PML area, $\epsilon_z=\epsilon_x=1$, the above coefficients can be simplified as 
To facilitate programming, we give the coefficients for each term
$$
\begin{align}
& p[i,j] \rightarrow c\frac{w^2}{k[i,j]} - \frac{a}{h^2}( b[i+\frac{1}{2},j]+ b[i-\frac{1}{2},j]+b[i,j+\frac{1}{2}]+b[i,j-\frac{1}{2}]) - \frac{1-a}{2h^2}(b[i+\frac{1}{2},j+\frac{1}{2}]+b[i-\frac{1}{2},j-\frac{1}{2}]+b[i-\frac{1}{2},j+\frac{1}{2}]+b[i+\frac{1}{2},j-\frac{1}{2}]) \\
& p[i-1,j-1] \rightarrow  \frac{1-c-4d}{4} \frac{w^2}{k[i-1,j-1]} +  \frac{1-a}{2h^2}b[i-\frac{1}{2},j-\frac{1}{2}] \\
& p[i,j-1]   \rightarrow d\frac{w^2}{k[i,j-1]} + \frac{\alpha}{h^2} b[i,j-\frac{1}{2}] \\
& p[i+1,j-1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i+1,j-1]} + \frac{1-a}{2h^2}b[i+\frac{1}{2},j-\frac{1}{2}] \\
& p[i-1,j]   \rightarrow d\frac{w^2}{k[i-1,j]} + \frac{\alpha}{h^2} b[i-\frac{1}{2},j] \\
& p[i+1,j]   \rightarrow d\frac{w^2}{k[i+1,j]} + \frac{\alpha}{h^2} b[i+\frac{1}{2},j] \\
& p[i-1,j+1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i-1,j+1]}+\frac{1-a}{2h^2}b[i-\frac{1}{2},j+\frac{1}{2}] \\
& p[i,j+1]   \rightarrow d\frac{w^2}{k[i,j+1]} + \frac{\alpha}{h^2} b[i,j+\frac{1}{2}] \\
& p[i+1,j+1] \rightarrow \frac{1-c-4d}{4} \frac{w^2}{k[i+1,j+1]} + \frac{1-a}{2h^2}b[i+\frac{1}{2},j+\frac{1}{2}] \\
\end{align}
$$

Use the same way of derivation, we can extend this method to forth order.


```julia

```

255a0e1d675592d0c07d35a85987bf0adabb3a60
ipynb
Jupyter Notebook
docs/.ipynb_checkpoints/Frequency_FWI-checkpoint.ipynb
joaquin-17/SeisAcoustic.jl
98d25280bcf88c471152b2ade977dec754494271
MIT
2019-07-01T15:02:34.000Z
2022-02-08T12:48:43.000Z
docs/.ipynb_checkpoints/Frequency_FWI-checkpoint.ipynb
joaquin-17/SeisAcoustic.jl
98d25280bcf88c471152b2ade977dec754494271
MIT
2022-01-26T22:31:41.000Z
2022-02-05T20:15:21.000Z
docs/.ipynb_checkpoints/Frequency_FWI-checkpoint.ipynb
joaquin-17/SeisAcoustic.jl
98d25280bcf88c471152b2ade977dec754494271
MIT
2020-04-25T02:21:19.000Z
2022-01-11T00:27:24.000Z
Qwen/Qwen-72B
1. YES
2. YES