In this sequence of computer problems you will implement a simple example of seismic imaging in the frequency domain.

These exercises are best solved using a single python script file, which I will call `cp4_sol.py`

(but you can name it whatever you wish).

Will will need access to the functions you wrote for the previous computer problems (`cp1_sol.py`

, `cp2_sol.py`

, `cp3_sol.py`

). You can either copy and paste them into your new file or you can use a python `import`

statement to get access to them in your new file.

I suggest you develop your solutions inside this file, and then run the script in an IPython console. This is accomplished by starting an IPython console, from an command prompt, with the command `ipython --pylab`

.

Then, from the IPython console, you can run your script with the command

`%run cp4_sol`

As before,yYou will need a wave speed model for these problems. These are provided in msri_tools.py and can be accessed by downloading the file to the directory you will work in. From an `ipython --pylab`

console you can access the model function by executing the command `from msri_tools import model`

. You can also put the line

`from msri_tools import model`

at the top of your python scripts.

In the frequency domain, the (time-shifted) Ricker wavelet is defined as

\[ \hat w(\nu) = \frac{2}{\sqrt{\pi}} \frac{\nu^2}{\nu_0^3}\exp\{-\frac{\nu^2}{\nu_0^2}\}\exp\{-i2\pi\nu t_0\} \]

where \(\nu\) is the frequency in Hz.

**Problem 1:** Implement a frequency version of your `ricker_wavelet`

function. Plot the absolute value of \(\hat w(\nu)\) for a reasonable range of \(\nu\) and verify that the peak occurs at \(\pm \nu_0\).

**Problem 2:** Derive the Helmholtz (frequency domain) differential equation equivalent to the time domain wave equation from the first problem set. How do the mass, attenuation, and stiffness matrices relate to this equation?

Be careful of the Fourier transform convention. The Ricker wavelet above uses \(\hat f(\nu) = \int_0^T f(t)\exp\{-i2\pi\nu t\}\textrm{d}t\) as the convetion. Then, \(f(t) = \int_{-\infty}^{\infty} \hat f(\nu)\exp\{i2\pi\nu t\}\textrm{d}\nu\) is the inverse transform.

**Problem 3:** Using your `point_source`

and `construct_matrices`

functions, write a function `helmholtz_solver(C, source, nu, config)`

which solves the Helmholtz equation at a single frequency `nu`

, that you developed in problem 2. Write a function `forward_operator_frequency`

which behaves like the time-domain `forward_operator`

function, except that it implements

\[ \hat{\mathcal{F}[m]} = \hat u. \]

Your `record_data`

function should also be useful.

**Problem 4:** Derive the equivalent adjoint state method for the frequency domain. What is the adjoint equation? What are the adjoint sources? Write a function `adjoint_operator_frequency`

which implements the adjoint state method, using your `helmholtz_solver`

,

\[ \hat{F^*}\hat r_\text{ext} = \delta m. \]

Be very careful with your conjugates.

**Problem 5:** Derive the equivalent linear forward model, and using your `helmholtz_solver`

function, implement a `linear_forward_operator_frequency`

function for applying

\[ \hat F\delta m = \hat u_1. \]

**Problem 6:** Verify that the adjoint relationship is satisfied for these operators. Again, be careful with your conjugates.

**Problem 7:** Implement gradient descent for the frequency domain. Why don't you get useful answers using a single frequency (this is exceptional for 1D)? What happens if you formulate the optimization such that you are using multiple frequencies at once? How do your results compare to those from the time domain?

**Bonus Problem 1:** Derive the frequecy domain version of the Ricker wavelet.

**Bonus Problem 2:** Derive how you can implement these same frequency domain operations using only time-domain solvers. Can you implement this? What extra 'tweaks' do you need to do to ensure that your time-based frequency operators pass the adjoint test?

**Bonus Problem 3:** Numerically verify that your time- and frequency-domain solvers satisfy the discrete Plancharel's theorem.