Skip to content

Instantly share code, notes, and snippets.

@bloyl
Last active November 30, 2016 17:11
Show Gist options
  • Save bloyl/84a2dc04d77ef0f300aef397f57201f2 to your computer and use it in GitHub Desktop.
Save bloyl/84a2dc04d77ef0f300aef397f57201f2 to your computer and use it in GitHub Desktop.
Compute integration points according to the last # formula in section 25.4.62 in the "Handbook of Mathematical Functions: # With Formulas, Graphs, and Mathematical Tables" edited by Abramowitz and Stegun. # http://people.math.sfu.ca/~cbm/aands/abramowitz_and_stegun.pdf
import numpy as np
def computeIntegrationPoints(h, acc, base=None, shape='circle'):
# These coil definitions make use of integration points according to the
# last formula in section 25.4.62 in the "Handbook of Mathematical
# Functions: With Formulas, Graphs, and Mathematical Tables"
# edited by Abramowitz and Stegun.
# http://people.math.sfu.ca/~cbm/aands/abramowitz_and_stegun.pdf
if (shape != 'circle'):
raise (Exception('Only circular coils are excepted'))
int_points = []
int_weights = []
if (acc == 0):
int_points.append([0, 0, 0])
int_weights.append(1)
elif (acc == 1):
int_points.append([h/2, h/2, 0])
int_weights.append(1/4.)
int_points.append([h/2, -h/2, 0])
int_weights.append(1/4.)
int_points.append([-h/2, h/2, 0])
int_weights.append(1/4.)
int_points.append([-h/2, -h/2, 0])
int_weights.append(1/4.)
elif (acc == 2):
int_points.append([0, 0, 0])
int_weights.append(1/4.)
x = np.sqrt(2/3.) * h
int_points.append([x, 0, 0])
int_weights.append(1/8.)
int_points.append([-x, 0, 0])
int_weights.append(1/8.)
x = np.sqrt(1/6.) * h
y = np.sqrt(2)/2. * h
int_points.append([x, y, 0])
int_weights.append(1/8.)
int_points.append([x, -y, 0])
int_weights.append(1/8.)
int_points.append([-x, y, 0])
int_weights.append(1/8.)
int_points.append([-x, -y, 0])
int_weights.append(1/8.)
if base and base > 0:
for w, p in zip(int_weights, int_points):
int_weights.append(-1 * w)
int_points.append([p[0], p[1], base])
return (int_weights, int_points)
def generate_coil_def(ident, senType, size, baseline, ori, name, coil_shape):
""" Generate coil_def lines for general Coil"""
# these are adapted from the top of coil_def.dat to match waht is actually
# in the file
hdr_format = '{} {:6d} {:3d} {:4d} {:10.3e} {:10.3e} "{}"\n'
pt_format = '{:8.4f} {:10.3e} {:10.3e} {:10.3e} {:6.3f} {:6.3f} {:6.3f}\n'
ret_string = ''
for acc in [0, 1, 2]:
ws, ps = computeIntegrationPoints(size/2., acc,
base=baseline, shape=coil_shape)
ret_string += hdr_format.format(senType, ident, acc, len(ws),
size, baseline, name)
for w, p in zip(ws, ps):
ret_string += pt_format.format(w, p[0], p[1], p[2],
ori[0], ori[1], ori[2])
return ret_string
def generate_artemis_coil_def():
"""Generate Coil_def lines for the artemis123 sensors.
Artemis123 has 3 sensor types (see Manual PDF)
1) circular Axial Gradiometers (meg sensors)
Coil diameter - 14.8590mm ( in)
Baseline - 57.4040mm ( in)
2) circular magnetometer (part of the 3 axis magnetometers)
Coil diameter - 14.8534mm ( in)
3) circular Reference Axial Gradiometers
Coil diameter - 14.8590mm ( in)
Baseline - 29.9974mm ( in)
Orientation for all these coils is set as the z axis. rotations are
handled in specifing the locs object in the measurement info.
"""
coil_def_str = ''
# orientation
ori = [0., 0., 1.]
coil_shape = 'circle'
# Axial grads first
ident = 7501
senType = 2
d = 14.8590E-3 # in meters
b = 57.4040E-3 # in meters
name = 'Artemis123 axial gradiometer ' + \
'size = {:0.3e} mm base = {:0.3e} mm'.format(d, b)
coil_def_str += generate_coil_def(ident, senType, d, b,
ori, name, coil_shape)
# Reference magnetometers
ident = 7502
senType = 1
d = 14.8534E-3 # in meters
b = 0 # in meters
name = 'Artemis123 reference magnetometer size = {:0.3e} mm'.format(d)
coil_def_str += generate_coil_def(ident, senType, d, b,
ori, name, coil_shape)
# Reference magnetometers
ident = 7503
senType = 2
d = 14.8590E-3 # in meters
b = 29.9974 # in meters
name = 'Artemis123 reference axial gradiometer ' + \
'size = {:0.3e} mm base = {:0.3e} mm'.format(d, b)
coil_def_str += generate_coil_def(ident, senType, d, b,
ori, name, coil_shape)
return coil_def_str
if __name__ == "__main__":
print (generate_artemis_coil_def())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment