MATLAB programming for Numerical differenciation in Petroleum Geoscience - Acoustic wave propagation

I am interested in simulating a seismic experiment in a simplest medium possible consisting only two layers separated by a horizontal interface. I have the systems of equilibrum equation and equation of motion fully decribing sound waves in an arbitrary medium produced by the pressure source. The free surface is described by the Dirichet boundary condition on constant pressure equals to zero. Since the scenerios do not have boundaries, the most efficient solution used in the industry is to surround the space by regions of finite width where waves are attenuated in a smart way. I have informations of height, velocity and Density for two different layers of three different models for the purpose of this project. The tasks to be accomplish include; Derive and implement a 2D finite difference solver for the problem. You are free to use whatever numerical differential operators you like based on finite differences. Hint: staggered grid approach both in time and space works very well in this situation. Do not jump onto the inhomogeneous media right away. Test your algorithm on a simple homogeneous medium with an absorbing region instead of free surface. Analyse the method’s order of approximation and stability. There is an analytical solution for the problem in a homogeneous medium. It works best in the frequency domain and is implemented for you in the MATLAB function (I will provide this on request) Study the solution. Choose the lowest wave speed from the proposed and compute the wavefields in the homogeneous model using your FD code. The lower the velocity, the shorter the wavelengths, and hence the stronger the numerical dispersion. Find the spatial sampling that provides the best fit with the analytical solution. Compare both zero and far-offset solutions. A record from one receiver is called a trace. Compare the traces. Remember that the discretized version of Dirac delta function should be scaled by the FD grid cell’s area to provide the correct amplitudes, i.e. the amplitudes of your numerical solution should not depend on the spacing. Do not worry if fit is not perfect, especially in the far-offser recording. Where do the noises come from? Choose the highest wave speed from the proposed and compute the wavefields in the homogenous space using the FD and analytical solution. Find the padding area width and the decay factor that provide the best attenuation for the fields. Describe the attenuation dependence on the time step. Choose the parameters of the top layer in Model 1 and benchmark the numerical solution against the analytical solution with the free surface on. How well is your free surface implemented? Compute the wavefields in Model 1 without the free surface.