Source code for bruges.rockphysics.elastic

Elastic impedance.

:copyright: 2015 Agile Geoscience
:license: Apache 2.0
import numpy as np

[docs]def elastic_impedance(vp, vs, rho, theta1, k=None, normalize=False, constants=None, use_sin=False, rho_term=False): """ Returns the elastic impedance as defined by Connolly, 1999; we are using the formulation reported in Whitcombe et al. (2001). Inputs should be shape m x 1, angles should be n x 1. The result will be m x n. Args: vp (ndarray): The P-wave velocity scalar or 1D array. vs (ndarray): The S-wave velocity scalar or 1D array. rho (ndarray): The bulk density scalar or 1D array. theta1 (ndarray): Incident angle(s) in degrees, scalar or array. k (float): A constant, see Connolly (1999). Default is None, which computes it from Vp and Vs. normalize (bool): if True, returns the normalized form of Whitcombe. constants (tuple): A sequence of 3 constants to use for normalization. If you don't provide this, then normalization just uses the means of the inputs. If these are scalars, the result will be the acoustic impedance (see Whitcombe et al., 2001). use_sin (bool): If True, use sin^2 for the first term, instead of tan^2 (see Connolly). rho_term (bool): If True, provide alternative form, with Vp factored out; use in place of density in generating synthetics in other software (see Connolly). In other words, the result can be multipled with Vp to get the elastic impedance. Returns: ndarray: The elastic impedance log at the specficied angle or angles. """ alpha = np.asanyarray(vp, dtype=float) beta = np.asanyarray(vs, dtype=float) rho = np.asanyarray(rho, dtype=float) op = np.sin if use_sin else np.tan k = np.mean(beta**2.0 / alpha**2.0) if k is None else k new_shape = [-1] + alpha.ndim * [1] theta1 = np.radians(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.") a = 1 + op(theta1)**2.0 b = -8 * k * np.sin(theta1)**2.0 c = 1 - 4 * k * np.sin(theta1)**2.0 ei = alpha**a * beta**b * rho**c if normalize: if constants is None: # Use the means; this will return acoustic impedance for scalars. alpha_0 = np.nanmean(vp) beta_0 = np.nanmean(vs) rho_0 = np.nanmean(rho) else: try: alpha_0, beta_0, rho_0 = constants except ValueError as e: raise ValueError("You must provide a sequence of 3 constants.") ei *= alpha_0**(1 - a) * beta_0**(-b) * rho_0**(1 - c) if rho_term: ei /= alpha if ei.size == 1: return np.asscalar(ei) else: return np.squeeze(ei)