# Source code for bruges.reflection.reflection

"""
Various reflectivity algorithms.

"""
from functools import wraps
from collections import namedtuple
import warnings

import numpy as np

from bruges.rockphysics import moduli
from bruges.rockphysics import anisotropy
from bruges.util import deprecated

[docs]def critical_angles(vp1, vp2, vs2=None):
"""Compute critical angle at an interface

Given the upper (vp1) and lower (vp2) velocities at an interface. If you want the PS-wave critical angle as well,
pass vs2 as well.

Args:
vp1 (ndarray): Upper layer P-wave velocity.
vp2 (ndarray): Lower layer P-wave velocity.
vs2 (ndarray): Lower layer S-wave velocity (optional).

Returns:
tuple: The first and second critical angles at the interface, in
degrees. If there isn't a critical angle, it returns np.nan.
"""
ca1 = ca2 = np.nan

if vp1 < vp2:
ca1 = np.degrees(np.arcsin(vp1/vp2))

if (vs2 is not None) and (vp1 < vs2):
ca2 = np.degrees(np.arcsin(vp1/vs2))

return ca1, ca2

[docs]def reflection_phase(reflectivity):
"""
Compute the phase of the reflectivity. Returns an array (or float) of
the phase, in positive multiples of 180 deg or pi rad. So 1 is opposite
phase. A value of 1.1 would be +1.1 \times \pi rad.

Args:
reflectivity (ndarray): The reflectivity, eg from zoeppritz().

Returns:
ndarray: The phase, strictly positive
"""
ph = np.arctan2(np.imag(reflectivity), np.real(reflectivity)) / np.pi
ph[ph == 1] = 0
ph[ph < 0] = 2 + ph[ph < 0]
return ph

[docs]def acoustic_reflectivity(vp=None, rho=None, axis=0, mode='same'):
"""
The acoustic reflectivity, given Vp and RHOB logs. You can also pass
either vp or rho and get velocity or density reflectivity; if you only
have an impedance log or dataset, pass it as vp.

For offset reflectivity, use reflectivity().

Args:
vp (ndarray): The P-wave velocity. At least one of Vp and Rho is required.
rho (ndarray): The bulk density. At least one of Vp and Rho is required.
If you only have Impedance, pass it in as Vp and leave this as None.
It will use a value of 1 to compute impedance and you'll get velocity
reflectivity.
axis (int): The dimension along which to compute the reflectivity.
mode (str): 'same' means the output will be the same shape. 'valid'
means that only strictly valid reflectivities are computed so
the result will be one sample shorted in the axis dimension.

Returns:
ndarray: The reflectivity coefficient series.
"""
if (vp is None) and (rho is None):
raise ValueError("You must pass either vp or rho.")
elif vp is None:
vp = np.ones_like(rho)
elif rho is None:
rho = np.ones_like(vp)

if axis < 0:
axis = vp.ndim + axis

if axis > 0:
imp = np.moveaxis(vp * rho, axis, 0)
else:
imp = vp * rho

if mode == 'same':
pad_width = [(0, 1)] + (imp.ndim - 1) * [(0, 0)]

rc = (imp[1:] - imp[:-1]) / (imp[1:] + imp[:-1])

if axis > 0:
return np.moveaxis(rc, 0, axis)
else:
return rc

[docs]def reflectivity(vp, vs, rho, theta=0, method='zoeppritz_rpp', axis=0, mode='same'):
"""
Offset reflectivity, given Vp, Vs, rho, and offset.

For acoustic reflectivity, use acoustic_reflectivity().

Computes 'upper' and 'lower' intervals from the three provided arrays,
then passes the result to the specified method to compute reflection
coefficients.

For acoustic reflectivity, either use the acoustic_reflectivity()
function, or call reflectivity() passing any log as Vs, e.g. just give
the Vp log twice (it won't be used anyway):

reflectivity(vp, vp, rho)

For anisotropic reflectivity, use either anisotropy.blangy() or
anisotropy.ruger().

Args:
vp (ndarray): The P-wave velocity; float or 1D array length m.
vs (ndarray): The S-wave velocity; float or 1D array length m.
rho (ndarray): The density; float or 1D array length m.
theta (ndarray): The incidence angle; float or 1D array length n.
axis (int): The dimension along which to compute the reflectivity.
mode (str): 'same' means the output will be the same shape. 'valid'
means that only strictly valid reflectivities are computed so
the result will be one sample shorted in the axis dimension.
method (str): The reflectivity equation to use; one of:

- 'scattering_matrix': scattering_matrix
- 'zoeppritz_element': zoeppritz_element
- 'zoeppritz': zoeppritz
- 'zoeppritz_rpp': zoeppritz_rpp
- 'akirichards': akirichards
- 'akirichards_alt': akirichards_alt
- 'fatti': fatti
- 'shuey': shuey
- 'bortfeld': bortfeld
- 'hilterman': hilterman

Notes:

- scattering_matrix gives the full solution
- zoeppritz_element gives a single element which you specify
- zoeppritz returns RPP element only; use zoeppritz_rpp instead
- zoeppritz_rpp is faster than zoeppritz or zoeppritz_element
- You can also pass a function with the same API as these built-in
functions.

Returns:
ndarray. The result of running the specified method on the inputs.
Will be a float (for float inputs and one angle), a 1 x n array
(for float inputs and an array of angles), a 1 x m-1 array (for
float inputs and one angle), or an m-1 x n array (for array inputs
and an array of angles).
"""
methods = {
'scattering_matrix': scattering_matrix,
'zoeppritz_element': zoeppritz_element,
'zoeppritz': zoeppritz,
'zoeppritz_rpp': zoeppritz_rpp,
'akirichards': akirichards,
'akirichards_alt': akirichards_alt,
'fatti': fatti,
'shuey': shuey,
'bortfeld': bortfeld,
'hilterman': hilterman,
}
func = methods.get(method.casefold(), method)

if axis < 0:
axis = vp.ndim + axis

if axis > 0:
vp_  = np.moveaxis(np.asanyarray(vp, dtype=float), axis, 0)
vs_  = np.moveaxis(np.asanyarray(vs, dtype=float), axis, 0)
rho_ = np.moveaxis(np.asanyarray(rho, dtype=float), axis, 0)
else:
vp_, vs_, rho_ = vp, vs, rho

if mode == 'same':
pad_width = [(0, 1)] + (vp.ndim - 1) * [(0, 0)]

rc = func(vp_[:-1], vs_[:-1], rho_[:-1], vp_[1:], vs_[1:], rho_[1:], theta)

if axis > 0:
return np.moveaxis(rc, 0, axis)
else:
return np.squeeze(rc)

[docs]def vectorize(func):
"""
Decorator to make sure the inputs are arrays. We also add a dimension
to theta to make the functions work in an 'outer product' way.

Takes a reflectivity function requiring Vp, Vs, and RHOB for 2 rocks
(upper and lower), plus incidence angle theta, plus kwargs. Returns
that function with the arguments transformed to ndarrays.
"""
@wraps(func)
def wrapper(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, **kwargs):

vp1 = np.asanyarray(vp1, dtype=float)
vs1 = np.asanyarray(vs1, dtype=float) + 1e-12  # Prevent singular matrix.
rho1 = np.asanyarray(rho1, dtype=float)
vp2 = np.asanyarray(vp2, dtype=float)
vs2 = np.asanyarray(vs2, dtype=float) + 1e-12  # Prevent singular matrix.
rho2 = np.asanyarray(rho2, dtype=float)

new_shape = [-1] + vp1.ndim * [1]
theta1 = theta1.reshape(*new_shape)
if (np.nan_to_num(theta1) > np.pi/2.).any():
raise ValueError("Incidence angle theta1 must be less than 90 deg.")

return func(vp1, vs1, rho1, vp2, vs2, rho2, theta1, **kwargs)
return wrapper

[docs]def preprocess(func):
"""
Decorator to preprocess arguments for the reflectivity equations.

Takes a reflectivity function requiring Vp, Vs, and RHOB for 2 rocks
(upper and lower), plus incidence angle theta, plus kwargs. Returns
that function with some arguments transformed.
"""
@wraps(func)
def wrapper(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, **kwargs):

# Interpret tuple for theta1 as a linspace.
if isinstance(theta1, tuple):
if len(theta1) == 2:
start, stop = theta1
theta1 = np.linspace(start, stop, num=stop+1)
elif len(theta1) == 3:
start, stop, step = theta1
steps = (stop / step) + 1
theta1 = np.linspace(start, stop, num=steps)
else:
raise TypeError("Expected 2 or 3 parameters for theta1 expressed as range.")

# Convert theta1 to radians and complex numbers.

return func(vp1, vs1, rho1, vp2, vs2, rho2, theta1, **kwargs)
return wrapper

[docs]@preprocess
@vectorize
def scattering_matrix(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0):
"""
Full Zoeppritz solution, considered the definitive solution.
Calculates the angle dependent p-wave reflectivity of an interface
between two mediums.

Originally written by: Wes Hamlyn, vectorized by Agile.

Returns the complex reflectivity.

Args:
vp1 (float): The upper P-wave velocity.
vs1 (float): The upper S-wave velocity.
rho1 (float): The upper layer's density.
vp2 (float): The lower P-wave velocity.
vs2 (float): The lower S-wave velocity.
rho2 (float): The lower layer's density.
theta1 (ndarray): The incidence angle; float or 1D array length n.

Returns:
ndarray. The exact Zoeppritz solution for all modes at the interface.
A 4x4 array representing the scattering matrix at the incident
angle theta1.
"""
theta1 *= np.ones_like(vp1)
p = np.sin(theta1) / vp1  # Ray parameter.
theta2 = np.arcsin(p * vp2)  # Trans. angle of P-wave.
phi1 = np.arcsin(p * vs1)    # Refl. angle of converted S-wave.
phi2 = np.arcsin(p * vs2)    # Trans. angle of converted S-wave.

# Matrix form of Zoeppritz equations... M & N are matrices.
M = np.array([[-np.sin(theta1), -np.cos(phi1), np.sin(theta2), np.cos(phi2)],
[np.cos(theta1), -np.sin(phi1), np.cos(theta2), -np.sin(phi2)],
[2 * rho1 * vs1 * np.sin(phi1) * np.cos(theta1),
rho1 * vs1 * (1 - 2 * np.sin(phi1) ** 2),
2 * rho2 * vs2 * np.sin(phi2) * np.cos(theta2),
rho2 * vs2 * (1 - 2 * np.sin(phi2) ** 2)],
[-rho1 * vp1 * (1 - 2 * np.sin(phi1) ** 2),
rho1 * vs1 * np.sin(2 * phi1),
rho2 * vp2 * (1 - 2 * np.sin(phi2) ** 2),
-rho2 * vs2 * np.sin(2 * phi2)]])

N = np.array([[np.sin(theta1), np.cos(phi1), -np.sin(theta2), -np.cos(phi2)],
[np.cos(theta1), -np.sin(phi1), np.cos(theta2), -np.sin(phi2)],
[2 * rho1 * vs1 * np.sin(phi1) * np.cos(theta1),
rho1 * vs1 * (1 - 2 * np.sin(phi1) ** 2),
2 * rho2 * vs2 * np.sin(phi2) * np.cos(theta2),
rho2 * vs2 * (1 - 2 * np.sin(phi2) ** 2)],
[rho1 * vp1 * (1 - 2 * np.sin(phi1) ** 2),
-rho1 * vs1 * np.sin(2 * phi1),
- rho2 * vp2 * (1 - 2 * np.sin(phi2) ** 2),
rho2 * vs2 * np.sin(2 * phi2)]])

M_ = np.moveaxis(np.squeeze(M), [0, 1], [-2, -1])
A = np.linalg.inv(M_)
N_ = np.moveaxis(np.squeeze(N), [0, 1], [-2, -1])
Z_ = np.matmul(A, N_)

return np.transpose(Z_, axes=list(range(Z_.ndim - 2)) + [-1, -2])

[docs]def zoeppritz_element(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, element='PdPu'):
"""
Returns any mode reflection coefficients from the Zoeppritz
scattering matrix. Pass in the mode as element, e.g. 'PdSu' for PS.

Wraps scattering_matrix().

Returns the complex reflectivity.

Args:
vp1 (float): The upper P-wave velocity.
vs1 (float): The upper S-wave velocity.
rho1 (float): The upper layer's density.
vp2 (float): The lower P-wave velocity.
vs2 (float): The lower S-wave velocity.
rho2 (float): The lower layer's density.
theta1 (ndarray): The incidence angle; float or 1D array length n.
element (str): The name of the element to return, must be one of:
'PdPu', 'SdPu', 'PuPu', 'SuPu', 'PdSu', 'SdSu', 'PuSu', 'SuSu',
'PdPd', 'SdPd', 'PuPd', 'SuPd', 'PdSd', 'SdSd', 'PuSd', 'SuSd'.

Returns:
ndarray. Array length n of the exact Zoeppritz solution for the
specified modes at the interface, at the incident angle theta1.
"""
elements = np.array([['PdPu', 'SdPu', 'PuPu', 'SuPu'],
['PdSu', 'SdSu', 'PuSu', 'SuSu'],
['PdPd', 'SdPd', 'PuPd', 'SuPd'],
['PdSd', 'SdSd', 'PuSd', 'SuSd']])

Z = scattering_matrix(vp1, vs1, rho1, vp2, vs2, rho2, theta1).T

return np.squeeze(Z[np.where(elements == element)].T)

[docs]def zoeppritz(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0):
"""
Returns the PP reflection coefficients from the Zoeppritz
scattering matrix. Wraps zoeppritz_element().

Returns the complex reflectivity.

Args:
vp1 (float): The upper P-wave velocity.
vs1 (float): The upper S-wave velocity.
rho1 (float): The upper layer's density.
vp2 (float): The lower P-wave velocity.
vs2 (float): The lower S-wave velocity.
rho2 (float): The lower layer's density.
theta1 (ndarray): The incidence angle; float or 1D array length n.

Returns:
ndarray. Array length n of the exact Zoeppritz solution for the
specified modes at the interface, at the incident angle theta1.
"""
return zoeppritz_element(vp1, vs1, rho1, vp2, vs2, rho2, theta1, 'PdPu')

[docs]@preprocess
@vectorize
def zoeppritz_rpp(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0):
"""
Exact Zoeppritz from expression.

This is useful because we can pass arrays to it, which we can't do to
scattering_matrix().

Dvorkin et al. (2014). Seismic Reflections of Rock Properties. Cambridge.

Returns the complex reflectivity.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.

Returns:
ndarray. The exact Zoeppritz solution for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
p = np.sin(theta1) / vp1  # Ray parameter
theta2 = np.arcsin(p * vp2)
phi1 = np.arcsin(p * vs1)  # Reflected S
phi2 = np.arcsin(p * vs2)  # Transmitted S

a = rho2 * (1 - 2 * np.sin(phi2)**2.) - rho1 * (1 - 2 * np.sin(phi1)**2.)
b = rho2 * (1 - 2 * np.sin(phi2)**2.) + 2 * rho1 * np.sin(phi1)**2.
c = rho1 * (1 - 2 * np.sin(phi1)**2.) + 2 * rho2 * np.sin(phi2)**2.
d = 2 * (rho2 * vs2**2 - rho1 * vs1**2)

E = (b * np.cos(theta1) / vp1) + (c * np.cos(theta2) / vp2)
F = (b * np.cos(phi1) / vs1) + (c * np.cos(phi2) / vs2)
G = a - d * np.cos(theta1)/vp1 * np.cos(phi2)/vs2
H = a - d * np.cos(theta2)/vp2 * np.cos(phi1)/vs1

D = E*F + G*H*p**2

rpp = (1/D) * (F*(b*(np.cos(theta1)/vp1) - c*(np.cos(theta2)/vp2)) \
- H*p**2 * (a + d*(np.cos(theta1)/vp1)*(np.cos(phi2)/vs2)))

return np.squeeze(rpp)

[docs]@preprocess
@vectorize
def akirichards(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
The Aki-Richards approximation to the reflectivity.

This is the formulation from Avseth et al., Quantitative seismic
interpretation, Cambridge University Press, 2006. Adapted for a 4-term
formula. See http://subsurfwiki.org/wiki/Aki-Richards_equation.

Returns the complex reflectivity.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.

Returns:
ndarray. The Aki-Richards approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta2 = np.arcsin(vp2/vp1*np.sin(theta1))
drho = rho2-rho1
dvp = vp2-vp1
dvs = vs2-vs1
meantheta = (theta1+theta2) / 2.0
rho = (rho1+rho2) / 2.0
vp = (vp1+vp2) / 2.0
vs = (vs1+vs2) / 2.0

# Compute the coefficients
w = 0.5 * drho/rho
x = 2 * (vs/vp1)**2 * drho/rho
y = 0.5 * (dvp/vp)
z = 4 * (vs/vp1)**2 * (dvs/vs)

# Compute the terms
term1 = w
term2 = -1 * x * np.sin(theta1)**2
term3 = y / np.cos(meantheta)**2
term4 = -1 * z * np.sin(theta1)**2

if terms:
fields = ['term1', 'term2', 'term3', 'term4']
AkiRichards = namedtuple('AkiRichards', fields)
return AkiRichards(np.squeeze([term1 for _ in theta1]),
np.squeeze(term2),
np.squeeze(term3),
np.squeeze(term4)
)
else:
return np.squeeze(term1 + term2 + term3 + term4)

[docs]@preprocess
@vectorize
def akirichards_alt(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
This is another formulation of the Aki-Richards solution.
See http://subsurfwiki.org/wiki/Aki-Richards_equation

Returns the complex reflectivity.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.

Returns:
ndarray. The Aki-Richards approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta2 = np.arcsin(vp2/vp1*np.sin(theta1))
drho = rho2-rho1
dvp = vp2-vp1
dvs = vs2-vs1
theta = (theta1+theta2)/2.0
rho = (rho1+rho2)/2.0
vp = (vp1+vp2)/2.0
vs = (vs1+vs2)/2.0

# Compute the three terms
term1 = 0.5 * (dvp/vp + drho/rho)
term2 = (0.5*dvp/vp-2*(vs/vp)**2*(drho/rho+2*dvs/vs)) * np.sin(theta)**2
term3 = 0.5 * dvp/vp * (np.tan(theta)**2 - np.sin(theta)**2)

if terms:
fields = ['term1', 'term2', 'term3']
AkiRichards = namedtuple('AkiRichards', fields)
return AkiRichards(np.squeeze([term1 for _ in theta1]),
np.squeeze(term2),
np.squeeze(term3)
)
else:
return np.squeeze(term1 + term2 + term3)

[docs]@preprocess
@vectorize
def fatti(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
Compute reflectivities with Fatti's formulation of the Aki-Richards
equation, which does not account for the critical angle. See Fatti et al.
(1994), Geophysics 59 (9). Real numbers only.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.

Returns:
ndarray. The Fatti approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta1 = np.real(theta1)

drho = rho2-rho1
rho = (rho1+rho2) / 2.0
vp = (vp1+vp2) / 2.0
vs = (vs1+vs2) / 2.0
dip = (vp2*rho2 - vp1*rho1)/(vp2*rho2 + vp1*rho1)
dis = (vs2*rho2 - vs1*rho1)/(vs2*rho2 + vs1*rho1)
d = drho/rho

# Compute the three terms
term1 = (1 + np.tan(theta1)**2) * dip
term2 = -8 * (vs/vp)**2 * dis * np.sin(theta1)**2
term3 = -1 * (0.5 * np.tan(theta1)**2 - 2 * (vs/vp)**2 * np.sin(theta1)**2) * d

if terms:
fields = ['term1', 'term2', 'term3']
Fatti = namedtuple('Fatti', fields)
return Fatti(np.squeeze(term1),
np.squeeze(term2),
np.squeeze(term3)
)
else:
return np.squeeze(term1 + term2 + term3)

[docs]@preprocess
@vectorize
def shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0,
terms=False,
"""
Compute Shuey approximation with 3 terms.
http://subsurfwiki.org/wiki/Shuey_equation

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.
return_gradient (bool): Whether to return a tuple of the intercept
and gradient (i.e. the second term divided by sin^2(theta)).

Returns:
ndarray. The Aki-Richards approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta1 = np.real(theta1)

drho = rho2-rho1
dvp = vp2-vp1
dvs = vs2-vs1
rho = (rho1+rho2)/2.0
vp = (vp1+vp2)/2.0
vs = (vs1+vs2)/2.0

# Compute three-term reflectivity
r0 = 0.5 * (dvp/vp + drho/rho)
g = 0.5 * dvp/vp - 2 * (vs**2/vp**2) * (drho/rho + 2 * dvs/vs)
f = 0.5 * dvp/vp

term1 = r0
term2 = g * np.sin(theta1)**2
term3 = f * (np.tan(theta1)**2 - np.sin(theta1)**2)

Shuey = namedtuple('Shuey', fields)
return Shuey(np.squeeze(r0), np.squeeze(g))
elif terms:
fields = ['R0', 'Rg', 'Rf']
Shuey = namedtuple('Shuey', fields)
return Shuey(np.squeeze([term1 for _ in theta1]),
np.squeeze(term2),
np.squeeze(term3)
)
else:
return np.squeeze(term1 + term2 + term3)

def shuey2(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0):
"""
Compute Shuey approximation with 2 terms. Wraps shuey(). Deprecated,
use shuey() instead.
"""
r, g, _ = shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta1=theta1, terms=True)
return r + g

def shuey3(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
Compute Shuey approximation with 3 terms. Wraps shuey(). Deprecated,
use shuey() instead.
"""
return shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta1=theta1)

[docs]@preprocess
@vectorize
def bortfeld(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
Compute Bortfeld approximation with three terms.
http://sepwww.stanford.edu/public/docs/sep111/marie2/paper_html/node2.html
Real numbers only.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.

Returns:
ndarray. The 3-term Bortfeld approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta1 = np.real(theta1)

drho = rho2-rho1
dvp = vp2-vp1
dvs = vs2-vs1
rho = (rho1+rho2)/2.0
vp = (vp1+vp2)/2.0
vs = (vs1+vs2)/2.0
k = (2 * vs/vp)**2
rsh = 0.5 * (dvp/vp - k*drho/rho - 2*k*dvs/vs)

# Compute three-term reflectivity
term1 = 0.5 * (dvp/vp + drho/rho)
term2 = rsh * np.sin(theta1)**2
term3 = 0.5 * dvp/vp * np.tan(theta1)**2 * np.sin(theta1)**2

if terms:
fields = ['term1', 'term2', 'term3']
Bortfeld = namedtuple('Bortfeld', fields)
return Bortfeld(np.squeeze([term1 for _ in theta1]),
np.squeeze(term2),
np.squeeze(term3)
)
else:
return np.squeeze(term1 + term2 + term3)

def bortfeld2(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
The 2-term Bortfeld approximation for ava analysis. Wraps shuey().
Deprecated, use bortfeld() instead.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.

Returns:
ndarray. The 2-term Bortfeld approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta2 = np.arcsin(vp2/vp1*np.sin(theta1))
term1 = 0.5 * np.log((vp2*rho2*np.cos(theta1)) / (vp1*rho1*np.cos(theta2)))
svp2 = (np.sin(theta1)/vp1)**2
dvs2 = (vs1**2-vs2**2)
term2 = svp2 * dvs2 * (2+np.log(rho2/rho1)/np.log(vs2/vs1))

if terms:
return term1, term2
else:
return (term1 + term2)

def bortfeld3(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
return bortfeld(vp1, vs1, rho1, vp2, vs2, rho2, theta1=theta1)

[docs]@preprocess
@vectorize
def hilterman(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, terms=False):
"""
Not recommended, only seems to match Zoeppritz to about 10 deg.

Hilterman (1989) approximation from Mavko et al. Rock Physics Handbook.
According to Dvorkin: "arguably the simplest and a very convenient
[approximation]." At least for small angles and small contrasts. Real
numbers only.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.
terms (bool): Whether or not to return a tuple of the terms of the
equation. The first term is the acoustic impedance.

Returns:
ndarray. The Hilterman approximation for P-P reflectivity at the
interface. Will be a float (for float inputs and one angle), a
1 x n array (for float inputs and an array of angles), a 1 x m
array (for float inputs and one angle), or an n x m array (for
array inputs and an array of angles).
"""
theta1 = np.real(theta1)

ip1 = vp1 * rho1
ip2 = vp2 * rho2
rp0 = (ip2 - ip1) / (ip2 + ip1)

pr2, pr1 = moduli.pr(vp2, vs2), moduli.pr(vp1, vs1)
pravg = (pr2 + pr1) / 2.
pr = (pr2 - pr1) / (1 - pravg)**2.

term1 = rp0 * np.cos(theta1)**2.
term2 = pr * np.sin(theta1)**2.

if terms:
fields = ['term1', 'term2']
Hilterman = namedtuple('Hilterman', fields)
return Hilterman(np.squeeze(term1), np.squeeze(term2))
else:
return np.squeeze(term1 + term2)

[docs]def blangy(vp1, vs1, rho1, vp2, vs2, rho2,
d1=0, e1=0, d2=0, e2=0,
theta1=0):
"""Implements the Blangy equation with the same interface as the other
reflectivity equations. Wraps bruges.anisotropy.blangy(), which you may
prefer to use directly.

Args:
vp1 (ndarray): The upper P-wave velocity; float or 1D array length m.
vs1 (ndarray): The upper S-wave velocity; float or 1D array length m.
rho1 (ndarray): The upper layer's density; float or 1D array length m.
vp2 (ndarray): The lower P-wave velocity; float or 1D array length m.
vs2 (ndarray): The lower S-wave velocity; float or 1D array length m.
rho2 (ndarray): The lower layer's density; float or 1D array length m.
d1 (ndarray): The upper delta; float or 1D array length m.
e1 (ndarray): The upper epsilon; float or 1D array length m.
d2 (ndarray): The lower delta; float or 1D array length m.
e2 (ndarray): The lower epsilon; float or 1D array length m.
theta1 (ndarray): The incidence angle; float or 1D array length n.

Returns:
ndarray. The Blangy approximation for P-P reflectivity at the
interface. Wraps anisotropy.blangy().
"""
_, anisotropic = anisotropy.blangy(vp1, vs1, rho1, d1, e1,  # UPPER
vp2, vs2, rho2, d2, e2,  # LOWER
theta1)
return anisotropic