''' A package with Python functions for forward and inverse calculations
of the EM-34 response. Written by Vincent Post, April 2007.'''
from numpy import array, sqrt
from scipy import optimize, hstack, r_, arange
def EMresponse(spacing, z, sigma, print_output = True, data = array([])):
'''Function that returns the theoretical response of the EM-34
assuming a horizontally-layered subsurface. Formulas from McNeill (1980).
Input:
spacing: array with loop spacings
z: array with depths to bottoms of layers
sigma: array with conductivities of the layers
print_output (optional): boolean to output results to screen
data (optional): array with the measured data to be printed
Output:
Array with 6 numbers. The first 3 numbers represent the theoretical
conductivity recorded by the EM-34 for loop spacings of 10, 20 and 40 m
for a vertical dipole (i.e., horizontal loops), the last 3 numbers
represent the theoretical conductivity recorded by the EM-34 for loop
spacings of 10, 20 and 40 m for a vertical dipole (i.e., horizontal
loops).'''
hloops = vloops = array([])
for s in spacing:
# For horizontal loops
z_s = z / s
z_contr = 1 / sqrt(1 + 4 * (z_s ** 2))
z_contr = hstack((1.0, z_contr, 0.0))
z_contr = z_contr[:-1] - z_contr[1:]
z_contr = z_contr * sigma
hloops = r_[hloops, sum(z_contr)]
#print 'hor. loops (vert. dipole)', s, 'm: ', sum(z_contr)
# For vertical loops
z_s = z / s
z_contr = sqrt(1 + 4 * (z_s ** 2)) - 2 * z_s
z_contr = hstack((1.0, z_contr, 0.0))
z_contr = z_contr[:-1] - z_contr[1:]
z_contr = z_contr * sigma
vloops = r_[vloops, sum(z_contr)]
#print 'vert. loops (hor. dipole)', s, 'm: ', sum(z_contr)
if print_output:
print_response(spacing, hloops, vloops)
if (len(data)):
print_data(spacing, data)
print_model(z, sigma)
return r_[hloops, vloops]
def collect_unknowns(z, sigma):
'''Function that looks for negative numbers in the z and sigma arrays. These are used as flags that indicate
they have to be optimized.
Input:
z: array with depths to bottoms of layers
sigma: array with conductivities of the layers
Output:
Array with the unknown depths followed by the unknown conductivities. These are passed on in the
optimize.leastsq function'''
return_value = array([])
for a in z:
if (a < 0):
return_value = r_[return_value, -a]
for a in sigma:
if (a < 0):
return_value = r_[return_value, -a]
return return_value
def EMresponse_wrapper(p, spacing, z, sigma):
'''This function is a wrapper that first updates the z and sigma arrays (or rather, local copies thereof) with the
numbers that are being optimized. It then calls EMresponse().
Input:
p: Array with the unknown depths followed by the unknown conductivities.
spacing: array with loop spacings
z: array with depths to bottoms of layers
sigma: array with conductivities of the layers
Output:
Same as EMresponse function. See above'''
z_tmp, sigma_tmp = update_model(p, z, sigma)
return EMresponse(spacing, z_tmp, sigma_tmp, False)
def update_model(p, z, sigma):
'''This function looks for negative numbers in the z and sigma arrays and inserts the numbers
that are being optimized (from the p array) at the appropriate positions.
Input:
p: Array with the unknown depths followed by the unknown conductivities.
z: array with depths to bottoms of layers
sigma: array with conductivities of the layers
Output:
Local copies of the updated z and sigma arrays'''
i = 0
# Note that we use local copies of z and sigma
z_tmp = z.copy()
for j in arange(len(z_tmp)):
if (z_tmp[j] < 0):
z_tmp[j] = p[i]
i = i + 1
sigma_tmp = sigma.copy()
for j in arange(len(sigma_tmp)):
if (sigma_tmp[j] < 0):
sigma_tmp[j] = p[i]
i = i + 1
return z_tmp, sigma_tmp
def fit_model(spacing, z, sigma, data):
'''Function that optimizes the layered conductivity model by minimizing the difference between the calculated
theoretical EM-34 response and the measured conductivities.
Input:
spacing: array with loop spacings
z: array with depths to bottoms of layers
sigma: array with conductivities of the layers
data: Array with 6 numbers. The first 3 numbers represent the EM-34 measured conductivity for loop spacings
of 10, 20 and 40 m for a vertical dipole (i.e., horizontal loops), the last 3 numbers represent the EM-34
measured conductivity for loop spacings of 10, 20 and 40 m for a vertical dipole (i.e., horizontal loops).
Output:
One array with the depths to bottom of the layers (optimized z) and one array with the conductivities (optimized
sigma) of the layers.'''
# First determine which dephts and conductivities need to be optimized and store in p0
p0 = collect_unknowns(z, sigma)
# Define the target function to which the field data will be fitted
fitfunc = lambda p, w, x, y: EMresponse_wrapper(p, w, x, y)
# Define the objective function
errfunc = lambda p, w, x, y, z: fitfunc(p, w, x, y) - z
# Perform the actual optimization using the optimize object
p1, success = optimize.leastsq(errfunc, p0, args = (spacing, z, sigma, data))
return update_model(p1, z, sigma)
def print_model(z, sigma):
'''Procedure thats prints the layered conductivity model as columns.
Input:
z: array with depths to bottoms of layers
sigma: array with conductivities of the layers'''
print ''
print 'Layered conductivity model'
print '=========================='
s = "%s\t%s\t%s" % ('Conduc', 'Resis', 'Depth')
s = "%s\n%s\t%s" % (s, 'tivity', 'tivity')
print s
print '---------------------'
i = 0
for v in sigma:
if (i < len(z)):
s = "%5.1f\t%7.1f\t%5.1f" % (v, 1000.0 / v, z[i])
else:
s = "%5.1f\t%7.1f" % (v, 1000.0 / v)
i = i + 1
print s
print ''
def print_response(spacing, hloops, vloops):
'''Procedure thats prints the calculated theoretical EM-34 response on the screen.
Input:
spacing: Array with loop spacing
hloops: Array with calculated EM-34 response for horizontal loops (vertical dipole)
vloops: Array with calculated EM-34 response for vertical loops (horizontal dipole)'''
print ''
print 'Calculated theoretical EM-34 response'
print '====================================='
s = "\t\t%s" % (' loop spacing (m) ')
print s
s = "\t\t%d\t%d\t%d" % (spacing[0], spacing[1], spacing[2])
print s
s = "\t\t%s" % ('---------------------')
print s
s = "%s\t%5.1f\t%5.1f\t%5.1f" % ('hor. loops', hloops[0], hloops[1], hloops[2])
s = "%s\n%s\t%5.1f\t%5.1f\t%5.1f" % (s, 'vert. loops', vloops[0], vloops[1], vloops[2])
print s
print ''
def print_data(spacing, data):
'''Procedure thats prints the calculated theoretical EM-34 response on the screen.
Input:
spacing: Array with loop spacing
hloops: Array with calculated EM-34 response for horizontal loops (vertical dipole)
vloops: Array with calculated EM-34 response for vertical loops (horizontal dipole)'''
print ''
print 'Measured EM-34 response'
print '====================================='
s = "\t\t%s" % (' loop spacing (m) ')
print s
s = "\t\t%d\t%d\t%d" % (spacing[0], spacing[1], spacing[2])
print s
s = "\t\t%s" % ('---------------------')
print s
s = "%s\t%5.1f\t%5.1f\t%5.1f" % ('hor. loops', data[0], data[1], data[2])
s = "%s\n%s\t%5.1f\t%5.1f\t%5.1f" % (s, 'vert. loops', data[3], data[4], data[5])
print s
print ''