# Forward Simulation for Straight Ray Tomography in 2D#

Here we module SimPEG.seismic.straight_ray_tomography to predict arrival time data for a synthetic velocity/slowness model. In this tutorial, we focus on the following:

• How to define the survey

• How to define the forward simulation

• How to predict arrival time data

## Import Modules#

```import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

from discretize import TensorMesh

from SimPEG import maps
from SimPEG.seismic import straight_ray_tomography as tomo
from SimPEG.utils import model_builder

save_file = False
```

## Defining the Survey#

Here, we define survey that will be used for the forward simulation. The survey consists of a horizontal line of point receivers at Y = 100 m and a horizontal line of point sources at Y = -100 m. The shot by each source is measured by all receivers.

```# Define the locations for the sources and receivers.
x = np.linspace(-100, 100, 11)
y_sources = -100 * np.ones(len(x))

source_locations = np.c_[x, y_sources]

# Define the list of receivers used by each source

# Define an empty list to store sources objects. Define each source and
# provide its corresponding receivers list
source_list = []
for ii in range(0, len(y_sources)):
source_list.append(
)

# Define they tomography survey
survey = tomo.Survey(source_list)
```

## Defining a Tensor Mesh#

Here, we create the tensor mesh that will be used to predict arrival time data.

```dh = 10.0  # cell width
N = 21  # number of cells in X and Y direction
hx = [(dh, N)]
hy = [(dh, N)]
mesh = TensorMesh([hx, hy], "CC")
```

## Model and Mapping on Tensor Mesh#

Here, we create the velocity model that will be used to predict the data. Since the physical parameter for straight ray tomography is slowness, we must define a mapping which converts velocity values to slowness values. The model consists of a lower velocity block within a higher velocity background.

```# Define velocity of each unit in m/s
background_velocity = 3000.0
block_velocity = 1500.0

# Define the model. Models in SimPEG are vector arrays.
model = background_velocity * np.ones(mesh.nC)

ind_block = model_builder.getIndicesBlock(np.r_[-50, 20], np.r_[50, -20], mesh.gridCC)
model[ind_block] = block_velocity

# Define a mapping from the model (velocity) to the slowness. If your model
# consists of slowness values, you can use *maps.IdentityMap*.
model_mapping = maps.ReciprocalMap()

# Plot Velocity Model
fig = plt.figure(figsize=(6, 5.5))

ax1 = fig.add_axes([0.15, 0.15, 0.65, 0.75])
mesh.plot_image(model, ax=ax1, grid=True, pcolor_opts={"cmap": "viridis"})
ax1.set_xlabel("x (m)")
ax1.set_ylabel("y (m)")
ax1.plot(x, y_sources, "ro")  # source locations

ax2 = fig.add_axes([0.82, 0.15, 0.05, 0.75])
norm = mpl.colors.Normalize(vmin=np.min(model), vmax=np.max(model))
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.viridis
)
``` ## Simulation: Arrival Time#

Here we demonstrate how to predict arrival time data for the 2D straight ray tomography problem using the 2D Integral formulation.

```# Define the forward simulation. To do this we need the mesh, the survey and
# the mapping from the model to the slowness values on the mesh.
simulation = tomo.Simulation(mesh, survey=survey, slownessMap=model_mapping)

# Compute predicted data for some model
dpred = simulation.dpred(model)
```

## Plotting#

```n_source = len(source_list)

fig = plt.figure(figsize=(8, 5))
obs_string = []

for ii in range(0, n_source):
ax.plot(x, dpred_plotting[:, ii])
obs_string.append("source {}".format(ii + 1))

ax.set_xlim(np.min(x), np.max(x))
ax.set_xlabel("x (m)")
ax.set_ylabel("arrival time (s)")
ax.set_title("Positions vs. Arrival Time")
ax.legend(obs_string, loc="upper right")
``` ```<matplotlib.legend.Legend object at 0x7f44a9aedaf0>
```

## Optional: Exporting Results#

Write the data and true model

```if save_file:
dir_path = os.path.dirname(tomo.__file__).split(os.path.sep)[:-3]
dir_path.extend(["tutorials", "seismic", "assets"])
dir_path = os.path.sep.join(dir_path) + os.path.sep

noise = 0.05 * dpred * np.random.rand(len(dpred))

data_array = np.c_[
np.kron(np.ones(n_source), x),
dpred + noise,
]

fname = dir_path + "tomography2D_data.obs"
np.savetxt(fname, data_array, fmt="%.4e")

output_model = model
fname = dir_path + "true_model_2D.txt"
np.savetxt(fname, output_model, fmt="%.4e")
```

Total running time of the script: (0 minutes 7.191 seconds)

Estimated memory usage: 8 MB

Gallery generated by Sphinx-Gallery