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

computeIntegrationPoints.py

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())