The file `cp7.py`

is provided to you at the computing resources website. This file implements a significant portion of the code you will need for today and tomorrow.

In `cp7.py`

, you will find the lines

```
plt.show()
assert False
```

When you get started, they appear right above '# Problem 2'. The 'sys.exit()' statement terminates the script so that the subsequent incomplete code does not raise errors. As you are progressing through the problem set, you should move these statements right at the end of the current problem you are working on.

Blocks of code that you will need to write are marked within the file by the text `#### FILLOUT ####`

(a single of these statements can mean a few lines of code)

This first session on geometrical optics aims at putting in evidence the relation between the Hamilton-Jacobi equation and the eikonal equation. Later on we will infer some properties of the solutions of the PDEs studied in this class (waves, Helmholtz) based on properties of the sound speed.

Recall the Hamilton-Jacobi equation:

\[ \dot{x} = c(x,z) \cos\theta, \\\\ \dot{z} = c(x,z) \sin\theta, \\\\ \dot\theta = - \cos\theta \partial_z c + \sin\theta \partial_x c = [\nabla c, \vec{\theta}], \\\\ (x,z,\theta)(0) = (x_0,z_0,\theta_0). \]

These equations describe the dynamics of a wave packet with position \((x(t),z(t))\) and velocity \(\vec{v}(t) = c(x(t),z(t))\vec{\theta}(t)\) and preserving the Hamiltonian \(H(X,\Xi) = c(X) |\Xi|\). \(\Xi\) represents the 'slowness covector', dual quantity to the velocity \(\vec{v}\).

**Problem 1:** Visualize the sound speeds of 'gradient' and 'lens' types and infer, based on the dynamical system above, what the trajectories might look like.

**Problem 2:** Using the function `c_dxc_dzc`

, which computes the values of the sound speed as well as its partial derivatives at given coordinates \((x,z)\), complete the functions `geodesic_move`

and `geodesic`

in the file `cp7.py`

to compute solutions of the Hamilton-Jacobi equation. Use explicit Euler scheme (see discussion).

`geodesic_move(x_in, z_in, th_in, dt, Cparams)`

implements one step increment of size`dt`

.`geodesic`

computes a whole geodesic, i.e., the \((x,z,\theta)\) coordinates, between 0 and T with stepsize`dt`

, starting from initial points \((x_0,z_0, \theta_0)\) at \(t=0\).

For further purposes, these functions should be able to handle multiple starting points, i.e. if the variables `x_in`

, `z_in`

and `th_in`

are vectors of size `n_points`

, the geodesic function should return three two-dimensional arrays, each row of which denotes one solution of the Hamilton-Jacobi equation.

When \((x_0,z_0)\) denotes the position of an initial point, you may create the effect of a "point source at \((x_0,z_0)\)" by feeding the `geodesic`

with the set of initial inputs below.

```
nthetas = 60
x0s = x0 * np.ones((nthetas))
z0s = z0 * np.ones((nthetas))
th0s = np.linspace(0.0, 2*np.pi, nthetas)
```

Important note: calling the `plt.plot(X,Y)`

function with `X,Y`

as two 2D arrays will plot the columns of `X`

and `Y`

separately. That is, it is equivalent to

```
for i in xrange(ncol):
plt.plot(X[:,i],Z[:,i])
```

Plot the geodesics on top of the sound speed map.

**Bonus Problem:** Replace your implementation of explicit Euler with a Runge-Kutta 2 scheme (given in class) or any other method of your choice (RK4, etc...).

We now focus on the eikonal equation. In our case, fixing some initial point \(X_0 = (x_0,z_0)\), we want to solve the following problem:

\[ c(X) |\nabla \tau(X) | = 1, \quad X\ne X_0, \\\\ \tau (X_0) = 0. \]

As seen in Prof. Symes' class, there are (at least) two methods to solve this equation:

- by solving the PDE directly,
- by using solutions of the Hamilton-Jacobi equation ('ray tracing'), which make a good approximation until caustic sets are encountered.

The PDE solver is provided and I will briefly explain in class how it is done. The second solver is given by your geodesic function coded above.

**Problem 3:** For the sound speed of 'gradient' type:

- Solve an eikonal problem with the same initial point \((x_0,z_0)\) as above and visualize the solution. Use the matplotlib
`plt.contour`

function for visualizing the \(\tau = \text{constant}\) curves. - Superimpose some geodesics emanating from \((x_0,z_0)\). What do you notice about the geodesics with respect to the iso-\(\tau\) curves?
In another plot, visualize at some fixed \(t\) the locus of all geodesic points at \(t\) ('iso-\(t\)'), and plot on the same figure the iso-\(\tau\) curves corresponding to the same \(t\) (you may look up the

`clabel`

function in the matplotlib doc). Do the experiment fort0s = (0.4, 0.8, 1.2, 1.6, 2.0).

**Problem 4:** Repeat the questions above for the sound speed of 'lens' type. Are *all* the iso-\(t\) and iso-\(\tau\) curves still in good agreement ? If not, what is responsible for the discrepancy ?