Source code for welleng.errors.iscwsa_mwd

import numpy as np
from numpy import sin, cos, tan, pi, sqrt
import yaml
import os
# import imp

# import welleng.error
from ..utils import NEV_to_HLA

# since this is running on different OS flavors
PATH = os.path.dirname(__file__)
FILENAME = os.path.join(
    '', *[PATH, 'tool_codes', 'error_codes.yaml']
)


[docs] class iscwsaMwd:
[docs] def __init__( self, error, model ): """ Class using the ISCWSA MWD (Rev4) model to determine well bore uncertatinty. Parameters ---------- error: an intitiated welleng.error.ErrorModel object Returns ------- errors: welleng.error.ErrorModel object A populated ErrorModel object for the selected error model. """ error.__init__ self.e = error self.errors = {} with open(FILENAME, 'r') as file: iscwsa_error_models = yaml.safe_load(file) self.em = iscwsa_error_models[model] if 'Default Tortusity (rad/m)' in self.em['header']: self.tortuosity = self.em['header']['Default Tortusity (rad/m)'] else: self.tortuosity = None if model == "iscwsa_mwd_rev5": assert self.tortuosity is not None, ( "No default tortuosity defined in model header" ) self._initiate_func_dict() for err in self.em['codes']: # func = self._get_the_func_out(err) func = self.em['codes'][err]['function'] mag = self.em['codes'][err]['magnitude'] propagation = self.em['codes'][err]['propagation'] self.errors[err] = ( self.call_func( code=err, func=func, error=self.e, mag=mag, propagation=propagation, tortuosity=self.tortuosity, ) ) self.cov_NEVs = np.zeros((3, 3, len(self.e.survey_rad))) for _, value in self.errors.items(): self.cov_NEVs += value.cov_NEV self.cov_HLAs = NEV_to_HLA(self.e.survey_rad, self.cov_NEVs)
def _get_the_func_out(self, err): if err in self.exceptional_funcs: func = self.exceptional_funcs[err] else: func = self.em['codes'][err]['function'] return func
[docs] def call_func(self, code, func, error, mag, propagation, **kwargs): """ Function for calling functions by mapping function labels to their functions. """ assert func in self.func_dict, f"no function for function {func}" return self.func_dict[func](code, error, mag, propagation, **kwargs)
def _initiate_func_dict(self): """ This dictionary will need to be updated if/when additional error functions are added to the model. """ self.func_dict = { 'ABXY_TI1': ABXY_TI1, 'ABXY_TI2': ABXY_TI2, 'ABZ': ABZ, 'AMIL': AMIL, 'ASXY_TI1': ASXY_TI1, 'ASXY_TI2': ASXY_TI2, 'ASXY_TI3': ASXY_TI3, 'ASZ': ASZ, 'DBH': DBH, 'AZ': AZ, 'DREF': DREF, 'DSF': DSF, 'DST': DST, 'MBXY_TI1': MBXY_TI1, 'MBXY_TI2': MBXY_TI2, 'MBZ': MBZ, 'MSXY_TI1': MSXY_TI1, 'MSXY_TI2': MSXY_TI2, 'MSXY_TI3': MSXY_TI3, 'MSZ': MSZ, 'SAG': SAG, 'XYM1': XYM1, 'XYM2': XYM2, 'XYM3': XYM3, 'XYM4': XYM4, 'SAGE': SAGE, 'XCL': XCL, # requires an exception 'XYM3L': XYM3L, # looks like there's a mistake in the ISCWSA model 'XYM4L': XYM4L, }
# error functions #
[docs] def DREF(code, error, mag=0.35, propagation='random', NEV=True, **kwargs): dpde = np.full((len(error.survey_rad), 3), [1., 0., 0.]) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def DSF( code, error, mag=0.00056, propagation='systematic', NEV=True, **kwargs ): dpde = np.full((len(error.survey_rad), 3), [1., 0., 0.]) dpde = dpde * np.array(error.survey_rad) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def DST( code, error, mag=0.00000025, propagation='systematic', NEV=True, **kwargs ): dpde = np.full((len(error.survey_rad), 3), [1., 0., 0.]) dpde[:, 0] = error.survey.tvd dpde = dpde * np.array(error.survey_rad) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def ABXY_TI1( code, error, mag=0.0040, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = -cos(error.survey_rad[:, 1]) / error.survey.header.G dpde[:, 2] = ( cos(error.survey_rad[:, 1]) * tan(error.survey.header.dip) * sin(error.survey.azi_mag_rad) ) / error.survey.header.G e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def ABXY_TI2( code, error, mag=0.004, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) with np.errstate(divide='ignore', invalid='ignore'): dpde[:, 2] = np.nan_to_num( ( ( tan(-(error.survey_rad[:, 1]) + (pi/2)) - tan(error.survey.header.dip) * cos(error.survey.azi_mag_rad) ) / error.survey.header.G ), posinf=0.0, neginf=0.0 ) e_DIA = dpde * mag sing = np.where( error.survey_rad[:, 1] < error.survey.header.vertical_inc_limit ) if len(sing[0]) < 1: return error._generate_error(code, e_DIA, propagation, NEV) else: e_NEV = error._e_NEV(e_DIA) n = np.array( 0.5 * error.drdp_sing['double_delta_md'] * -sin(error.drdp_sing['azi2']) * mag ) / error.survey.header.G e = np.array( 0.5 * error.drdp_sing['double_delta_md'] * cos(error.drdp_sing['azi2']) * mag ) / error.survey.header.G v = np.zeros_like(n) e_NEV_sing = np.vstack( ( np.zeros((1, 3)), np.stack((n, e, v), axis=-1), np.zeros((1, 3)) ) ) if error.error_model.split('_')[-1] == 'rev5': e_NEV_sing[1, 1] = ( ( error.survey.md[2] + error.survey.md[1] - 2 * error.survey.md[0] ) / 2 * mag * cos(error.survey.azi_true_rad[1]) / error.survey.header.G ) e_NEV[sing] = e_NEV_sing[sing] e_NEV_star = error._e_NEV_star(e_DIA) n = np.array( 0.5 * error.drdp_sing['delta_md'] * -sin(error.drdp_sing['azi2']) * mag ) / error.survey.header.G e = np.array( 0.5 * error.drdp_sing['delta_md'] * cos(error.drdp_sing['azi2']) * mag ) / error.survey.header.G v = np.zeros_like(n) e_NEV_star_sing = np.vstack( ( np.zeros((1, 3)), np.stack((n, e, v), axis=-1), np.zeros((1, 3)) ) ) if error.error_model.split('_')[-1] == 'rev5': e_NEV_star_sing[1, 1] = ( (error.survey.md[1] - error.survey.md[0]) * mag * ( cos(error.survey.azi_true_rad[1]) / error.survey.header.G ) ) e_NEV_star[sing] = e_NEV_star_sing[sing] return error._generate_error( code, e_DIA, propagation, NEV, e_NEV, e_NEV_star )
[docs] def ABZ(code, error, mag=0.004, propagation='systematic', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = -sin(np.array(error.survey_rad)[:, 1]) / error.survey.header.G dpde[:, 2] = ( sin(np.array(error.survey_rad)[:, 1]) * tan(error.survey.header.dip) * sin(error.survey.azi_mag_rad) ) / error.survey.header.G e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def ASXY_TI1( code, error, mag=0.0005, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = sin( np.array(error.survey_rad)[:, 1] ) * cos(np.array(error.survey_rad)[:, 1]) / sqrt(2) dpde[:, 2] = ( sin(np.array(error.survey_rad)[:, 1]) * -tan(error.survey.header.dip) * cos(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) ) / sqrt(2) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def ASXY_TI2( code, error, mag=0.0005, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = sin( np.array(error.survey_rad)[:, 1] ) * cos(np.array(error.survey_rad)[:, 1]) / 2 dpde[:, 2] = ( sin(np.array(error.survey_rad)[:, 1]) * -tan(error.survey.header.dip) * cos(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) ) / 2 e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def ASXY_TI3( code, error, mag=0.0005, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( sin(np.array(error.survey_rad)[:, 1]) * tan(error.survey.header.dip) * cos(error.survey.azi_mag_rad) - cos(np.array(error.survey_rad)[:, 1])) / 2 e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def ASZ(code, error, mag=0.0005, propagation='systematic', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = ( -sin(np.array(error.survey_rad)[:, 1]) * cos(np.array(error.survey_rad)[:, 1]) ) dpde[:, 2] = ( sin(np.array(error.survey_rad)[:, 1]) * tan(error.survey.header.dip) * cos(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MBXY_TI1( code, error, mag=70.0, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( -cos(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) ) / (error.survey.header.b_total * cos(error.survey.header.dip)) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MBXY_TI2( code, error, mag=70.0, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( cos(error.survey.azi_mag_rad) / ( error.survey.header.b_total * cos(error.survey.header.dip) ) ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MBZ(code, error, mag=70.0, propagation='systematic', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( -sin(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) ) / (error.survey.header.b_total * cos(error.survey.header.dip)) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MSXY_TI1( code, error, mag=0.0016, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( sin(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) * ( tan(error.survey.header.dip) * cos(np.array(error.survey_rad)[:, 1]) + sin(np.array(error.survey_rad)[:, 1]) * cos(error.survey.azi_mag_rad) ) / sqrt(2) ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MSXY_TI2( code, error, mag=0.0016, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( sin(error.survey.azi_mag_rad) * ( tan(error.survey.header.dip) * sin(np.array(error.survey_rad)[:, 1]) * cos(np.array(error.survey_rad)[:, 1]) - cos(np.array(error.survey_rad)[:, 1]) * cos(np.array(error.survey_rad)[:, 1]) * cos(error.survey.azi_mag_rad) - cos(error.survey.azi_mag_rad) ) / 2 ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MSXY_TI3( code, error, mag=0.0016, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( cos(np.array(error.survey_rad)[:, 1]) * cos(error.survey.azi_mag_rad) * cos(error.survey.azi_mag_rad) - cos(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) * sin(error.survey.azi_mag_rad) - tan(error.survey.header.dip) * sin(np.array(error.survey_rad)[:, 1]) * cos(error.survey.azi_mag_rad) ) / 2 e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def MSZ( code, error, mag=0.0016, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = -( sin(np.array(error.survey_rad)[:, 1]) * cos(error.survey.azi_mag_rad) + tan(error.survey.header.dip) * cos(np.array(error.survey_rad)[:, 1]) ) * sin(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def AZ(code, error, mag=0.00628, propagation='systematic', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = 1 e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def DBH( code, error, mag=np.radians(5000), propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = 1 / ( error.survey.header.b_total * cos(error.survey.header.dip) ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def DBHR( code, error, mag=np.radians(3000), propagation='random', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = 1 / ( error.survey.header.b_total * cos(error.survey.header.dip) ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def AMIL(code, error, mag=220.0, propagation='systematic', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = ( -sin(np.array(error.survey_rad)[:, 1]) * sin(error.survey.azi_mag_rad) / (error.survey.header.b_total * cos(error.survey.header.dip)) ) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def SAG( code, error, mag=0.00349, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = sin(np.array(error.survey_rad)[:, 1]) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def SAGE( code, error, mag=0.00175, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = sin(np.array(error.survey_rad)[:, 1]) ** 0.25 e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def XYM1( code, error, mag=0.00175, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = np.absolute(sin(np.array(error.survey_rad)[:, 1])) e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def XYM2( code, error, mag=0.00175, propagation='systematic', NEV=True, **kwargs ): propagation = 'systematic' # incorrect in the rev5 model tab dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 2] = -1 e_DIA = dpde * mag return error._generate_error(code, e_DIA, propagation, NEV)
[docs] def XYM3( code, error, mag=0.00175, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = ( np.absolute(cos(np.array(error.survey_rad)[:, 1])) * cos(error.survey.azi_true_rad) ) with np.errstate(divide='ignore', invalid='ignore'): dpde[:, 2] = np.nan_to_num( -( np.absolute(cos(np.array(error.survey_rad)[:, 1])) * sin(error.survey.azi_true_rad) ) / sin(np.array(error.survey_rad)[:, 1]), posinf=0.0, neginf=0.0 ) e_DIA = dpde * mag sing = np.where( error.survey_rad[:, 1] < error.survey.header.vertical_inc_limit ) if len(sing[0]) < 1: return error._generate_error(code, e_DIA, propagation, NEV) else: e_NEV = error._e_NEV(e_DIA) n = np.array(0.5 * error.drdp_sing['double_delta_md'] * mag) e = np.zeros(len(error.drdp_sing['double_delta_md'])) v = np.zeros_like(n) e_NEV_sing = np.vstack( ( np.zeros((1, 3)), np.stack((n, e, v), axis=-1), np.zeros((1, 3)) ) ) e_NEV[sing] = e_NEV_sing[sing] e_NEV_star = error._e_NEV_star(e_DIA) n = np.array(0.5 * error.drdp_sing['delta_md'] * mag) e = np.zeros(len(error.drdp_sing['delta_md'])) v = np.zeros_like(n) e_NEV_star_sing = np.vstack( ( np.zeros((1, 3)), np.stack((n, e, v), axis=-1), np.zeros((1, 3)) ) ) e_NEV_star[sing] = e_NEV_star_sing[sing] return error._generate_error( code, e_DIA, propagation, NEV, e_NEV, e_NEV_star )
[docs] def XYM4( code, error, mag=0.00175, propagation='systematic', NEV=True, **kwargs ): dpde = np.zeros((len(error.survey_rad), 3)) dpde[:, 1] = np.absolute( cos(np.array(error.survey_rad)[:, 1]) ) * sin(error.survey.azi_true_rad) with np.errstate(divide='ignore', invalid='ignore'): dpde[:, 2] = np.nan_to_num( ( np.absolute(np.cos(np.array(error.survey_rad)[:, 1])) * cos(error.survey.azi_true_rad) ) / sin(np.array(error.survey_rad)[:, 1]), posinf=0.0, neginf=0.0 ) e_DIA = dpde * mag sing = np.where( error.survey_rad[:, 1] < error.survey.header.vertical_inc_limit ) if len(sing[0]) < 1: return error._generate_error(code, e_DIA, propagation, NEV) else: e_NEV = error._e_NEV(e_DIA) n = np.zeros(len(error.drdp_sing['double_delta_md'])) e = np.array(0.5 * error.drdp_sing['double_delta_md'] * mag) v = np.zeros_like(n) e_NEV_sing = np.vstack( ( np.zeros((1, 3)), np.stack((n, e, v), axis=-1), np.zeros((1, 3)) ) ) e_NEV[sing] = e_NEV_sing[sing] e_NEV_star = error._e_NEV_star(e_DIA) n = np.zeros(len(error.drdp_sing['delta_md'])) e = np.array(0.5 * error.drdp_sing['delta_md'] * mag) v = np.zeros_like(n) e_NEV_star_sing = np.vstack( ( np.zeros((1, 3)), np.stack((n, e, v), axis=-1), np.zeros((1, 3)) ) ) e_NEV_star[sing] = e_NEV_star_sing[sing] return error._generate_error( code, e_DIA, propagation, NEV, e_NEV, e_NEV_star )
[docs] def XCL(code, error, mag=0.0167, propagation='random', NEV=True, **kwargs): """ Dummy function to manage the ISCWSA workbook not correctly defining the weighting functions. """ tortuosity = kwargs['tortuosity'] if code == "XCLA": return XCLA( code, error, mag=mag, propagation=propagation, NEV=NEV, tortuosity=tortuosity ) else: return XCLH( code, error, mag=mag, propagation=propagation, NEV=NEV, tortuosity=tortuosity )
[docs] def XCLA(code, error, mag=0.0167, propagation='random', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) def manage_sing(error, kwargs): temp = np.absolute( sin(error.survey.inc_rad[1:]) * ((( error.survey.azi_true_rad[1:] - error.survey.azi_true_rad[:-1] + pi ) % (2 * pi)) - pi) ) temp[np.where( error.survey.inc_rad[:-1] < error.survey.header.vertical_inc_limit )] = 0 return temp dpde[1:, 0] = ( (error.survey.md[1:] - error.survey.md[0:-1]) * np.amax(np.stack(( manage_sing(error, kwargs), ( kwargs['tortuosity'] * (error.survey.md[1:] - error.survey.md[0:-1]) ) ), axis=-1), axis=-1) * -sin(error.survey.azi_true_rad[1:]) ) dpde[1:, 1] = ( (error.survey.md[1:] - error.survey.md[0:-1]) * np.amax(np.stack(( manage_sing(error, kwargs), ( kwargs['tortuosity'] * (error.survey.md[1:] - error.survey.md[0:-1]) ) ), axis=-1), axis=-1) * cos(error.survey.azi_true_rad[1:]) ) e_DIA = dpde * mag return error._generate_error( code, e_DIA, propagation, NEV, e_NEV=e_DIA, e_NEV_star=e_DIA )
[docs] def XCLH(code, error, mag=0.0167, propagation='random', NEV=True, **kwargs): dpde = np.zeros((len(error.survey_rad), 3)) dpde[1:, 0] = ( (error.survey.md[1:] - error.survey.md[0:-1]) * np.amax(np.stack(( np.absolute( (error.survey.inc_rad[1:] - error.survey.inc_rad[:-1]) ), ( kwargs['tortuosity'] * (error.survey.md[1:] - error.survey.md[0:-1]) ) ), axis=-1), axis=-1) * cos(error.survey.inc_rad[1:]) * cos(error.survey.azi_true_rad[1:]) ) dpde[1:, 1] = ( (error.survey.md[1:] - error.survey.md[0:-1]) * np.amax(np.stack(( np.absolute( (error.survey.inc_rad[1:] - error.survey.inc_rad[:-1]) ), ( kwargs['tortuosity'] * (error.survey.md[1:] - error.survey.md[0:-1]) ) ), axis=-1), axis=-1) * cos(error.survey.inc_rad[1:]) * sin(error.survey.azi_true_rad[1:]) ) dpde[1:, 2] = ( (error.survey.md[1:] - error.survey.md[0:-1]) * np.amax(np.stack(( np.absolute( (error.survey.inc_rad[1:] - error.survey.inc_rad[:-1]) ), ( kwargs['tortuosity'] * (error.survey.md[1:] - error.survey.md[0:-1]) ) ), axis=-1), axis=-1) * -sin(error.survey.inc_rad[1:]) ) e_DIA = dpde * mag return error._generate_error( code, e_DIA, propagation, NEV, e_NEV=e_DIA, e_NEV_star=e_DIA )
[docs] def XYM3L(code, error, mag=0.0167, propagation='random', NEV=True, **kwargs): coeff = np.ones(len(error.survey.md) - 1) coeff = np.amax(np.stack(( coeff, sqrt( 10 / (error.survey.md[1:] - error.survey.md[:-1]) ) ), axis=-1), axis=-1) dpde = np.zeros((len(error.survey_rad), 3)) dpde[1:, 1] = np.absolute( cos(error.survey.inc_rad[1:]) * cos(error.survey.azi_true_rad[1:]) * coeff ) dpde[0, 1] = dpde[1, 1] with np.errstate(divide='ignore', invalid='ignore'): dpde[1:, 2] = np.nan_to_num( ( -np.absolute( cos(error.survey.inc_rad[1:]) ) * ( sin(error.survey.azi_true_rad[1:]) / sin(error.survey.inc_rad[1:]) ) * coeff ), posinf=0, neginf=0 ) dpde[0, 2] = dpde[1, 2] e_DIA = dpde * mag sing = np.where( error.survey_rad[:, 1] < error.survey.header.vertical_inc_limit ) if len(sing[0]) < 1: return error._generate_error(code, e_DIA, propagation, NEV) else: e_NEV = error._e_NEV(e_DIA) e_NEV_sing = np.zeros_like(e_NEV) e_NEV_sing[1:-1, 0] = ( coeff[:-1] * ( error.survey.md[2:] - error.survey.md[:-2] ) / 2 * mag ) e_NEV_sing[1, 0] = ( coeff[1] * ( error.survey.md[2] + error.survey.md[1] - 2 * error.survey.md[0] ) / 2 * mag ) e_NEV_sing[-1, 0] = ( coeff[-1] * ( error.survey.md[-1] - error.survey.md[-2] ) / 2 * mag ) e_NEV[sing] = e_NEV_sing[sing] e_NEV_star = error._e_NEV_star(e_DIA) e_NEV_star_sing = np.zeros_like(e_NEV) e_NEV_star_sing[1:, 0] = ( ( error.survey.md[1:] - error.survey.md[:-1] ) / 2 * mag ) e_NEV_star[sing] = e_NEV_star_sing[sing] return error._generate_error( code, e_DIA, propagation, NEV, e_NEV, e_NEV_star )
[docs] def XYM4L(code, error, mag=0.0167, propagation='random', NEV=True, **kwargs): propagation = 'random' coeff = np.ones(len(error.survey.md)) coeff[1:] = np.amax(np.stack(( coeff[1:], sqrt( 10 / (error.survey.md[1:] - error.survey.md[:-1]) ) ), axis=-1), axis=-1) dpde = np.zeros((len(error.survey_rad), 3)) with np.errstate(divide='ignore', invalid='ignore'): dpde[:, 2] = np.nan_to_num( np.absolute( cos(error.survey.inc_rad) * cos(error.survey.azi_true_rad) / sin(error.survey.inc_rad) * coeff ), posinf=0, neginf=0, ) dpde[:, 1] = ( np.absolute( cos(error.survey.inc_rad) ) * ( sin(error.survey.azi_true_rad) ) * coeff ) e_DIA = dpde * mag sing = np.where( error.survey_rad[:, 1] < error.survey.header.vertical_inc_limit ) if len(sing[0]) < 1: return error._generate_error(code, e_DIA, propagation, NEV) else: e_NEV = error._e_NEV(e_DIA) e_NEV_sing = np.zeros_like(e_NEV) e_NEV_sing[1:-1, 1] = ( coeff[1:-1] * ( error.survey.md[2:] - error.survey.md[:-2] ) / 2 * mag ) e_NEV_sing[1, 1] = ( coeff[1] * ( error.survey.md[2] + error.survey.md[1] - 2 * error.survey.md[0] ) / 2 * mag ) e_NEV_sing[-1, 1] = ( coeff[-1] * ( error.survey.md[-1] - error.survey.md[-2] ) / 2 * mag ) e_NEV[sing] = e_NEV_sing[sing] e_NEV_star = error._e_NEV_star(e_DIA) e_NEV_star_sing = np.zeros_like(e_NEV) e_NEV_star_sing[1:, 1] = ( ( error.survey.md[1:] - error.survey.md[:-1] ) / 2 * mag ) e_NEV_star_sing[1, 1] = ( ( error.survey.md[1] - error.survey.md[0] ) * mag ) e_NEV_star[sing] = e_NEV_star_sing[sing] return error._generate_error( code, e_DIA, propagation, NEV, e_NEV, e_NEV_star )