Cartesian staggered grid FDTD in Matlab

I have a problem with a technique with particular focus on non-reflecting boundary conditions and handling non-Cartesian surfaces. The balance between numerical efficiency and accuracy, in particular with regard to the approximation of the non-Cartesian surface have to be explored.


My problem is the following two-dimensional problem :
An omnidirectional point source (line in 3D) is located 0.25a above an impedance plane described by a real impedance Z=3c. The impedance plane is bent down as shown above. Left and right of the source the plane is flat over a distance 5a. Then the plane bends upward following a cosine profile over a distance 5a to a height d. The point of interest is located 15a from this cosine profile at 2a from the impedance plane.
The area above the impedance plane is filled with a medium where the wave speed is c.

This area is unbounded upward and towards the left and right hand side, but the impedance plane continues to infinity at the left and at the right.
The frequencies of interest range up to k.a=10 (lower limit k.a=0.1) where k is the wave number =2π/λ.

How can I implement a Cartesian staggered grid FDTD in Matlab to obtain the ratio of the p-field at the point of interest to the p-field at the same distance in open space as a function of frequency for d=0 (flat surface), and d=a?

Would be usefull the FDTD-2D matlab codes of the forum or this is different problem?

Thanks!