<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A quick wedge model
Let's make a quick wedge model and its associated synthetic.
We'll do a zero-offset (aka normal incidence) synthetic, and a full offset one.
Zero offset
We can produce a simple wedge model just by calling the wedge() function.
Step1:  You can then use this integer model to index into an array of rock properties
Step2:  We can use these to make vp and rho earth models. We can use NumPy’s fancy indexing by passing our array of indicies to access the rock properties (in this case acoustic impedance) for every element at once.
Step3:  Each of these new arrays is the shape of the model, but is filled with a rock property
Step4:  Now we can create the reflectivity profile
Step5:  Then make a wavelet and convolve it with the reflectivities
Step6:  The easiest way to check everything worked is probably to plot it.
Step7:  Alternative workflow
In the last example, we made an array of integers, then used indexing to place rock properties in the array, using the index as a sort of look-up.
But we could make the impedance model directly, passing rock properties in to the wedge() function via teh strat argument. It just depends how you want to make your models. 
The strat argument was the default [0, 1, 2] in the last example. Let's pass in the rock properties instead.
Step8:  And look at the result
Step9:  Now the wedge contains rock properties, not integer labels.
Offset reflectivity
Let's make things more realistic by computing offset reflectivities, not just normal incidence (acoustic) reflectivity. We'll need Vs as well
Step10:  We need the model with integers like 0, 1, 2 again
Step11:  Index to get the property models
Step12:  Compute the reflectivity for angles up to 45 degrees
Step13:  The result is three-dimensional
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
import bruges as bg
w, top, base, ref = bg.models.wedge()
plt.imshow(w, interpolation='none')
plt.axvline(ref, color='k', ls='--')
plt.plot(top, 'r-', lw=4)
plt.plot(base, 'r-', lw=4)
plt.show()
Explanation: A quick wedge model
Let's make a quick wedge model and its associated synthetic.
We'll do a zero-offset (aka normal incidence) synthetic, and a full offset one.
Zero offset
We can produce a simple wedge model just by calling the wedge() function.
End of explanation
import numpy as np
vps = np.array([2320, 2350, 2350])
rhos = np.array([2650, 2600, 2620])
Explanation: You can then use this integer model to index into an array of rock properties:
End of explanation
vp = vps[w]
rho = rhos[w]
Explanation: We can use these to make vp and rho earth models. We can use NumPy’s fancy indexing by passing our array of indicies to access the rock properties (in this case acoustic impedance) for every element at once.
End of explanation
vp.shape
vp[:5, :5]
Explanation: Each of these new arrays is the shape of the model, but is filled with a rock property:
End of explanation
rc = bg.reflection.acoustic_reflectivity(vp, rho)
Explanation: Now we can create the reflectivity profile:
End of explanation
ricker, _ = bg.filters.ricker(duration=0.064, dt=0.001, f=40)
syn = bg.filters.convolve(rc, ricker)
syn.shape
Explanation: Then make a wavelet and convolve it with the reflectivities:
End of explanation
fig, axs = plt.subplots(figsize=(17, 4), ncols=5,
                        gridspec_kw={'width_ratios': (4, 4, 4, 1, 4)})
axs[0].imshow(w)
axs[0].set_title('Wedge model')
axs[1].imshow(vp * rho)  
axs[1].set_title('Impedance')
axs[2].imshow(rc)
axs[2].set_title('Reflectivity')
axs[3].plot(ricker, np.arange(ricker.size))
axs[3].axis('off')
axs[3].set_title('Wavelet')
axs[4].imshow(syn)
axs[4].set_title('Synthetic')
axs[4].plot(top, 'w', alpha=0.5)
axs[4].plot(base, 'w', alpha=0.5)
plt.show()
Explanation: The easiest way to check everything worked is probably to plot it.
End of explanation
vps = np.array([2320, 2350, 2350])
rhos = np.array([2650, 2600, 2620])
impedances = vps * rhos
w, top, base, ref = bg.models.wedge(strat=impedances)
Explanation: Alternative workflow
In the last example, we made an array of integers, then used indexing to place rock properties in the array, using the index as a sort of look-up.
But we could make the impedance model directly, passing rock properties in to the wedge() function via teh strat argument. It just depends how you want to make your models. 
The strat argument was the default [0, 1, 2] in the last example. Let's pass in the rock properties instead.
End of explanation
plt.imshow(w, interpolation='none') 
plt.axvline(ref, color='k', ls='--')
plt.plot(top, 'r-', lw=4)
plt.plot(base, 'r-', lw=4)
plt.colorbar()
plt.show()
Explanation: And look at the result:
End of explanation
vps = np.array([2320, 2350, 2350])
vss = np.array([1150, 1250, 1200])
rhos = np.array([2650, 2600, 2620])
Explanation: Now the wedge contains rock properties, not integer labels.
Offset reflectivity
Let's make things more realistic by computing offset reflectivities, not just normal incidence (acoustic) reflectivity. We'll need Vs as well:
End of explanation
w, top, base, ref = bg.models.wedge()
Explanation: We need the model with integers like 0, 1, 2 again:
End of explanation
vp = vps[w]
vs = vss[w]
rho = rhos[w]
Explanation: Index to get the property models:
End of explanation
rc = bg.reflection.reflectivity(vp, vs, rho, theta=range(46))
rc.shape
Explanation: Compute the reflectivity for angles up to 45 degrees:
End of explanation
plt.imshow(rc.real[:, :, 50].T)
Explanation: The result is three-dimensional: the angles are in the first dimension. So the zero-offset reflectivities are in w[0] and 30 degrees is at w[30].
Or, you can slice this cube in another orientation and see how reflectivity varies with angle:
End of explanation
<END_TASK>