Geometrical Optics

Resources and Code Setup

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)

Part 1: Hamilton-Jacobi Equations and the Eikonal Equation

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.

Hamilton-Jacobi equations

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).

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...).

The Eikonal Equation

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:

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:

  1. 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.
  2. Superimpose some geodesics emanating from \((x_0,z_0)\). What do you notice about the geodesics with respect to the iso-\(\tau\) curves?

  3. 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 for

    t0s = (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 ?

Accessibility